# Forecasting the Contribution of Polarized Extragalactic Radio Sources in CMB
observations

G. Puglisi11affiliationmark:  22affiliationmark:  , V.
Galluzzi33affiliationmark:  55affiliationmark:  , L.
Bonavera44affiliationmark:  , J. Gonzalez-Nuevo44affiliationmark:  , A.
Lapi11affiliationmark:  , M. Massardi33affiliationmark:  66affiliationmark:  ,
F. Perrotta 11affiliationmark:  , C. Baccigalupi11affiliationmark:
22affiliationmark:  , A. Celotti11affiliationmark:  22affiliationmark:
77affiliationmark:  , L. Danese11affiliationmark:  1SISSA- International
School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy  2INFN-
National Institute for Nuclear Physics, Via Valerio 2, I-34127 Trieste, Italy
3INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, I-40129 Bologna,
Italy  4Departamento de Física, Universidad de Oviedo, C. Federico García
Lorca 18, 33007 Oviedo, Spain  5Dipartimento di Fisica e Astronomia,
Università di Bologna, via Gobetti 93/2, I-40129 Bologna, Italy  6Italian Alma
Regional Centre, Istituto di Radioastronomia, Via Piero Gobetti 101, I-40129
Bologna, Italy  7INAF, Osservatorio Astronomico di Brera, via Bianchi 46,
23807 Merate (LC<EMAIL_ADDRESS>
###### Abstract

We combine the latest datasets obtained with different surveys to study the
frequency dependence of polarized emission coming from Extragalactic Radio
Sources (ERS). We consider data over a very wide frequency range starting from
$1.4$ GHz up to $217$ GHz. This range is particularly interesting since it
overlaps the frequencies of the current and forthcoming Cosmic Microwave
Background (CMB) experiments. Current data suggest that at high radio
frequencies, ($\nu\geq 20$ GHz) the fractional polarization of ERS does not
depend on the total flux density. Conversely, recent datasets indicate a
moderate increase of polarization fraction as a function of frequency,
physically motivated by the fact that Faraday depolarization is expected to be
less relevant at high radio-frequencies. We compute ERS number counts using
updated models based on recent data, and we forecast the contribution of
unresolved ERS in CMB polarization spectra. Given the expected sensitivities
and the observational patch sizes of forthcoming CMB experiments about $\sim
200$ ( up to $\sim 2000$ ) polarized ERS are expected to be detected. Finally,
we assess that polarized ERS can contaminate the cosmological B-mode
polarization if the tensor-to-scalar ratio is $<0.05$ and they have to be
robustly controlled to _de-lens_ CMB B-modes at the arcminute angular scales.

###### Subject headings:

Cosmology: Cosmic Microwave Background – Radio Sources– observations

††thanks:  Corresponding authors: Giuseppe Puglisi<EMAIL_ADDRESS>
##  1\. Introduction

The Cosmic Microwave Background (CMB) is a relic radiation generated at the
decoupling of matter and radiation as the temperature of the Universe dropped
below ${3000}$ K. Its temperature and polarization anisotropies can be
exploited to probe the early stages of the Universe when an exponential
expansion, the so called _inflation_ might have occurred (Guth, 1981;
Starobinsky, 1982).

Since last decades, several experiments have tried to measure the CMB
polarized signal in order to find the imprints on its polarized anisotropies
of a stochastic background of primordial gravitational waves (PGW) that might
have been produced during the inflationary phase. Polarization anisotropies
are commonly decomposed into two scalar quantities called _E-_ and _B-modes_
(Seljak & Zaldarriaga, 1997; Hu & White, 1997), and to date, lots of efforts
have been made to observe the latter since their amplitude at degree scale is
expected to come mainly from PGW.

On one hand, E-mode photons get deflected via gravitational interaction by
intervening matter of large scale structures during the path toward us,
producing the so called _lensing B-modes_ at arcminute scale. Lensing B-modes
have been observed since four years (The Polarbear Collaboration: P. A. R. Ade
et al., 2014; Louis et al., 2017; Keisler et al., 2015; The Polarbear
Collaboration et al., 2017) with better and better accuracy and they represent
a powerful tool to probe the large scale structure of our Universe. On the
contrary, the primordial B-mode amplitude is unknown and is quantified by the
_tensor-to-scalar ratio_ , $r$, that relates the amplitude of tensor
perturbations of the space time metric, e.g. PGW, with respect to the scalar
perturbations. The joint collaboration of BICEP2 and Planck yielded so far the
latest upper limit on $r<0.07$ at $95\%$ confidence level (BICEP2/Keck and
Planck Collaborations et al., 2015). Meaning that the primordial B-mode
amplitude could be even lower than the lensing one.

To date, several challenges have prevented to detect primordial B-modes mostly
because of the diffuse polarized radiation coming from the Milky Way, known as
_Galactic Foregrounds_. The list of Galactic foregrounds is long and includes
anything emitting at sub-millimeter wavelengths between us and the CMB:
thermal dust, synchrotron radiation, free-free and several molecular line
emissions (Planck Collaboration et al., 2016b). All these emissions are
partially polarized and the main contribution comes from synchrotron and dust
(both polarized up to $20\%$ level Planck Collaboration et al., 2016e, d). At
high-frequency ($\nu>90$ GHz), such a large polarization degree is produced by
thermal dust grains aligning along the Galactic magnetic field lines. At low
frequencies ($\nu\lesssim 70$ GHz), cosmic electrons spiralling into the
Galactic magnetic field produce synchrotron radiation. Molecular lines are
expected to be polarized at lower levels $\lesssim 1\%$ (Goldreich & Kylafis,
1981; Puglisi et al., 2017), whereas free-free emission can be essentially
considered unpolarized. This is the justification of the recent efforts aimed
at observing the CMB polarization in a very wide range of frequencies and at
accurately characterizing both the spatial and frequency distribution of each
Galactic polarized foreground. Moreover, such an investigation allows to
design algorithms known as _component separation_ or _foreground cleaning_
techniques to extract B-modes out of a multi-frequency experimental setup.

For these reasons, (i) more focal plane pixels in multiple telescopes are
needed to increase sensitivity and (ii) multi-band polarization measurements
are required to recover the cosmic signal from the Galactic one via component
separation. As the focal plane will encode larger and larger number of
detectors, the next _stages_ in CMB experiment sensitivity will be achieved by
more accurately measuring $r$. To date, several ground based experiments are
updating their focal planes to a step forward from the so called CMB-Stage 2
(CMB-S2) to Stage 3 (CMB-S3 Arnold et al., 2014; Henderson et al., 2016;
Benson et al., 2014), including up to $10,000$ detectors observing up to $7\%$
of the sky. The ultimate step for a B-mode detection from the ground is
represented by CMBStage 4 experiments (CMB-S4 Abazajian et al., 2016), which
will account for up to $100,000$ detectors, observing half of the sky. CMB-S4
aims at measuring $r$ with the target accuracy $\sigma(r)\sim 0.0005$ .

At smaller scales the Extragalactic Radio Sources (ERS) and star-forming dusty
galaxies are the major contaminants (Tucci et al., 2011), although the latter
can also largely contribute to large angular scales due to clustering (De
Zotti et al., 2015). In this work, we mostly focus on the polarized emission
of ERS. To date, a few studies have been conducted regarding polarization of
ERS at the frequencies of CMB experiments (see Galluzzi & Massardi (2016) or
Bonavera et al. (2017a)) since polarization observations in the mm wavelength
bands are more challenging than at cm bands (at $1.4\div 20$ GHz) and
extrapolations are very common in this field of research (Tucci & Toffolatti,
2012).

The mechanism behind the polarized emission of radio sources is mostly due to
synchrotron radiation sourced by an Active Galactic Nucleus (AGN), where a
central super-massive black hole ($10^{6}\div 10^{9}M_{\odot}$) is hosted.
Most of the energy of an AGN comes from the gravitational potential energy of
the material located in a thin surrounding accretion disc, released as the
matter falls into the central black hole. Another component is constituted by
jets (usually paired) of material ejected toward the polar directions from the
black hole. Jets are observed to be very collimated and can travel very large
distances. Therefore, radio-galaxies seldom present double structures referred
as _lobes_ constantly fed by the jets of new energetic particles and magnetic
energy.

Depending on which components dominates the emission, such complex objects can
obviously appear with different morphologies and therefore be grouped in
different observational categories. One of the most important distinction is
related to the different orientations an AGN can be observed with respect to
the line of sight (see De Zotti et al. (2010) for a wide review). If edge-on,
the torus obscures the core and the inner disc, so that the emission is
dominated by the optically thin radio lobes presenting a steep spectral index
$\alpha$ at low frequencies $1\div 5$ GHz111The radio-source flux is described
by a power law $S_{\nu}\propto\nu^{-\alpha}$, and the threshold between flat
and steep spectral behaviour is commonly fixed at $\alpha=0.5$.. Objects with
$\alpha>0.5$ are commonly referred as _Steep Spectrum Radio Quasars_ (SSRQs)
and, generally, their optical counterpart is an elliptical galaxy. If seen
pole-on, the brightness is dominated by the approaching jet, the emission
looks compact and it is mostly Doppler boosted since particles move at
relativistic speeds. The emission is optically thick, does not contain many
optical features in the continuum but is characterized by a _flat_ spectrum
($\alpha<0.5$). Similar sources are called _Flat Spectrum Radio Quasars_
(FSRQs).

However, each source presents both the components, i.e. a flat-spectrum core
and extended steep-spectrum lobes, and it can be easily understood that a
simple-power law cannot be applied to resemble the large radio frequency range
(Massardi et al., 2011; Bonaldi et al., 2013). External and self-absorption,
from free-free and synchrotron, may affect and change the dependence of
$S_{\nu}$, so that the spectrum could increase as a function of frequency
(Galluzzi et al., 2017).

There is an increasing interest on polarization of ERS at high-radio
frequencies not only to better understand the physics behind the emitting
system, e.g. the degree of ordering of the magnetic field, the direction of
its field lines (Tucci et al., 2011), but also because polarized ERS will be
largely detected by forthcoming CMB experiments. Furthermore, the ERS
contaminating signal in the polarization power spectra cannot be neglected to
assess the power spectrum of lensing B-modes. This is the reason why recent
works in the literature can be found addressing this issue: De Zotti et al.
(2015, 2016) predicted the contribution in polarization both for ERS and dusty
galaxies at frequency channels of the Cosmic ORigin Explorer (CORE) satellite;
Curto et al. (2013) estimated for future CMB missions the contamination
produced by radio and far-Infrared sources at the level of bispectrum
considering different shapes of the primordial non-Gaussianity parameter,
$f_{nl}$.

In section 2 we describe the datasets we combine in order to determine the
polarization dependence as a function of frequency, discussed in section 4. In
section 3, we present the models for number counts adopted in this analysis.
In section 5, we show the results of a forecast package we developed to assess
the contamination of polarized ERS in terms of CMB power spectra given the
nominal specifics of current and forthcoming CMB experiments. Finally, we
devote section 6 to discuss and summarize our results.

##  2\. Data

In this section, we present the data collected from publicly available
catalogues. The data, summarized in Table 1, have been used to characterize
the polarization fraction of ERS in about two orders of magnitude in frequency
range (i.e. from $1.4$ to $217$ GHz).

Table 1Summary of the catalogues that we use Sec.5.  | Frequency [GHz] | Sky Region | FWHM | Detect. flux |  $90\%$ Compl. |  $\\#$ Sources
---|---|---|---|---|---|---
NVSS | 1.4 | $\delta>-40^{\circ}$ | $45^{\prime\prime}$ |  $0.29$ mJy$/\mathrm{beam}$ |  $2.3$ mJy | $1.8\times 10^{6}$
S-PASS | 2.3 | $\delta<-1^{\circ}$ | $8.9^{\prime}$ |  $1$ mJy$/\mathrm{beam}$ |  $420$ mJy | $533$
JVAS | 8.4 | $\delta\geq 0^{\circ},\,|b|\geq 2.5^{\circ}$ | $0.2{{}^{\prime\prime}}$ |  $50$ mJy |  $200$ mJy | 2720
CLASS | 8.4 | $0\geq\delta\geq 70^{\circ}$ | $0.2{{}^{\prime\prime}}$ |  $20$ mJy |  $30$ mJy | 16503
AT20G |  $4.8$, $8.6$, $20$ | $\delta<0^{\circ},\,|b|<1.5^{\circ}$ | $10{{}^{\prime\prime}},6{{}^{\prime\prime}},11{{}^{\prime\prime}}$ |  $40$ mJy |  $100$ mJy$/\mathrm{beam}$ | 5890
VLA |  $4.8$, $8.5$, | $\delta>-15^{\circ}$ |  $12^{\prime\prime}$, $6^{\prime\prime}$, |  $0.7$, $0.3$, |  $40$ mJy | $159$
$22.5$, $43.5$ |  |  $4^{\prime\prime}$,$2^{\prime\prime}$ |  $0.9$, $1.2$ mJy$/\mathrm{beam}$ |  |
PACO | $20$ |  Ecl. lat.$<-65^{\circ}$ | $11{{}^{\prime\prime}}$ |  $40$ mJy |  $200$ mJy | $104$
XPOL-IRAM | $86$ | $\delta>30^{\circ}$ | $28^{\prime\prime}$ |  $0.5$ Jy |  $1$ Jy | 145
PCCS2 |  $30$, $44$, |  |  $32.4^{\prime}$, $27.1^{\prime}$, |  $117$,$229$, |  $427$,$692$, |  $1560$,$934$,
$70$, $100$, | Full sky |  $13.3^{\prime}$, $9.7^{\prime}$, |  $225$, $106$, |  $501$,$269$, |  $1296$,$1742$,
$143$, $217$ |  |  $7.3^{\prime}$, $5.0^{\prime}$ |  $75$,$81$ mJy |  $177$,$152$ mJy |  $2160$,$2135$

###  2.1. The S-PASS/NVSS joint catalogue

The S-band Polarization All Sky Survey (S-PASS) survey observed the Southern
sky with declination $\delta<-1^{\circ}$ at 2.3 GHz with full width at half
maximum (FWHM) of 8.9 arcmin both in total intensity and polarization using
the $64$ m Parkes Radio Telescope. Lamee et al. (2016) cross-matched it with
the NRAO/VLA Sky Survey, Condon et al. (NVSS 1998), at $1.4$GHz (45 arcsec
(FWHM) and rms total brightness fluctuations of $\sim
0.29\,\mathrm{mJy\,beam^{-1}}$). Lamee et al. (2016) aimed at generating a
novel and independent polarization
catalogue222http://vizier.cfa.harvard.edu/viz-bin/VizieR?-source=J/ApJ/829/5
enclosing 533 bright ERS at $2.3$ GHz with polarized flux-density stronger
than $420\,\mathrm{mJy}$.

###  2.2. The JVAS/CLASS 8.4 GHz catalogue

We used the data from the JVAS/CLASS 8.4-GHz catalogue Jackson et al.
(2007)333http://vizier.cfa.harvard.edu/viz-bin/VizieR?-source=J/MNRAS/376/371,
which combined data taken from the Jodrell-VLA Astrometric Survey (JVAS) and
the Cosmic Lens All-Sky Survey (CLASS) both observing at $8.4$ GHz. The former
detected 2720 sources stronger than $200$ mJy in total intensity at $5$ GHz
and $\delta\geq 0^{\circ}$, masking the Galactic mid-plane at Galactic
latitude $|b|\geq 2.5^{\circ}$. To complement JVAS, CLASS consisted of all
sources with a fainter $5$ GHz flux, i.e. $S>30$ mJy observed in a sky region
between $0^{\circ}\leq\delta\leq 70^{\circ}$. Combining the two surveys, a
sample of $16\,503$ FSRQ intensity fluxes has been collected.

Jackson et al. (2010) were able to assess polarized fluxes only for a few
objects from the 133 sources observed by WMAP at $22$ and $43$ GHz Wright et
al. ($S>1$ Jy 2009) with counterpart in the JVAS/CLASS catalogues. For the
purposes of our work this sample was not large enough to be included in the
following analysis.

However, we exploit the data selection described by Pelgrims & Hutsemékers
(2015) that considered all the sources with polarized flux $\geq 1$ mJy in
order to obtain an unbiased sample of 3858 NED identified sources. We selected
2829 sources classified by Pelgrims & Hutsemékers (2015) as _QSOs_ and _Radio
Sources_. For a complete description of the catalogue and the surveys refer to
Jackson et al. (2007).

###  2.3. The AT20G Survey

The Australia Telescope 20 GHz (AT20G) Survey observed blindly the Southern
sky ($\delta<0^{\circ}$ excluding the Galactic plane strip at
$|b|<{1.5}^{\circ}$) at 20 GHz with the Australia Telescope Compact Array
(ATCA) from 2004 to 2009, (Murphy et al., 2010). The detected sources were
followed up almost simultaneously at 4.8 and 8.6 GHz. The AT20G source
catalogue444http://vizier.cfa.harvard.edu/viz-
bin/VizieR?-source=J/MNRAS/402/2403 includes 5890 sources at 20 GHz above the
total intensity detection limit of $40$ mJy, of which 3332 were detected at
all the observing frequencies. Averaged on the whole area of the survey, the
catalogue is $91\%$ complete above $100\,\mathrm{mJy\,beam^{-1}}$ (Murphy et
al., 2010). Polarization of sources was considered detected if the following
criteria were satisfied: polarized flux density $P>6$ mJy or at least three
times larger than its rms error, and polarized fraction above $1$ per cent.
Massardi et al. (2011) presented an analysis to characterize the radio
spectral properties of the whole sample both in total intensity and
polarization, involving $768$ sources detected at $20$ GHz ($467$ of them were
also detected in polarization at $4.8$ and/or at $8.6$ GHz). Given the goal of
this work, we include polarized flux densities from $3332$ sources, $2444$ of
them presenting a flat spectrum in total intensity, $|\alpha_{5}^{8}|<0.5$ and
the remaining $888$ a steep-spectrum sources ($|\alpha_{5}^{8}|>0.5$).

###  2.4. The VLA observations

Sajina et al. (2011) presented measurements555http://vizier.u-strasbg.fr/viz-
bin/VizieR?-source=J/ApJ/732/45 in flux densities and polarization of 159 ERSs
detected with the Very Large Array (VLA) at four frequency channels,
$4.86,\,8.46,\,22.46,\,43.34$ GHz. This sample was selected from the AT20G one
(Murphy et al., 2010; Massardi et al., 2011) by requiring a flux density
$S>40$ mJy in the equatorial field of the Atacama Cosmology Telescope (ACT)
survey on a region at declination north of $-15^{\circ}$ and excluding the
Galactic plane. The aim of this program was firstly to characterize the
spectra and variability both in total intensity and polarization of high-
frequency-selected radio sources and to improve the estimation of the ERS
contamination at high-frequency for CMB experiments.

In $40\%$ of the whole sample they detected polarized flux density in all the
bands, and observed an increasing trend of polarization fraction as a function
of frequency, more evident for SSRQs.

###  2.5. PACO with ATCA and ALMA

The Planck-ATCA Coeval Observations (PACO) project detected 464 sources
selected from the AT20G catalogue during 65 epochs between July 2009 and
August 2010, at frequencies ranging from 5.5 to 39 GHz with the ATCA. The
sources were simultaneously observed (within 10 days) by the Planck satellite
(Massardi et al., 2011; Bonavera et al., 2011). The project aimed at
characterizing, together with Planck data, the variability and spectral
behavior of sources over a wide frequency range (up to 857 GHz for some
sources), in total intensity only. The catalogue includes a complete sample of
159 sources selected to be brighter than $200$ mJy at $\delta<30^{\circ}$
(excluding the Galactic midplane $|b|<5^{\circ}$). A sub-sample of 104 of
these sources with ecliptic latitude $<-65^{\circ}$ (which coincides to one of
the _deep_ patches most frequently scanned by the Planck satellite scanning
strategy) has been re-observed with high sensitivity in polarization with ATCA
in 2014 and 2016 in the $1.1-39$ GHz frequency range (Galluzzi et al., 2017).
32 of them have been also followed up at 95 GHz onto 3 circular regions
($10^{\circ}$ of diameter) at ecliptic latitude $<-75^{\circ}$ with the
Atacama Large Millimeter Array (ALMA) to better characterize the polarization
properties of ERS at the frequencies of many CMB experiments and allowing an
accurate study of few reference targets which could be exploited for
calibration and validation of cosmological results. Further details will be
described in a companion paper (Galluzzi et al. 2018, in prep.). Data from
both 20 and 95 GHz have been included in this analysis.

###  2.6. First 3.5 mm Polarimetric Survey

Agudo et al. (2010) presented for the first time polarimetric data at $86$ GHz
of a sample of 145 flat spectrum radio galaxies at different epochs (from 2005
July to 2009 October)666http://vizier.u-strasbg.fr/viz-
bin/VizieR?-source=J/ApJS/189/1.. The measurements have been performed by
means of the XPOL polarimeter of the IRAM 30 m telescope, by selecting the
sources observed from 1978 to 1994 at $\delta>30^{\circ}$ whose total
intensity was above $\gtrsim 1$ Jy. They detected above $>3\sigma$ level
$1.5\%$ linear and $0.3\%$ circular polarization degree respectively for
$76\%$ and $6\%$ of the whole sample. Remarkably, they found a factor of $\sim
2$ excess in the polarization fraction at 86 GHz with respect to that measured
at 15 GHz.

###  2.7. The Second Planck Catalogue of Compact Sources

We exploit data from latest Planck Catalogue of Compact Sources (PCCS2 Planck
Collaboration, 2015)777http://pla.esac.esa.int/pla/ including polarimetric
detection of sources between $30$ and $353$ GHz from August 2009 to August
2013. The total intensity $90\%$ completeness ranges from $177$ to $692$ mJy
in this regime of frequencies, allowing to detect thousands of sources matched
both internally (between neighbor Planck channels) and with external
catalogues. On the contrary, the instrumental noise in polarization and the
presence of polarized Galactic foregrounds limited the number of polarized
sources to a few tens (with the exception of the $30$ GHz channel where 113
polarized sources were detected).

It is straightforward to state that only sources with high fractional
polarization have been detected by Planck and thus the statistics of ERS
polarization can be biased upward. Bonavera et al. (2017a) recently proposed a
methodology to cope with this issue by means of applying a _stacking
technique_ to Planck data. They used as main sample the 30 GHz catalogue,
consisting of 1560 sources above $S>427$ mJy at $90\%$ completeness level and
then followed the sample at higher Planck frequency maps. They further
distinguished sources inside and outside the Galactic plane defined by the
Planck Galactic mask GAL060 $(f_{sky}\approx 60\%)$ and the exclusion of the
Small and Large Magellanic clouds. This technique has been already applied by
Stil et al. (2014) to NVSS dataset to study the faint polarized signal of ERS
detected in total intensity: the signal from many weak sources is co-added to
achieve a statistical detection. Bonavera et al. (2017a) found that the ERS
polarization fraction is approximately constant with frequency over the Planck
frequency range. An alternative approach that attempts to overcome some of the
intrinsic statistical limitations of the stacking technique have been recently
exploited by Trombetti et al. (2017) obtaining results comparable both with
Bonavera et al. (2017a, b) and with other ground based observations.

We used both data coming from the PCCS2 catalogue and from Bonavera et al.
(2017a).

##  3\. Model for Number Counts

Figure 1.— Euclidean differential number counts at (top) 20 and (bottom) 95
GHz. Thick dotted, dashed, dot-dashed and solid lines are respectively the
number counts of BL Lacs, FSRQs, SSRQs and their total contribution predicted
by the D05 model (de Zotti et al., 2005). The thick solid gray line show the
number counts prediction from the C2Ex model (Tucci et al., 2011). Thinner
lines follow the same color scheme as the thick ones and refer to polarization
number counts, computed via a convolution with a log-normal distribution
function fitted from the data. Number counts estimates from several surveys
are also shown. (top) The circle data point in the upper curves are data from
AT20G (Massardi et al., 2008), whereas upper triangles are from WMAP5-yr
survey (K-band Massardi et al., 2009); in lower curves polarization number
counts from a resampling of PACO data (Galluzzi et al., 2018, , circles) and
from WMAP polarization point source catalogue (Lopez-Caniego et al., 2009, ,
upper triangles). (bottom panel) Diamonds in upper curves are number counts
from SPT (Mocanu et al., 2013), squares are from Planck ERCSC catalogue
(Planck Collaboration et al., 2011a); the lower triangles have been obtained
from a bootstrap resampling of 32 polarized fluxes detected with PACO at 95
GHz.

We adopted the evolutionary model proposed by de Zotti et al. (2005,
hereafter, D05) that describes the population properties of ERSs and dusty
galaxies above $\nu\gtrsim 5$ GHz. The model assumes a simple analytic
luminosity evolution in order to fit the available data on local luminosity
functions (LF), source counts888Available online
http://w1.ira.inaf.it/rstools/srccnt/srccnt_tables.html.  and redshift
distributions for sources down to few mJy. It determines the epoch-dependent
LF starting from local LFs for several source populations. For each population
the model adopts a different evolution laws estimating a set of free
parameters from available data. Recently, Bonato et al. (2017) and Mancuso et
al. (2017) improved the predictions of D05 model by updating the LF and
redshift evolution with state of art data of radio-emitting star-forming
galaxies and AGNs.

The D05 model assumes a power-law spectrum for each considered population of
ERS and each one is described by one (or at most two) constant spectral index.
This simple assumptions could not hold anymore when large frequency ranges are
taken into account. Departures from single power law spectra are expected
because of (i) electron ageing (ii) transition from optically thick to
optically thin regime, (iii) different components yielding different spectral
contributions at different frequencies. Therefore, this _simplified_ model
requires adjusting when source counts measurement are observed at frequencies
$>40$ GHz.

Tucci et al. (2011) showed that radio spectra in AGN cores can differ from a
single power-law when large frequency intervals are considered. In particular,
they focused on the blazar spectra for which a steepening of the spectral
index from $0.5$ to $1.2$ has been observed (Planck Collaboration et al.,
2011a, b) due to the transition from optically-thick to optically-thin
synchrotron emission of AGN jets (Kellermann, 1966; Blandford & Koenigl,
1979). Therefore, Tucci et al. (2011) proposed the so called _C2Ex_ model that
assumes a spectral break and different parameters for BL Lacs and FSRQs and
allows to properly fit the number counts especially at high-frequency
($\nu\geq 100$ GHz). Furthermore, Planck Collaboration (2015, XXVI) found that
all radio sources observed at the Low Frequency Instrument (LFI) channels
present flat and narrow spectral index distribution with $\alpha_{LFI}\lesssim
0.2$, whereas sources in the High Frequency Instrument (HFI) catalogues have a
broader distribution showing a steeper spectral index, $\alpha_{HFI}\gtrsim
0.5$ and these findings supports the scenario of BL Lac transition happening
at larger frequencies $\nu>100$ GHz with respect to the FSRQ one (at
$10<\nu<100$ GHz).

We plot in Figure 1 as thicker curves. the differential number counts, $n(S)$,
predicted with D05 and C2Ex models respectively as blue and grey thick solid
lines. The top (bottom) panel refers to number counts at $20\,(95)$
GHz999Source number counts for a wider range of frequencies are shown in
Figure A.1.. We also plot the contributions estimated by the D05 model for BL
Lacs, FSRQs, SSRQs respectively as dotted, dashed, dot-dashed lines. To
compare the quantities with those expected in a Euclidean Universe, counts are
normalized by a factor of $S^{5/2}$. The data points shown are number counts
as measured by AT20G survey (Massardi et al., 2008, blue circles), from South
Pole Telescope (SPT Vieira et al., 2010; Mocanu et al., 2013, blue diamonds),
from the Wilkinson Microwave Anisotropy Probe (WMAP Massardi et al., 2009,
yellow upper triangles) and from Planck (Planck Collaboration et al., 2011a,
2013, yellow squares).

The lower thinner curves in Figure 1 are Euclidean normalized differential
polarized emission number counts, $P^{5/2}n(P)$, computed from polarized flux-
density measurements and will be discussed in Sec. 4.

By comparing the predictions from the two models, we find that both are in a
reasonable agreement, with differences well below the uncertainties at $20$
GHz. However, as discussed above and shown in the bottom panel of Figure 1,
number counts estimated with D05 are systematically a factor of $\sim 2$
higher than the C2Ex ones at larger fluxes $100$ mJy, consistently with the
findings of Planck Collaboration et al. (2011a).

In the following, we make use of both D05 and C2Ex models to assess
respectively _conservative_ and _realistic_ estimates of polarized ERS to CMB
measurements.

##  4\. Statistical properties of ERS polarization fraction

Figure 2.—  Distribution function of polarization fraction for data at 20 GHz
(left) and at 95 GHz (right). The best fit values of log-normal parameters are
shown. The reduced $\tilde{\chi}^{2}$ estimated from the fit is 0.13 and 0.15
respectively for left and right panels.

Polarization number counts have to be assessed to know how many sources can be
detected at a certain polarized flux density, $P=\sqrt{Q^{2}+U^{2}}$, with $Q$
and $U$ being the linear polarization _Stokes_ parameters. Polarization
measurements at mm-wavelengths are scarce because of the faintness of the
polarized signal, so that both high sensitivity and robust estimates of
systematic effects are required. Furthermore, completeness is very hard to be
achieved with polarized samples. This is the reason why, to date,
extrapolations from low frequency observations ($1.4\div 5$ GHz) are commonly
adopted though the uncertainties due to intra-beam effects and bandwidth
depolarization may seriously affect the estimation.

To encompass this issue, several works in the literature (Battye et al., 2011;
Tucci & Toffolatti, 2012; Massardi et al., 2013; Bonavera et al., 2017a)
considered the probability function $\mathcal{P}(\Pi)$ of the polarization
fraction, $\Pi=P/S$. Differential polarization number counts can be defined as

$\displaystyle n(P)=$ $\displaystyle
N\int_{S_{0}=P}^{\infty}\mathcal{P}(P,S)dS=N\int_{S_{0}=P}^{\infty}\mathcal{P}(\Pi,S)\frac{dS}{S}$
$\displaystyle=$
$\displaystyle\int_{S_{0}=P}^{\infty}\mathcal{P}(\Pi)n(S)\frac{dS}{S},$ (1)

where $N$ is the total number of sources with $S\geq S_{0}$,
$\mathcal{P}(P,S)$ and $\mathcal{P}(\Pi,S)$ are the probability functions of
finding a source with flux $S$ and polarized flux $P$ or polarization fraction
$\Pi$ and both can be constrained from observations.

Notice that, in the last equation of (1), we assume that $\Pi$ and $S$ are
statistically independent. On one hand, recent results at low frequencies
indicate that this might not be the case: Stil et al. (2014) found that
fainter sources ($\sim$ 1 mJy) of NVSS catalogue present a higher median
fractional polarization. These results were confirmed by Lamee et al. (2016)
with S-PASS: they found indications of a possible correlation between the
polarization fraction and total intensity of steep-spectrum sources ranging
from $0.42$ to $10$ Jy, whereas the correlation disappears when FSRQs are
involved. On the other hand, at higher frequencies (above $20$ GHz), Massardi
et al. (2008) and Tucci & Toffolatti (2012) did not find a clear correlation
between $\Pi$ and $S$ (at fluxes above $500$ mJy) for both FSRQs and SSRQs,
but they found fractional polarization correlating at frequencies between
$4.8$ and $20$ GHz.

To date, surveys at high-frequencies have not been sensitive enough to probe
fainter polarized fluxes in order to seek whether this assumption holds or
not. Tucci et al. (2004) further argued that at higher frequencies we observe
two possible effects: (i) depolarization from Faraday rotation is essentially
negligible at frequencies above $\nu\gtrsim 10$ GHz, (ii) by observing compact
objects (i.e. FSRQs) at increasing frequency, we probe inner and inner
regions, closer to the nucleus where the magnetic field is expected to be
highly ordered. Consequently if this is the case, _the polarization fraction
may increase with the frequency_.

Given the goals of our work and the fact that frequencies above $10$ GHz are
involved in the forecast analysis, we assume polarized fraction and flux-
density uncorrelated and statistically independent but we look for some
eventual dependence of $\Pi$ as function of frequency.

Following Battye et al. (2011), we model $\mathcal{P}(\Pi)$ by means of a log-
normal distribution, i.e.

$\mathcal{P}\left(\Pi\right)=\frac{A}{\sqrt{2\pi\sigma^{2}}\Pi}\exp\left[-\frac{\left(\ln(\Pi/\mu)\right)^{2}}{2\sigma^{2}}\right]$
(2)

where $\mu$ and $\sigma$ are respectively the median and the standard
deviation in log. Notice that eq. (2) holds only if $0\leq\Pi<\infty$.
Although an infinite value of $\Pi$ does not have any physical meaning
(synchrotron emission can be polarized up to $75\%$), the values of $\mu$ and
$\sigma$ are orders of magnitude smaller. Thus $\Pi$ can be effectively
assumed to range up to a large value. This allows us to write a good
approximation of the fractional polarization by a combination of the log-
normal parameters101010For further details refer to Battye et al. (2011).

$\displaystyle\langle\Pi\rangle$ $\displaystyle\approx$ $\displaystyle\mu
e^{\frac{1}{2}\sigma^{2}},$ (3) $\displaystyle\langle\Pi^{2}\rangle$
$\displaystyle\approx$ $\displaystyle\mu^{2}e^{2\sigma^{2}},$ (4)
$\displaystyle\Pi_{med}$ $\displaystyle\approx$ $\displaystyle\mu.$ (5)

We derive polarization fraction distribution by using a _bootstrap-resampling_
method outlined in Austermann et al. (2009). This generates $N_{resamp}$
simulations of the catalogue and values for unpolarized and polarized flux
densities are randomly assigned for each source, from a normal distribution
$\mathcal{N}(\mu_{src},\sigma_{src})$ peaking at the observed value
$\mu_{src}$ and with a width $\sigma_{src}$ equal to the flux uncertainty. In
the case of upper limits, a random number is extracted from a normal
distribution centred on 0 and with width $\sigma_{src}$. For each resampling
we compute the polarization fraction and the values are distributed across
bins (ranging from 5 to 15 bins depending on the number of data collected in
each catalogue). The final distribution is thus given by the mean value within
each bin and vertical error bars computed by means of Poisson statistics, at
$68\%$ of confidence level (CL, Gehrels, 1986), counting the observed sources
in each polarization fraction bin. Finally, a log-normal distribution function
(2) is fitted from each dataset and $\langle\Pi\rangle$,
$\langle\Pi^{2}\rangle$ and $\Pi_{med}$ are then estimated from the log-normal
parameters $\mu$ and $\sigma$ as in (3),(4),(5).

Table 2Values of log-normal parameters obtained by fitting data from each
catalogue.  Flat-spectrum sources
---
$\nu$[GHz] | $N_{src}$ | $A$ | $\mu$ | $\sigma$ | $\langle\Pi\rangle$ | $\Pi_{med}$ | $\langle\Pi^{2}\rangle^{1/2}$ | Reference
1.4 | 82 | $0.54\pm 0.08$ | $1.73\pm 0.24$ | $1.05\pm 0.09$ | $2.98\pm 0.64$ | $1.72\pm 0.24$ | $5.15\pm 1.53$ | Lamee et al. (2016)
2.3 | 82 | $0.53\pm 0.07$ | $1.51\pm 0.13$ | $1.05\pm 0.07$ | $2.64\pm 0.36$ | $1.52\pm 0.13$ | $4.59\pm 0.87$ | Lamee et al. (2016)
4.8 | 2335 | $1.57\pm 0.07$ | $2.36\pm 0.02$ | $0.75\pm 0.01$ | $3.14\pm 0.03$ | $2.37\pm 0.02$ | $4.16\pm 0.08$ | Murphy et al. (2010)
8.6 | 2335 | $1.55\pm 0.02$ | $2.46\pm 0.01$ | $0.73\pm 0.01$ | $3.21\pm 0.03$ | $2.46\pm 0.01$ | $4.20\pm 0.06$ | Murphy et al. (2010)
8.6 | 2827 | $0.52\pm 0.01$ | $2.41\pm 0.05$ | $0.76\pm 0.01$ | $3.23\pm 0.08$ | $2.41\pm 0.05$ | $4.31\pm 0.14$ | Pelgrims & Hutsemékers (2015)
4.8 | 109 | $0.60\pm 0.06$ | $2.02\pm 0.13$ | $0.84\pm 0.05$ | $2.89\pm 0.28$ | $2.02\pm 0.13$ | $4.12\pm 0.55$ | Sajina et al. (2011)
8.6 | 109 | $0.74\pm 0.14$ | $2.12\pm 0.24$ | $0.84\pm 0.09$ | $3.01\pm 0.51$ | $2.12\pm 0.23$ | $4.27\pm 1.02$ | Sajina et al. (2011)
22 | 155 | $1.36\pm 0.09$ | $3.1\pm 0.10$ | $0.88\pm 0.03$ | $4.57\pm 0.22$ | $3.10\pm 0.09$ | $6.74\pm 0.46$ | Sajina et al. (2011)
43 | 111 | $2.59\pm 0.08$ | $4.48\pm 0.06$ | $1.00\pm 0.03$ | $7.42\pm 0.17$ | $4.47\pm 0.06$ | $12.32\pm 0.41$ | Sajina et al. (2011)
20 | 104 | $0.73\pm 0.05$ | $1.73\pm 0.16$ | $0.98\pm 0.06$ | $2.91\pm 0.42$ | $1.73\pm 0.16$ | $4.89\pm 0.99$ | Galluzzi et al. (2018)
89 | 145 | $1.20\pm 0.06$ | $2.86\pm 0.10$ | $0.64\pm 0.03$ | $3.52\pm 0.17$ | $2.86\pm 0.10$ | $4.32\pm 0.28$ | Agudo et al. (2010)
95 | 32 | $1.09\pm 0.21$ | $2.13\pm 0.23$ | $0.97\pm 0.09$ | $3.20\pm 0.60$ | $2.07\pm 0.24$ | $4.94\pm 1.32$ | This work
30 | 114 | $1.51\pm 0.23$ | $2.05\pm 0.36$ | $1.08\pm 0.08$ | $3.69\pm 0.92$ | $2.06\pm 0.37$ | $6.61\pm 2.19$ | Planck Collaboration (2015)
44 | 30 | $2.63\pm 0.26$ | $2.72\pm 0.26$ | $0.77\pm 0.11$ | $3.69\pm 0.66$ | $2.73\pm 0.26$ | $5.00\pm 1.32$ | Planck Collaboration (2015)
70 | 34 | $3.91\pm 0.55$ | $2.52\pm 0.05$ | $0.58\pm 0.06$ | $2.97\pm 0.15$ | $2.51\pm 0.05$ | $3.52\pm 0.30$ | Planck Collaboration (2015)
100 | 20 | $2.18\pm 0.28$ | $5.15\pm 0.69$ | $0.80\pm 0.10$ | $7.19\pm 1.59$ | $5.17\pm 0.73$ | $9.99\pm 3.07$ | Planck Collaboration (2015)
143 | 25 | $3.13\pm 0.10$ | $5.98\pm 0.16$ | $0.80\pm 0.04$ | $8.39\pm 0.35$ | $6.02\pm 0.13$ | $11.69\pm 0.80$ | Planck Collaboration (2015)
217 | 11 | $3.44\pm 0.32$ | $3.74\pm 0.29$ | $0.88\pm 0.11$ | $5.47\pm 0.34$ | $3.70\pm 0.27$ | $8.09\pm 1.09$ | Planck Collaboration (2015)
Steep-spectrum sources
1.4 | 388 | $1.12\pm 0.08$ | $1.47\pm 0.11$ | $1.05\pm 0.08$ | $2.56\pm 0.35$ | $1.47\pm 0.11$ | $4.45\pm 0.95$ | Lamee et al. (2016)
2.3 | 388 | $1.78\pm 0.07$ | $1.93\pm 0.06$ | $0.80\pm 0.05$ | $2.66\pm 0.14$ | $1.93\pm 0.06$ | $3.66\pm 0.31$ | Lamee et al. (2016)
4.8 | 952 | $2.07\pm 0.07$ | $2.83\pm 0.07$ | $0.81\pm 0.04$ | $3.92\pm 0.15$ | $2.84\pm 0.08$ | $5.43\pm 0.37$ | Murphy et al. (2010)
8.4 | 952 | $3.02\pm 0.03$ | $2.13\pm 0.12$ | $1.13\pm 0.03$ | $4.85\pm 0.05$ | $3.02\pm 0.05$ | $7.79\pm 0.18$ | Murphy et al. (2010)
20 | 952 | $4.55\pm 0.12$ | $6.98\pm 0.12$ | $0.55\pm 0.01$ | $8.10\pm 0.18$ | $6.98\pm 0.12$ | $9.41\pm 0.27$ | Murphy et al. (2010)
4.8 | 39 | $2.72\pm 0.65$ | $2.35\pm 0.46$ | $1.07\pm 0.43$ | $4.19\pm 1.42$ | $2.35\pm 0.42$ | $7.49\pm 5.61$ | Sajina et al. (2011)
8.6 | 39 | $1.94\pm 0.14$ | $3.39\pm 0.31$ | $1.04\pm 0.10$ | $5.90\pm 1.10$ | $3.41\pm 0.32$ | $10.23\pm 2.95$ | Sajina et al. (2011)
22 | 38 | $2.51\pm 0.10$ | $5.76\pm 0.19$ | $0.82\pm 0.05$ | $8.08\pm 0.44$ | $5.73\pm 0.17$ | $11.40\pm 0.99$ | Sajina et al. (2011)
43 | 15 | $4.74\pm 0.08$ | $9.89\pm 0.20$ | $0.73\pm 0.02$ | $12.43\pm 0.29$ | $9.62\pm 0.13$ | $16.06\pm 0.55$ | Sajina et al. (2011)

We show in Figure 2 the polarization fraction distributions from PACO-ATCA at
20 GHz and PACO-ALMA at 95 GHz (the best fit parameters of the other datasets
used in this analysis are summarized in Table 2). We show in top panel of
Figure 1 the polarization number counts computed by Galluzzi et al. (2018) at
$20$ GHz (blue circles) as a result of the convolution of total intensity
number counts with the log-normal distribution $\mathcal{P}(\Pi)$ as in eq.
(1). We further overlap the predicted total counts from both the D05 (solid
thin blue) and C2Ex (solid thin gray) models convolved with the distribution
function. As already stated in Sec. 3, at $20$ GHz both models are equivalent
even for polarized number counts.

In bottom panel of Figure 1 are shown the polarized number counts at 95 GHz
coming from the PACO-ALMA sample of 32 sources as lower green triangles. Given
the paucity of this sample, we re-sample it by means of 1,000 bootstrap-
resampling. The resampled source counts (shown as green lower triangles in
Figure 1) are then computed in a similar manner as for the $20$ GHz
observations and are summarized in the companion paper by Galluzzi et al.
(2018, in prep.). The error bar estimation of each data point include the
Poissonian $68\%$ CL uncertainties (Gehrels, 1986) plus the error derived from
the uncertainties of log-normal parameters
$\delta_{A},\delta_{\mu},\delta_{\sigma}$ (summarized in Table 2). This error
has been assessed by means of differencing the number counts convolved with an
upper and a lower log-normal function, respectively estimated at maximum and
minimum values of log-normal parameters.

We would like to stress that this is the first time that number counts from
the PACO-ALMA sample have been computed and exploited for this kind of
analysis. Notice that the data are very well fitted by both predictions.

The estimated values of $\langle\Pi\rangle$, $\Pi_{med}$ and
$\langle\Pi^{2}\rangle^{1/2}$ for FSRQ(left panel) and SSRQ (right panel) are
shown in Figure 3. By comparing the two panels, we note that the SSRQ
fractional polarizations increases with frequency. Although this could be
simply related to observational bias (at higher frequencies, steep-spectrum
sources contributes at fainter fluxes), such frequency dependence of $\Pi$ for
SSRQs has been already discussed in Tucci & Toffolatti (2012). On the
contrary, the fractional polarization measured for the FSRQ remains almost
constant during the frequency range studied. To quantify this dependence, we
estimate a linear fit on $\langle\Pi^{2}\rangle^{1/2}$ as a function of a wide
(around 2 orders of magnitude) range of frequencies. This choice is mainly due
to the fact that $\langle\Pi^{2}\rangle$ values are needed to estimate B-mode
angular power spectrum of polarized ERSs and we include in the linear fit also
the values of $\langle\Pi^{2}\rangle^{1/2}$ estimated by Bonavera et al.
(2017a) between 30 and 217 GHz. They were derived assuming a log-normal
distribution as in this work. In particular, for the best fit, we retain only
fractional polarization from the FSRQs and BL Lacs since their contribution
dominates number counts at larger fluxes and at frequencies $>20$ GHz (see
Figures 1 and A.1). The linear fit involves the data for which the estimation
of $\mu$ and $\sigma$ are reliable (filled symbols in Figure 3). Open symbols
indicate data that have not been included to the fit, mainly because of the
poor statistics in fitting the log-normal distribution (e.g. less than 20
polarized sources have been detected in polarization in the Planck HFI
channels, see Table2).

We find a negligible frequency dependence of $\langle\Pi^{2}\rangle^{1/2}$:

$\displaystyle\langle\Pi^{2}\rangle^{1/2}(\nu)$
$\displaystyle=\left({0.005}\pm 0.006\mathrm{GHz^{-1}}\right)\nu$ (6)
$\displaystyle+\left({4.170}\pm 0.22\right).$

In the top left panel of Figure 3 we show the linear fit as a gray solid line
with darker and lighter shaded areas resembling respectively the $1\sigma$ and
$2\sigma$ uncertainties on best fit parameters. Notice that for $\nu>20$ GHz,
we found $\langle\Pi^{2}\rangle^{1/2}\sim 4\%$, in agreement with the value
found by Tucci & Toffolatti (2012) and consistent with the expectations of
Tucci et al. (2004) and Stil et al. (2014).

At l $\nu<20$ GHz, SSRQs have to be taken into account to forecast the
contribution of ERS to CMB observations. Thus, we perform the same linear fit
by including SSRQs for all the datasets at frequencies smaller than $20$ GHz,
shown in Figure 3 (top right panel). The best fit equation changes into

$\displaystyle\langle\Pi^{2}\rangle^{1/2}(\nu)$ $\displaystyle=\left(-0.015\pm
0.009)\mathrm{\,GHz}^{-1}\right)\nu$ $\displaystyle+\left(5.43\pm
0.23\right).$

Nonetheless the slope is still negligible, the presence of SSRQs enhances the
average polarization fraction of sources at frequencies $\nu\lesssim 20$ GHz
and, as one can notice in Figure 3, this is consistently observed in
${\langle\Pi\rangle}$ as well.

We would like to stress that selection effects could bias our results towards
larger values of $\Pi$, especially where few tens of polarized sources have
been detected, see Table 2. This is the reason why we excluded PCCS2 HFI data
(magenta diamonds) in Figure 3 and we considered the ones from Bonavera et al.
(2017a) (gray pentagons). To this regard, the stacking technique helps because
it includes the faint sources to the statistical estimate of $\Pi$ even if
those sources are not directly detectable.

Figure 3.—  Values of ${\langle\Pi^{2}\rangle}^{1/2}$ (top),
$\langle\Pi\rangle$ (center) and $\Pi_{med}$ (bottom) derived from best fit
lognormal parameters (4), (3), (5). Open symbols are data which have not been
taken into account for the linear fit in (6). We distinguished FSRQs (left)
from the SSRQs (right). In top left panel, a linear function is fitted from
the data to provide a scaling of polarization fraction as a function of
frequency. Light and dark shaded area are respectively the $1\sigma$ and
$2\sigma$ uncertainties.

##  5\. Forecasts for forthcoming CMB ground-based experiment

In this section we present the forecast analysis for current and forthcoming
CMB surveys performed with a Python package Point Source ForeCast (PS4C) made
publicly available111111https://gitlab.com/giuse.puglisi/PS4C. PS4C is a user
friendly platform which allows to forecast the contribution of radio point
sources both in total intensity and polarized flux-densities given the nominal
specifics of a CMB experiment. In Table3 we summarize the specifics of 5 CMB
experiments with whom we forecast the ERS contribution with PS4C:

  * •

the Q-U-I JOint TEnerife López-Caniego et al. (QUIJOTE 2014) CMB experiment
designed to observe the polarized emissions from the CMB, our Galaxy and the
extra-galactic sources at four frequencies in the range between 10 and 20 GHz
and at FWHM resolution of $\sim 1^{\circ}$. Observations started observing in
November 2012, covering $18,000\,\mathrm{deg}^{2}$ of the Northern hemisphere,
and achieved the sensitivity of $1800{\,\mu}$K arcmin in polarization;

  * •

a generic CMB-S2 experiment observing at $95,\,150$ GHz within a patch
including $2\%$ of the sky at the resolution of $3.5$ arcmin, at $25\div
30\,\mu\mathrm{K\,arcmin}$ sensitivity;

  * •

a CMB-S3 ground based experiment with the so-called _strawman_ configuration,
as it has been defined in Abazajian et al. (2016), for the “measuring-r”
survey. It consists of an array of small-aperture (SA, $\sim 1$ m) telescopes
and one large-aperture (LA, $\sim 5$ m) telescope, observing at the accessible
atmospheric windows in the sub-millimeter range (at about 30, 40, 90, 150
GHz). The sensitivities at these frequencies are targeted to be about $1\div
10\,\mu$K arcmin.

  * •

the Lite satellite for the studies of B-mode polarization and Inflation from
cosmic Background Radiation Detection (LiteBIRD Matsumura et al., 2016) is a
satellite mission proposed to JAXA aimed at measuring the CMB polarized signal
at degree angular scale. Its goal is to characterize the measurement of $r$
with an uncertainty $\sigma(r)<0.001$. In order to achieve such high accuracy,
the target detector sensitivity is $2\,\mu$K arcmin observing over a wide
range of frequencies (from $40$ to $320$ GHz). The current effort aims to
launch in 2025;

  * •

the Cosmic ORigin Explorer (Delabrouille et al., 2017, CORE) is a next
generation space-borne experiment and it has been proposed as a Medium-size
ESA mission opportunity. It has been designed as the Planck satellite
successor, planned to have better angular resolution and sensitivity than
Planck. We consider the _CORE150_ configuration: a satellite involving a $1.5$
m telescope, observing over a wide range of frequency channels (up to $800$
GHz) with sensitivities ranging from $\sim 10$ to $5\,\mu\mathrm{K\,arcmin}$.
In this work, we restrict our analysis to a selection of frequency channels,
(see the last row of Table 3) to compare the expectations with the ones
previously obtained by De Zotti et al. (2016).

Table 3Nominal specifics of CMB experiments described in Sec.5.  | Frequency [GHz] |  Sensitivity $\left[\mu\mathrm{K\,arcmin}\right]$ | FWHM | $f_{sky}$
---|---|---|---|---
QUIJOTE | 11,13,17,19 | 1800 | $1^{\circ}$ | $50\%$
CMB-S2 | 95, 150 | 25,30 | $3.5^{\prime}$ | $5\%$
CMB-S3 SA | 30, 40, 95,150 | $8,\,6,\,1,\,2\,$ | $1^{\circ}$ | $20\%$
CMB-S3 LA | 30, 40, 95,150 | $8,\,6,\,1,\,2\,$ | $10^{\prime},\,7^{\prime},\,3^{\prime},\,2^{\prime}$ | $20\%$
LiteBIRD | 40, 50, 60, 68, 78 | $53,32,25,19,15$ | $1^{\circ}$ | $100\%$
89, 100,119, 140,166 | $12,15.6,12.6,8.3,8.7$ | $1^{\circ}$ | $100\%$
CORE150 | $60,100,145$ | $10.6,7.1,5.1$ | $14^{\prime},8^{\prime},6^{\prime}$ | $100\%$
Table 4 Number of polarized ERS detected above the $P_{3\sigma}$ flux density detection limit in polarization, by current and forthcoming CMB ground based experiments. Counts are estimated both from the D05 and the C2Ex predictions (in brackets). | CMB -S2 | CMB -S3
---|---|---
|  |  | SA |  | LA |
$\nu$[GHz] |  $P_{3\sigma}$ [mJy] | $N_{3\sigma}$ |  $P_{3\sigma}$ [mJy] | $N_{3\sigma}$ |  $P_{3\sigma}$ [mJy] | $N_{3\sigma}$
30 | … | … | 15 | 236 (191) | 1.5 | 2329 (2278)
40 | … | … | 15 | 215 (156) | 1.5 | 1867 (1810)
95 | 100 | 3 (2) | 10 | 355 (222) | 1 | 2432 (2136)
150 | 100 | 3 (1) | 15 | 146 (74) | 1.5 | 1145 (867)
Table 5 Number of sources detected above the $\geq S_{\lim}$ and $\geq P_{lim}$ flux densities limit by the QUIJOTE experiment, assuming the nominal and conservative values for sensitivity. Values are estimated using D05 and C2Ex models (ins brackets).  $\nu$[GHz] |  $S_{lim}$ [Jy] | $N_{src}$ |  $P_{lim}$ [Jy] | $N_{src}$
---|---|---|---|---
11 | 0.5 | 694 (673) | 0.5 | 6 (4)
1 | 347 (340) | 1 | 2 (1)
13 | 0.5 | 445 (434) | 0.5 | 2 (1)
1 | 210 (205) | 1 | 0 (0)
17 | 1 | 201 (197) | 1 | 0 (0)
2 | 86 (83) | 2 | 0 (0)
19 | 1 | 128 (125) | 1 | 0 (0)
2 | 52 (51) | 2 | 0 (0)

Although most of the frequency channels of future experiments range up to
$350$ GHz, we forecast up to $150$ GHz. This is because at higher frequencies
the contribution coming from dusty galaxies and Cosmic Infrared Background
cannot be neglected121212 We have already planned to include into the package
the contribution from dusty galaxies and forecasts with PS4C will be presented
in a future release that will be described in a future paper. (Negrello et
al., 2013; De Zotti et al., 2016). Bonavera et al. (2017b) estimated the
polarized contribution of dusty galaxies by stacking about 4700 sources
observed by Planck at $143,\,217,\,353$ GHz HFI channels. They estimated the
polarized contribution of dusty galaxies to B-mode power spectra and found
that at frequencies larger than $217$ GHz these population of sources might
remarkably contaminate the primordial B-modes.

Table 6Number of sources observed above $3\sigma_{det}$ limit in terms of polarized flux density $P_{3\sigma}$ by the LiteBIRD experiment. Bracketed values are estimated using the C2Ex model.  $\nu$ [GHz] |  $P_{3\sigma}$ [mJy] | $N_{3\sigma}$
---|---|---
40 | 450 | 4 (3)
50 | 240 | 11 (8)
60 | 210 | 9 (6)
68 | 300 | 4 (3)
78 | 240 | 6 (4)
89 | 210 | 12 (8)
100 | 240 | 10 (7)
119 | 210 | 14 (10)
140 | 270 | 8 (4)
166 | 270 | 7 (4)

We compute one realization of CMB power spectra by means of the CAMB package
(Lewis et al., 2000) by assuming the Planck best fit cosmological parameters
(Planck Collaboration et al., 2016c) and a tensor to scalar ratio $r=0.05$
(slightly below the current upper limits).

To assess the contribution of ERS to the power spectrum level, we assume their
distribution in the sky to be Poissonian, since the contribution of clustering
starts to be relevant for $S<10$ mJy (González-Nuevo et al., 2005; Toffolatti
et al., 2005). The power spectrum of temperature fluctuations coming from a
Poissonian distribution of sources is expected to be a constant contribution
at all multipoles. In particular, we consider as _masked_ all sources whose
flux-density is above $3\sigma$ the detection limit $S_{cut}=3\sigma_{det}$
and we do not include them into power spectrum estimate

$C_{\ell}^{T}=\left(\frac{dB}{dT}\right)^{-2}N\langle
S^{2}\rangle=\left(\frac{dB}{dT}\right)^{-2}\int_{0}^{S_{cut}}n(S)S^{2}dS,$
(7)

where $n(S)$ and $N$ are respectively the differential and the integral number
counts per steradian, and $dB/dT$ is the conversion factor from brightness to
temperature, being

$\left(\frac{dB}{dT}\right)^{-1}\approx
10^{-2}\frac{(e^{x}-1)^{2}}{x^{4}e^{x}}\frac{\mu\mathrm{K}}{\mathrm{Jy\,sr^{-1}}},$

with $x=\nu/57$ GHz. Tucci et al. (2004) found that it possible to relate the
ERS polarization power spectrum to the intensity one (7) as follows

$\displaystyle C_{\ell}^{Q}=$
$\displaystyle\left(\frac{dB}{dT}\right)^{-2}N\langle Q^{2}\rangle$
$\displaystyle=$ $\displaystyle\left(\frac{dB}{dT}\right)^{-2}N\langle
S^{2}\Pi^{2}\cos^{2}2\phi\rangle$ $\displaystyle=$
$\displaystyle\left(\frac{dB}{dT}\right)^{-2}N\langle
S^{2}\rangle\langle\Pi^{2}\rangle\langle\cos^{2}2\phi\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{dB}{dT}\right)^{-2}\langle\Pi^{2}\rangle
C_{\ell}^{T},$ (8)

where the $1/2$ factor comes from the average value of $\cos^{2}2\phi$, if the
polarization angle $\phi$ is uniformly distributed. The value for
$\langle\Pi^{2}\rangle$ is derived at each frequency from eq.(6). Since we do
expect point sources to equally contribute on average both to $Q$ and $U$, and
thus to the $E$\- and $B$\- modes, we can approximate $C_{\ell}^{B}\simeq
C_{\ell}^{E}\simeq C_{\ell}^{U}\simeq C_{\ell}^{Q}$. In the following, B-mode
power spectra are normalized by the usual normalization factor
$\mathcal{D}_{\ell}=\ell(\ell+1)C_{\ell}/2\pi$.

To forecast the number of sources that will be detected in intensity and
polarized flux-density above a given detection limit, we integrate the
differential number counts, $n(S)$ and $n(P)$ as

$\displaystyle N(>S)=$ $\displaystyle\int_{S_{cut}}^{\infty}n(S)dS,$ (9)
$\displaystyle N(>P)=$ $\displaystyle\int_{P_{cut}}^{\infty}n(P)dP.$ (10)

Finally, to compare the level of contamination produced by the ERS with the
Galactic foreground one, we rescale the Galactic foreground emission at a
given $f_{sky}$, frequency $\nu$ and multipole order $\ell$ as in Planck
Collaboration et al. (2016b),

$\displaystyle\mathcal{D}^{FG}(\ell,\nu,f_{sky})=$
$\displaystyle\frac{Var\left[\mathrm{Sync},f_{sky}\right]}{Var\left[\mathrm{Sync},f_{sky,0}\right]}q_{s}\left(\frac{\ell}{80}\right)^{\alpha_{s}}\frac{s_{s}(\nu)}{s_{s}(\nu_{s})}+$
$\displaystyle\frac{Var\left[\mathrm{Dust},f_{sky}\right]}{Var\left[\mathrm{Dust},f_{sky,0}\right]}q_{d}\left(\frac{\ell}{80}\right)^{\alpha_{d}}\frac{s_{d}(\nu)}{s_{d}(\nu_{d})}.$
(11)

with $s,d$ referring respectively to synchrotron and dust. For all the
parameters entering in (11), we use the best fit values quoted in Planck
Collaboration et al. (table 11 2016b) estimated outside the Galactic plane in
the UPB77 mask (Planck Collaboration et al., 2016a, defined in section 4.2).
The mask has been computed considering a common foreground mask after
component separation analysis with $1^{\circ}$ apodization scale. Therefore,
to rescale the estimate in eq.(11)to a patch with a smaller fraction of sky,
$f_{sky}$, we need to compute the variance of both synchrotron and thermal
dust template maps within the considered patch and within the Planck region
with $f_{sky,0}=73\%$. The rescaled foreground power spectra are shown in
Figure 4 as dotted lines.

###  5.1. PS4C with current and forthcoming CMB ground based experiments

Figure 4 shows our PS4C forecasts of foreground contamination to the recovery
of the CMB B-mode for the different experiments in the different panels: we
plot the expected spectrum in polarization of Galactic (dotted lines) and ERS
(dashed lines) emissions at the different frequencies available for each
experiment and the total CMB B-mode power spectrum (black solid line). The
black dot-dashed lines show the primordial ($r=0.05$) and lensed B-mode power
spectra separately. The power spectra are computed in the region outside the
UPB77  Planck mask (in order to exclude the Galactic plane and the ERS whose
flux density is below the $3\sigma$ detection limit). The Galactic foreground
turns out to be the most contaminating emission in the B-mode recovery. The
different colors for the Galactic and ERS spectra are for different
frequencies, going from purple to yellow as the frequency increases. It should
be commented that there exists several _component separation_ and _foreground
cleaning_ algorithms that can recover CMB intensity and polarization signals
with great accuracy (Planck Collaboration et al., 2016b). In addition, multi-
frequency observations and joint analyses from different experiments
(BICEP2/Keck and Planck Collaborations et al., 2015) can improve the
foreground cleaning. So, even if in our work we are considering the most
conservative cases, it should be stressed that such contamination could be
lowered ( at sub-percentage level Stompor et al., 2016; Errard et al., 2011)
by applying such foreground removal algorithms.

In particular, Figure 4 shows our forecasts for the QUIJOTE (top left) and
CMB-S2 (top right) experiments. As for QUIJOTE, the Galactic emission is much
higher than the CMB one and higher than the contribution from undetected ERS,
except at small angular scales where the ERS start to be dominant. Since the
QUIJOTE experiment ranges from $10$ to $20$ GHz, we need to take into account
the contribution from both FSRQs and SSRQs, with the resulting increase in the
average fractional polarization and number counts (see Figure 3 and Figure
A.1). Table 5 summarizes the total number of sources in total intensity (third
column) and polarization (fourth column) that QUIJOTE would detect
(frequencies are given in the first column), assuming nominal and conservative
sensitivity values (flux density limits in total intensity and polarization
are listed in columns two and three respectively). We found $694$, $445$,
$201$ and $128$ sources in total intensity at $11$, $13$, $17$, $19$ GHz
respectively. In polarization only a few of them would be detected and just in
the $11$ and $13$ GHz channels.

Figure 4.— Forecasts of foreground contamination with PS4C. In all panels, the
black dot-dashed lines show the primordial ($r=0.05$) and lensed CMB B-mode
power spectra and the black solid line is the the total CMB B-mode power
spectrum. The dotted (dashed) lines are the power spectrum of the polarized
Galactic emission (ERS emission) at the different frequencies available for
each experiment, the color scale is such that the colors go from purple to
yellow as the frequency increases. The power spectra depend are estimated
using eq.(11) in the region outside the UPB77  Planck mask (in order to
exclude the Galactic plane and ERS above the $3\sigma$ detection limit). The
different panels corresponds to predictions for different experiments. From
top to bottom and from left to right: QUJOTE ($11$, $13$, $17$, $19$ GHz),
CMB-S2 ($95$ and $150$ GHz), CMB-S3 observing with small and large aperture
telescopes ($30$, $40$, $95$, $150$ GHz), LiteBIRD (frequencies between
$40-166$ GHz) and CORE150 ($60,\,100,\,145$ GHz.

For the CMB-S2 experiment whose frequencies are greater than 95 GHz, the
Galactic emission (mostly thermal dust emission) is the most contaminating up
to $\ell\sim 350$, while the ERS are important at small angular scales. Unlike
the previous case, at these frequencies the CMB B-mode spectrum is comparable
to the one of undetected ERS.

In Figure 5 the triangles show the $C_{\ell}^{BB}$ of undetected ERS estimated
using eq.(8). The detection limits are given by the CMB-S2 sensitivities. The
$C_{\ell}^{BB}$ of the CMB B-mode are also plotted: the cyan dashed line is
for the case $\ell\approx 80$ and $r=0.05$ and the orange dashed line is for
$\ell\approx 1000$. Figure 5 shows what is the contamination due to undetected
ERS and consequently the level of source detection required to detect
primordial or lensing $B$-mode signal. In CMB-S2 the undetected ERS level of
the power spectrum is comparable to the lensing B-mode one. In this case,
given the experiment sensitivity and the size of the observed region, $\sim
150$ sources would be detected in total intensity and only few of them in
polarization at a $3\sigma$ level.

Among the experiments studied in this work, the CMB-S3 is the one with the
greatest sensitivity and best resolution. The results are shown in the central
panels of Figure 4 and in the left panel of Figure 5 with circles and
diamonds. As summarized in Table 4, the maximum number of polarized sources
detected above a $3\sigma$ level and using the large aperture telescope is
$2329$ with flux density $P_{lim}\gtrsim 1$ mJy. When using a smaller aperture
telescope, this number drops to a few hundreds with polarized flux densities
$P_{lim}\gtrsim 10$ mJy.

The contribution in polarization of undetected ERS is very small at high
frequencies($\nu\gtrsim 90$) and at low multipoles $\ell\lesssim 2000$. At
lower frequencies, undetected ERS still can contaminate and they have to be
taken into account to de-lens, lensing $B$-modes to get the primordial ones
for $r\lesssim 0.05$.

Figure 5.— Power spectra in polarization of undetected ERS in current and
future CMB experiments. Left panel: CMB-S2 (triangles) and CMB-S3 (circles for
the small aperture telescope and diamonds for the large aperture telescope).
Right panel: LiteBIRD (circles) and CORE150 (squares). The dotted lines are
the B-mode power spectra at the acoustic scale ($\ell=80$) and at the lensing
B-modes peak scale ($\ell\approx 1000$).

###  5.2. PS4C with future space missions

The results for the LiteBIRD experiment are shown in the left bottom panel of
Figure 4 and the filled circles in the right panel of Figure 5. On the whole,
the most contaminating contribution is the Galactic one, except at small
angular scales ($l\sim 400$) and high frequencies ($\nu>70$ GHz) where the ERS
contribution is comparable to the Galactic one. The ERS contribution, although
generally lower than the Galactic one, is also important being higher than the
CMB B-mode level even at large scales ($l\gtrsim 7$) and $\nu<70$ GHz (dashed
purple and blue lines). Moreover, at $\nu>80$ GHz and $l\gtrsim 70$ the ERS
contribution is comparable to the B-mode power spectrum. The number of sources
that would be detected in polarization above the $3\sigma$ level with this
experiment are listed in Table 6 and they range from 4 at 10 and 68 GHz to 14
at 119 GHz. The first column is the frequency in GHz, the second is the
polarized flux density limit in mJy and the third column is the number of
sources that would be detected by LiteBIRD (values in the brackets are
estimated from the C2Ex model).

Our findings for CORE are shown in the right bottom panel of Figures 4 and in
the right panel of Figure 5 (squares). Galactic emission is the most
contaminating for B-mode detection. Undetected ERS are important only at 60
GHz, where their power spectrum is comparable to the one of the B-mode due to
lensing. CORE would be able to detect up to 200 sources per steradian,
implying a lower contamination for the CMB B-mode power spectrum with respect
to LiteBIRD.

Table 7 compares the surface densities (i.e. number of sources per steradian,
last two columns) at CORE frequencies (first column) of the polarized ERS
above the $P_{4\sigma}$ flux density limit (second column) estimated by De
Zotti et al. (2016)(DZ16) and PS4C (values in the brackets are for C3Ex
estimate). In this comparison we use a $4\sigma$ flux density limit in order
to be consistent with the estimates by De Zotti et al. (2016). Above $100$
GHz, we find a discrepancy between D05 and DZ16 that could be due to two
effects that become more important at higher frequencies: (i) the D05
predictions tend to over-estimate the polarized source number counts (see Sec.
3) and (ii) at $\nu>100$ the polarization fraction is expected to suffer a
slight increase (from $\sim 4\%$ to $\sim 5\%$ from $100$ to $150$ GHz) as can
be seen in eq. (6) and Figure 3.

On one hand, at $100$ GHz, we find that accounting solely for the observation
in (ii), i.e. a $20\%$ increase of $\Pi$ to a value of $4.67\%$, the D05
forecasts predict source counts that are $20\%$ larger than DZ16131313 For
this estimate, we assume that differential source counts are described by a
power law with spectral index $>1$. On the other hand, at $150$ GHz, the
surface density estimated with PS4C with D05 model is $\sim 65\%$ larger than
the value referred by DZ16. Even accounting for the $25\%$ fractional increase
of $\Pi$ to $4.92\%$ from eq.(6), this is not enough to compensate the
observed discrepancy. We thus argue that the discrepancy at $150$ GHz is
caused by both (i) and (ii).

Contrary to the D05 forecasts, the C2Ex model is in reasonable agreement with
De Zotti et al. (2016), meaning that the C2Ex predictions are more robust than
the D05 ones at least at higher frequencies.

Table 7 Comparison of surface densities of polarized ERSs brighter than $P_{4\sigma}$ estimated by De Zotti et al. (2016)(DZ16) and by PS4C. Values in brackets refer to C2Ex estimates.  $\nu$ [GHz] |  $P_{4\sigma}$ [mJy] |  $N_{4\sigma}\,[\,\mathrm{sr^{-1}}\,]$
---|---|---
|  | DZ16 | PS4C
60 | 5.2 | 212 | 214 (198)
100 | 5.2 | 184 | 229 (164)
145 | 4.6 | 165 | 271 (142)

##  6\. Summary and conclusions

We describe and present the state-of-the-art observations on polarization of
ERS over a wide frequency range, namely from $1.4$ to $217$ GHz. We exploit
for the first time the polarization number counts at 95 GHz from a sample of
32 polarized sources detected with ALMA. The characterization of these sources
and their spectral behaviour in frequencies ranging from $1$ to $95$ GHz are
described in a companion paper by Galluzzi et al. (2018, in prep.)

By collecting polarization flux densities from 10 catalogues, we are able to
derive a relation of the average fractional polarization as a function of
frequency and to avoid extrapolations that have been commonly adopted to
forecast the average polarization fraction from low- ($\lesssim 20$ GHz where
enough data have been collected), to high-frequency ($\gtrsim 70$ GHz where
still few polarization measurements have been performed). Therefore, we fit a
linear function on data from several surveys, including Planck measurements
from both detection and stacking, and we find a mild dependence of
$\langle\Pi^{2}\rangle^{1/2}$ as a function of $\nu$.

This relation allows us to forecasts the contribution of ERSs to polarization
B-mode power spectrum given the nominal sensitivities of current and
forthcoming CMB experiments, by means of predictions of ERS counts coming from
two models D05 and C2Ex. The whole forecast suite is fully integrated into a
Python package, PS4C, made publicly available with online documentation and
tutorials.

We discuss the reasons why we do not assume a correlation between the level of
fractional polarization and the total intensity flux. Although still
controversial and not observed at high-radio frequencies (Galluzzi et al. 2018
in prep., Galluzzi et al., 2018, 2017; Massardi et al., 2013), deeper surveys
in polarization are critical to further proof the validity of this assumption,
not only at higher frequencies but also at fainter flux density levels.

Future CMB experiments could shed light on this interesting aspect: in fact,
we have shown that they are going to observe an increasing number of polarized
ERS (they are foreseen to detect up to $\sim 2000$ polarized ERS) because
their sensitivity will increasingly improve in the future.

A further potentiality of future CMB experiments is that they can be largely
exploited by the community as _wide global surveys_ to measure polarized flux
density of sources at very high-radio frequencies (Partridge et al., 2017).
Programs aimed at observing ERSs at higher resolution can thus benefit of CMB
large area surveys in an extremely wide range of frequencies, from $20$ up to
$300$ GHz.

Moreover, since in this work we mostly focus on blazar statistical
polarization, as it is the main bright source population at frequencies $<150$
GHz, we restrict our forecast analysis up to this frequency limit. At higher
frequencies, the far-IR dusty star forming galaxies constitute the majority of
extra-galactic sources (see Planck Collaboration (fig.25 2015)) and, similarly
to the ERSs, their polarized emission contaminates B-mode power
spectra141414In addition to the Poissonian contribution, an extra-term coming
from clustering has to be considered when dusty sources are involved. (De
Zotti et al., 2015). Recent works from Bonavera et al. (2017b); De Zotti et
al. (2016) have already shown statistical polarization properties of dusty
sources and forecasted their contribution for future CMB experiments.
Therefore, we plan to include those estimates within the PS4C package in a
future development.

As a final remark, we stress that ERSs below the detection flux limit may
introduce a bias at all the angular scales and at frequencies $\nu<50$ GHz:
their synchrotron emission is still strong enough to contaminate polarization
measurements even at low flux densities, namely $P<1$ mJy. At larger
frequencies, ERS polarization power spectra have to be assessed as long as
smaller angular scales are involved to estimate the CMB power spectrum at
multipoles around the lensing peak or to estimate the primordial B-mode power
spectrum at lower multipoles ($\ell<800$) by means of de-lensing algorithms.

_Acknowledgments._ We are pleased to thank Matteo Bonato, Bruce Partridge,
Vincent Pelgrims and the referee for their helpful comments and suggestions.
This work was supported by the RADIOFOREGROUNDS grant of the European Union’s
Horizon 2020 research and innovation programme (COMPET-05-2015, grant
agreement number 687312) and by by National Institute of Nuclear Physics
(INFN) INDARK. LB and JGN acknowledge financial support from the I+D 2015
project AYA2015-65887-P (MINECO/FEDER). J.G.N also acknowledges financial
support from the Spanish MINECO for a ‘Ramon y Cajal’ fellowship
(RYC-2013-13256). VG and MM acknowledge financial support by the Italian
Ministero dell’Istruzione, Università e Ricerca through the grant Progetti
Premiali 2012-iALMA (CUP C52I13000140001). We further acknowledge support from
the ASI-COSMOS Network151515www.cosmosnet.it.

## References

  * Abazajian et al. (2016) Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv e-prints, arXiv:1610.02743
  * Agudo et al. (2010) Agudo, I., Thum, C., Wiesemeyer, H., & Krichbaum, T. P. 2010, The Astrophysical Journal Supplement Series, 189, 1
  * Arnold et al. (2014) Arnold, K., Stebor, N., Ade, P. A. R., et al. 2014, in Proc. SPIE, Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, 91531F
  * Austermann et al. (2009) Austermann, J. E., Aretxaga, I., Hughes, D. H., et al. 2009, MNRAS, 393, 1573
  * Battye et al. (2011) Battye, R. A., Browne, I. W. A., Peel, M. W., Jackson, N. J., & Dickinson, C. 2011, MNRAS, 413, 132
  * Benson et al. (2014) Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in Proc. SPIE, Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, 91531P
  * BICEP2/Keck and Planck Collaborations et al. (2015) BICEP2/Keck and Planck Collaborations, Ade, P. A. R., Aghanim, N., et al. 2015, Phys. Rev. Lett., 114, 101301
  * Blandford & Koenigl (1979) Blandford, R. D., & Koenigl, A. 1979, Astrophys. Lett., 20, 15
  * Bonaldi et al. (2013) Bonaldi, A., Bonavera, L., Massardi, M., & De Zotti, G. 2013, MNRAS, 428, 1845
  * Bonato et al. (2017) Bonato, M., Negrello, M., Mancuso, C., et al. 2017, Monthly Notices of the Royal Astronomical Society, Volume 469, Issue 2, p.1912-1923, 469, 1912
  * Bonavera et al. (2017a) Bonavera, L., González-Nuevo, J., Argüeso, F., & Toffolatti, L. 2017a, MNRAS, 469, 2401
  * Bonavera et al. (2017b) Bonavera, L., González-Nuevo, J., De Marco, B., Argüeso, F., & Toffolatti, L. 2017b, MNRAS, 472, 628
  * Bonavera et al. (2011) Bonavera, L., Massardi, M., Bonaldi, A., et al. 2011, MNRAS, 416, 559
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, The Astronomical Journal, 8065, 1693
  * Curto et al. (2013) Curto, A., Tucci, M., González-Nuevo, J., et al. 2013, MNRAS, 432, 728
  * De Zotti et al. (2010) De Zotti, G., Massardi, M., Negrello, M., & Wall, J. 2010, Radio and millimeter continuum surveys and their astrophysical implications, 0908.1896
  * de Zotti et al. (2005) de Zotti, G., Ricci, R., Mesa, D., et al. 2005, A&A, 431, 893
  * De Zotti et al. (2015) De Zotti, G., Castex, G., González-Nuevo, J., et al. 2015, JCAP, 6, 018
  * De Zotti et al. (2016) De Zotti, G., Gonzalez-Nuevo, J., Lopez-Caniego, M., et al. 2016, ArXiv e-prints, arXiv:1609.07263
  * Delabrouille et al. (2017) Delabrouille, J., de Bernardis, P., Bouchet, F. R., et al. 2017, ArXiv e-prints, arXiv:1706.04516
  * Errard et al. (2011) Errard, J., Stivoli, F., & Stompor, R. 2011, Phys. Rev. D, 84, 069907
  * Galluzzi & Massardi (2016) Galluzzi, V., & Massardi, M. 2016, International Journal of Modern Physics D, 25, 1640005
  * Galluzzi et al. (2017) Galluzzi, V., Massardi, M., Bonaldi, A., et al. 2017, Monthly Notices of the Royal Astronomical Society, 465, 4085
  * Galluzzi et al. (2018) Galluzzi, V., Massardi, M., Bonaldi, A., et al. 2018, MNRAS, 475, 1306
  * Gehrels (1986) Gehrels, N. 1986, The Astrophysical Journal, 303, 336
  * Goldreich & Kylafis (1981) Goldreich, P., & Kylafis, N. D. 1981, The Astrophysical Journal, 243, L75
  * González-Nuevo et al. (2005) González-Nuevo, J., Toffolatti, L., & Argüeso, F. 2005, ApJ, 621, 1
  * Guth (1981) Guth, A. H. 1981, Phys. Rev. D, 23, 347
  * Henderson et al. (2016) Henderson, S. W., Allison, R., Austermann, J., et al. 2016, Journal of Low Temperature Physics, 184, 772
  * Hu & White (1997) Hu, W., & White, M. J. 1997, New Astron., 2, 323
  * Jackson et al. (2007) Jackson, N., Battye, R. A., Browne, I. W. A., et al. 2007, Monthly Notices of the Royal Astronomical Society, 376, 371
  * Jackson et al. (2010) Jackson, N., Browne, I. W. A., Battye, R. A., Gabuzda, D., & Taylor, A. C. 2010, MNRAS, 401, 1388
  * Keisler et al. (2015) Keisler, R., Hoover, S., Harrington, N., et al. 2015, The Astrophysical Journal, 807, 151
  * Kellermann (1966) Kellermann, K. I. 1966, ApJ, 146, 621
  * Lamee et al. (2016) Lamee, M., Rudnick, L., Farnes, J. S., et al. 2016, ApJ, 829, 5
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, Astrophys. J., 538, 473
  * Lopez-Caniego et al. (2009) Lopez-Caniego, M., Massardi, M., Gonzalez-Nuevo, J., et al. 2009, The Astrophysical Journal, Volume 705, Issue 1, pp. 868-876 (2009)., 705, 868
  * López-Caniego et al. (2014) López-Caniego, M., Rebolo, R., Aguiar, M., et al. 2014, ArXiv e-prints, arXiv:1401.4690
  * Louis et al. (2017) Louis, T., Grace, E., Hasselfield, M., et al. 2017, JCAP, 6, 031
  * Mancuso et al. (2017) Mancuso, C., Lapi, A., Prandoni, I., et al. 2017, ApJ, 842, 95
  * Marriage et al. (2011) Marriage, T. A., Baptiste Juin, J., Lin, Y.-T., et al. 2011, ApJ, 731, 100
  * Massardi et al. (2011) Massardi, M., Bonaldi, A., Bonavera, L., et al. 2011, MNRAS, 415, 1597
  * Massardi et al. (2009) Massardi, M., López-Caniego, M., González-Nuevo, J., et al. 2009, MNRAS, 392, 733
  * Massardi et al. (2008) Massardi, M., Ekers, R. D., Murphy, T., et al. 2008, Mon. Not. R. Astron. Soc, 384, 775
  * Massardi et al. (2011) —. 2011, Monthly Notices of the Royal Astronomical Society, 412, 318
  * Massardi et al. (2013) Massardi, M., Burke-Spolaor, S. G., Murphy, T., et al. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2915
  * Matsumura et al. (2016) Matsumura, T., Akiba, Y., Arnold, K., et al. 2016, Journal of Low Temperature Physics, 184, 824
  * Mocanu et al. (2013) Mocanu, L. M., Crawford, T. M., Vieira, J. D., et al. 2013, ApJ, 779, 61
  * Murphy et al. (2010) Murphy, T., Sadler, E. M., Ekers, R. D., et al. 2010, Monthly Notices of the Royal Astronomical Society, 402, 2403
  * Negrello et al. (2013) Negrello, M., Clemens, M., Gonzalez-Nuevo, J., et al. 2013, MNRAS, 429, 1309
  * Partridge et al. (2017) Partridge, B., Bonavera, L., López-Caniego, M., et al. 2017, Galaxies, 5, 47
  * Pelgrims & Hutsemékers (2015) Pelgrims, V., & Hutsemékers, D. 2015, Monthly Notices of the Royal Astronomical Society, 450, 4161
  * Planck Collaboration (2015) Planck Collaboration. 2015, Astronomy & Astrophysics, 594, A26
  * Planck Collaboration et al. (2011a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011a, Astronomy & Astrophysics, 536, A13
  * Planck Collaboration et al. (2011b) Planck Collaboration, Aatrokoski, J., Ade, P. A. R., et al. 2011b, Astronomy & Astrophysics, 536, A15
  * Planck Collaboration et al. (2013) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, Astronomy & Astrophysics, 550, A133
  * Planck Collaboration et al. (2016a) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016a, A&A, 594, A9
  * Planck Collaboration et al. (2016b) —. 2016b, A&A, 594, A10
  * Planck Collaboration et al. (2016c) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016c, A&A, 594, A13
  * Planck Collaboration et al. (2016d) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016d, Astronomy & Astrophysics, 594, A25
  * Planck Collaboration et al. (2016e) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016e, A&A, 586, A133
  * Puglisi et al. (2017) Puglisi, G., Fabbian, G., & Baccigalupi, C. 2017, MNRAS, 469, 2982
  * Sajina et al. (2011) Sajina, A., Partridge, B., Evans, T., et al. 2011, The Astrophysical Journal, 732, 45
  * Seljak & Zaldarriaga (1997) Seljak, U., & Zaldarriaga, M. 1997, Physical Review Letters, 78, 2054
  * Starobinsky (1982) Starobinsky, A. A. 1982, Physics Letters B, 117, 175
  * Stil et al. (2014) Stil, J. M., Keller, B. W., George, S. J., & Taylor, A. R. 2014, The Astrophysical Journal, 787, 99
  * Stompor et al. (2016) Stompor, R., Errard, J., & Poletti, D. 2016, Phys. Rev. D, 94, 083526
  * The Polarbear Collaboration et al. (2017) The Polarbear Collaboration, Ade, P. A. R., Aguilar, M., et al. 2017, ApJS, 1705.02907
  * The Polarbear Collaboration: P. A. R. Ade et al. (2014) The Polarbear Collaboration: P. A. R. Ade, Akiba, Y., Anthony, A. E., et al. 2014, ApJ, 794, 171
  * Toffolatti et al. (2005) Toffolatti, L., Negrello, M., González-Nuevo, J., et al. 2005, A&A, 438, 475
  * Trombetti et al. (2017) Trombetti, T., Burigana, C., De Zotti, G., Galluzzi, V., & Massardi, M. 2017, ArXiv e-prints, arXiv:1712.08412
  * Tucci et al. (2004) Tucci, M., Martínez-González, E., Toffolatti, L., González-Nuevo, J., & De Zotti, G. 2004, MNRAS, 349, 1267
  * Tucci & Toffolatti (2012) Tucci, M., & Toffolatti, L. 2012, Advances in Astronomy, 2012, 624987
  * Tucci et al. (2011) Tucci, M., Toffolatti, L., De Zotti, G., & Martínez-González, E. 2011, Astronomy & Astrophysics, 533, A57
  * Vieira et al. (2010) Vieira, J. D., Crawford, T. M., Switzer, E. R., et al. 2010, The Astrophysical Journal, 719, 763
  * Waldram et al. (2010) Waldram, E. M., Pooley, G. G., Davies, M. L., Grainge, K. J. B., & Scott, P. F. 2010, MNRAS, 404, 1005
  * Waldram et al. (2003) Waldram, E. M., Pooley, G. G., Grainge, K. J. B., et al. 2003, MNRAS, 342, 915
  * Wright et al. (2009) Wright, E. L., Chen, X., Odegard, N., et al. 2009, ApJS, 180, 283

##  Appendix A

Figure A.1.— Euclidean normalized differential number counts on a wide range
of frequencies. The dotted, dashed, dot-dashed and solid lines are
respectively the number counts of BL Lacs, FSRQs, SSRQs and their total
contribution predicted by the D05 model. The thick solid grey line are the
number counts coming from the C2Ex model. Number counts obtained with several
experiments observing at similar frequency channels are also shown: open
circles at 15 and 20 GHz are respectively from Waldram et al. (2003, 2010) and
Massardi et al. (2008); open circles at 30 and 44 GHz resemble counts from the
PACO dataset (Bonavera et al., 2011); upper triangles from 20 to 70 GHz are
data from WMAP5-yr survey (Massardi et al., 2009)); diamonds at 100, 143, 217
GHz are from SPT (Mocanu et al., 2013); stars at 143 GHz are ACT counts
(Marriage et al., 2011); squares at 30, 44, 70, 100, 143, 217 GHz are based on
data from Planck Collaboration et al. (2011a).