Robust paths to primordial $B$ modes in the presence of polarised point source contamination

Introduction

Current CMB experiments are mapping the microwave sky with unprecedented sensitivity, with the primary goal of detecting or further constraining the primordial $B$-mode signal. Since primordial $B$-mode polarisation can only be generated by tensor perturbations in the metric, its detection would provide a powerful probe of primordial gravitational waves. The amplitude of the tensor contribution is commonly parametrised by the tensor-to-scalar ratio, $r$, for which the current upper bound is $r < 0.032$ at 95% confidence level (Tristram et al. 2022). Current and upcoming experiments aim to further improve the sensitivity on $r$, targeting $\sigma(r) \lesssim0.001$ (The Simons Observatory Collaboration 2026; LiteBIRD et al. 2023).

However, a detection of the primordial $B$-mode CMB polarisation would need to be demonstrated to be robust against the contamination by polarised astrophysical foreground signals. These foregrounds arise mainly due to two emission mechanisms: thermal emission from interstellar dust and synchrotron radiation from charged particles accelerated by astrophysical processes. Both types of emission originate both in our own Galaxy as well as in other galaxies. From an observational perspective, however, the Galactic and extragalactic contaminants appear fundamentally different because of our location within the Milky Way.

The foregrounds of Galactic origin appear diffuse across the sky: while they have higher amplitude along the Galactic plane, they also extend to higher Galactic latitudes and are present in any region of the sky. In particular, even at higher Galactic latitudes, the polarised foreground emission is expected to be stronger than the CMB primordial $B$-mode signal (Krachmalnicoff et al. 2016). For these reasons, foreground removal is an unavoidable and delicate step in the data analysis aiming at detecting the primordial $B$ modes. Considerable efforts have therefore been dedicated to characterising the polarised Galactic foreground contamination (e.g. Krachmalnicoff et al. 2016; Akrami et al. 2020; Krachmalnicoff et al. 2018a; Watts et al. 2024) and a wide range of component separation techniques have been proposed to clean the CMB signal from these polarised foregrounds. Including both techniques that aim at modelling the foreground emissions, leveraging on their spectral dependence, either modelling and performing the separation at map-level [87] (Brandt et al. 1994; Eriksen et al. 2006; Stompor et al. 2009; Rizzieri et al. 2026) or at power spectrum-level [60] (Ade et al. 2015; Wolz et al. 2024; Azzoni et al. 2021; Azzoni et al. 2023), as well as methods that do not directly attempt at modelling the foreground signal, for instance leveraging on statistical properties of the different foreground components (Ichiki et al. 2018; Carones et al. 2023; Leloup et al. 2023; Morshed et al. 2024; Tran et al. 2026).

Turning instead to extragalactic polarised foregrounds for CMB observations, two types of sources are relevant: extragalactic radio sources (ERSs), whose polarised emission is dominated by synchrotron radiation (Toffolatti et al. 1998; Tucci et al. 2011; Planck [59] (Rubiño-Martín et al. 2023); and dusty star-forming galaxies (DSFGs), whose polarised emission may arise from thermal radiation by aligned dust grains in the interstellar medium (Bonavera et al. 2017; Shablovinskaia et al. 2025). Over the frequency range and angular scales most relevant to the primordial $B$-mode search, the polarised extragalactic signal is expected to be generally subdominant to Galactic foreground emission and manifest in the CMB maps with a clearly different spatial morphology. In particular, the emission associated with a single galaxy is spatially localised and, depending on the angular resolution and sensitivity of the telescope, it may or not be individually detected and resolved. Because these objects are compact relative to the angular resolution of typical CMB observations, they are commonly referred to as point sources (PSs). In principle, there are two contributions to the angular power spectrum of the emission from each population: a Poisson noise term, as the sources are assumed to be randomly and independently distributed in the sky; and a clustering 1-halo term (note that no 2-halo term is expected, as there is no reason to expect correlations among polarisation of galaxies within different haloes).

The standard treatment of extragalactic sources in CMB data-analyses focuses around two aspects. Sources detected above a given significance or flux threshold are identified, through source-detection algorithms, compiled into catalogues (e.g. [59] (Meyers et al. 2017; Perrott et al. 2021; Vargas et al. 2026)), and masked from the maps. In the analysis of smaller angular scales, where extragalactic source emission becomes increasingly significant relative to the CMB and the other foreground signals, the residual contribution from sources below the masking threshold is then modelled at the power spectrum level [60] (Dunkley et al. 2013; George et al. 2015; Choi et al. 2020; Beringue et al. 2025). Extending this strategy to larger angular scale polarisation measurements is, however, non-trivial. Direct source detection in polarisation is considerably more challenging than in total intensity, owing to the typically low polarisation fractions of extragalactic sources and the correspondingly low signal-to-noise ratio of individual polarised detections (López-Caniego et al. 2009; Argüeso et al. 2011; Herranz et al. 2021). Moreover, the extent to which the spectral power law modelling commonly adopted in small scale analyses remains adequate for primordial $B$-mode searches has not yet been fully established. In this regime, the cosmological signal of interest is extremely faint, and even a subdominant residual extragalactic contribution may become relevant for the inference of cosmological and foreground parameters.

Several studies have forecasted the number counts, polarisation fractions, and angular power spectra of polarised PSs relevant in CMB experiments (e.g. TToffolatti et al. 1998; De Zotti et al. 2005; Tucci et al. 2011; Tucci & Toffolatti 2012; Puglisi et al. 2018; Lagache et al. 2020). These studies generally predict larger polarised contamination from ERS than from DSFG. And find that, depending on the observing frequencies, source masking thresholds, and instrumental configuration, the residual ERS power can become comparable to the primordial $B$-mode signal around the recombination peak for sufficiently small values of $r$, while remaining less relevant at the reionisation peak. Previous studies have also highlighted the potential impact of extragalactic sources on delensing performance (Puglisi et al. 2018; Qu et al. 2025). These conclusions, however, have largely been obtained by comparing the predicted PS power spectra directly with the primordial $B$-mode power spectrum, or by expressing the contamination in terms of an equivalent tensor-to-scalar ratio. Such comparisons provide a useful indication of the characteristic amplitude of the contamination, but do not capture its propagation through a realistic multi-frequency component separation and cosmological parameter inference pipeline. A robust assessment of this propagation is still lacking, as is the development of concrete strategies for the removal of resolved polarised sources, and for the modelling of residual contamination.

In this work we revised such forecasts in view of ongoing and future experiments, such as the Simons Observatory (SO, The Simons Observatory Collaboration 2019), and we perform a full component separation procedure and inference to study the impact of the ERS and DSFG contamination in the data analysis and on the various astrophysical and cosmological parameters. The paper is organized as follows. In Sect. 2, we discuss the model we implement to forecast the number counts and power spectra for the two types of point sources: the general formalism, the masking strategy we adopt, the model for the ERSs and for the DSFGs. In Sect. 3 we discuss the analysis framework, and in particular how we simulate the total sky signal and covariance, and present the component separation framework and parameter inference. In Sect. 4, we present our results, focusing notably on the simulated polarised ERS signal (Sect. 4.1) and its impact on the analysis (Sect. 4.2), as well as the simulated DSFG polarised signal (Sect. 4.3). Finally, we summarise our conclusions in Sect. 5.

Point sources forecasting

Number counts and power spectra

Our goal is to have an estimate for the polarised power spectra due to the point sources that we would measure at all the auto and cross frequencies of the instrument, as well as for the number of sources that we would mask. We are hence interested in forecasting two quantities in particular related to the PS contamination: the number of detected point sources that will be masked, and the power spectrum of the residual point source contamination. We derive expressions for both quantities below. In doing so, we must consider three different frequencies:

• The reference frequency $\nu_*$ at which we compute the differential number counts $n'_t$ of point sources of type $t$,

  $$\frac{dn_t}{dz} \equiv\frac{d^3N_t}{dz\,d\Omega\,dS_{\nu_*}}. \tag*{(1)}$$

  (i.e. the angular number density of sources of type $t$ within an interval of specific flux $S_{\nu_*}$ and redshift $z$ at frequency $\nu_*$).

• The masking frequency $\nu_m$, at which sources are detected and masked. The final point source mask, used at all frequencies, is defined from these sources (see § 2.2).

• The measurement frequencies $\nu$, at which we make maps of the CMB which we use to infer constraints on primordial $B$ modes. In all cases explored here the point source mask is defined at one of these frequency channels (i.e. $\nu_m$ is also one of the measurement frequencies), although this need not be the case.

Number counts. The total number of PS per steradian of type $t$ with a flux above a certain limit $S_{\rm lim}$ in the masking frequency band $\nu_m$ may be calculated as

where $P(S_{\nu_m} > S_{\rm lim}\mid S_{\nu_*},t,z)$ is the probability for a source at redshift $z$ and of type $t$ to have a flux at frequency $\nu_m$ above the threshold, given its flux at frequency $\nu_*$. For a source with a spectral energy distribution (SED) determined by a set of parameters $\boldsymbol{\theta}_t$, this probability is simply a Heaviside function (1 when the flux $S_{\nu_m}$ extrapolated from $S_{\nu_*}$ using its SED is larger than $S_{\lim}$, and 0 otherwise), which we write as $\Theta(S_{\nu_m}(S_{\nu_*},\boldsymbol{\theta}_t,t,z)>S_{\lim})$. Our point source model will assume that different sources within a given type $t$ have spectra characterised by different parameters $\boldsymbol{\theta}_t$, with a distribution $p(\boldsymbol{\theta}_t)$. An example of such parameters is the source’s spectral index. The total probability for type $t$ is given by the average over this distribution:

The total number counts is then the sum of that of individual populations. In some cases, the differential number counts of a given type $t$ may be parametrised as

where $n(S_{\nu_*},z)$ is the total differential number counts at redshift $z$, and $P_{\mathrm{mix}}(t\mid S_{\nu_*})$ is the probability for a source of flux $S_{\nu_*}$ to be of type $t$.

Power spectra. The polarised auto and cross-frequency power spectra due to the residual PS contamination is generally dominated by the Poisson contribution from these sources, with equal contributions in $E$ and $B$ modes. In the simplest case, the power spectrum auto-correlation at frequency $\nu$ is (see e.g. De Zotti et al. 1999; Trombetti et al. 2018; Puglisi et al. 2018; Lagache et al. 2020):

where $\Pi$ is the polarisation fraction with squared mean $\langle\Pi^2\rangle$. We assume that $\Pi$ is independent of frequency and of the flux density $S_\nu$. While data suggest that this may not be the case at frequencies $\lesssim20$ GHz, there is no significant evidence yet for it to be the case at higher frequency (Puglisi et al. 2018; De Zotti et al. 2018; Gupta et al. 2019; Datta et al. 2019). The simple expression above is not sufficiently general for our purposes, as we need to calculate the general expression for cross-frequency correlations, accounting for populations of different types as well as the masking strategy used. Writing explicitly the redshift dependence and accounting for the SEDs of different source populations, as well as the contribution from imperfectly masked sources, the general expression for cross-correlations is:

where $f_\nu(S_\nu,t,\boldsymbol{\theta}_t)$ quantifies the fraction of the flux at frequency $\nu$ for a source of type $t$ with SED characterised by $\boldsymbol{\theta}_t$ that has been missed by our masking criterion. This is described in the next section. Note that the power spectra thus computed correspond to maps provided in units of specific intensity (e.g. Jansky per steradian). We then convert these to thermodynamic temperature units ($\mu\mathrm{K}_{\mathrm{CMB}}$) using the CMB blackbody spectrum for a temperature $T_{\mathrm{CMB}}=2.7255\ \mathrm{K}$. In addition to the Poisson noise, a clustering term may also be in principle relevant, see Sect. 2.4.3.

Note that the flux limit associated with both of these quantities is not given explicitly as an integration limit. Instead, it is encoded in the masking criterion defining $f_\nu$, in the case of the $C_\ell$'s, and in the selection criterion $(S_{\nu_m}>S_{\lim})$ in the case of the number counts of masked sources.

Masking strategy

We assume that the point source mask used in the analysis is based on the positions of the sources detected at a frequency $\nu_m$. Once defined, the same mask is then used to analyse maps at all other frequencies, and to derive $B$-mode constraints after component separation. This is an important aspect of $B$-mode analysis in the context of component separation, as using the same mask across frequencies ensures that the same foreground spatial modes are present in all cross-frequency power spectra. Besides the position of the masked sources, the other element defining the mask is the masking radius, defining the disk removed around each source. We parametrise this radius as

where $\theta_{\mathrm{FWHM}}(\nu_m)$ is the full-width at half-maximum aperture of the instrumental beam at the masking frequency $\nu_m$, and $f_{\mathrm{FWHM}}$ is a parameter quantifying the fraction of the FWHM aperture removed by the mask. The strictness of the mask is then defined by two parameters: $S_{\lim}$ determines the number of sources to mask, and $f_{\mathrm{FWHM}}$ the amount of area removed by masking.

Assuming a Gaussian source profile, the residual flux of a masked source at frequency $\nu$ is

Hence, the residual fraction at frequency $\nu$, $f_\nu$, entering the PS power spectrum (Eq. (6)) is given by:

This quantity depends implicitly on the SED parameters $\boldsymbol{\theta}_t$ through the relation between $S_\nu$ and $S_{\nu_m}$.

Model for polarised extragalactic radio sources

Polarised PS emission from extragalactic radio sources is caused by partially polarised synchrotron radiation sourced by active galactic nuclei (AGN) hosting a central super-massive black hole (SMBH). These objects present narrow collimated jets of material ejected from the SMBH and travelling far from the hosting galaxy. Lobes of material, constantly fed by the jets with new energetic particles, form around the jets and extend up to a few megaparsecs from the central SMBH. The relativistic jets emit Doppler-boosted synchrotron radiation with spectrum $S_\nu\propto\nu^\alpha$ which, in the frequency range $1-5$ GHz is approximately flat, with $\alpha>-0.5$, while the lobes emit in radio with a steeper spectral index $\alpha<-0.5$. Depending on each source’s viewing angle, one of these two contributions dominates: edge-on ERSs are typically dominated by the lobes and classified as steep-spectrum sources, while face-on sources are dominated by the jets and classified as flat-spectrum. Additionally, a less abundant sample of ERSs exists exhibiting spectra with a positive tilt, and is often referred to as the “inverted-spectrum” population. More details can be found in De Zotti et al. (2010); Jackson & Wall (1999) and references cited therein.

For our purposes, we are mostly interested in the spectral properties in the frequency range $20\ \mathrm{GHz}\lesssim\nu\lesssim300\ \mathrm{GHz}$, covered by most CMB experiments. We hence consider three ERS types, depending on the value of their effective spectral index in the frequency range $1$-$5$ GHz, $\alpha_1^5$: “steep” ($\alpha_1^5<-0.5$), “flat” ($-0.5\leq\alpha_1^5<0.3$), and “inverted” (if $\alpha_1^5\geq0.3$).

Because of the scarcity of data on polarised ERSs, and as done in previous analyses (e.g. Puglisi et al. 2018; Lagache et al. 2020), we do not consider models for the abundance of polarised ERS number counts directly. Instead, we use models derived from number counts in intensity, and produce predictions for the residual PS polarised power spectra assuming a given polarisation fraction probability density function $P(\Pi)$. Similarly, and following the approach of Tucci et al. (2011); Lagache et al. (2020), due to the wealth of data at low frequencies, we will base our calculations on models of the differential number counts of ERSs at $\nu_\ast= 5\,\mathrm{GHz}$. The emission at the higher frequency range of interest will then be predicted using redshift-independent models for the radio spectra of the different source populations.

Differential number counts

We use the model for the number counts of sources from all radio populations (steep+flat+inverted) at $\nu_\ast= 5\,\mathrm{GHz}$ described in Toffolatti et al. (1998). We extrapolate this model linearly in log-log space towards large fluxes ($S_{\nu_\ast} > 10\,\mathrm{Jy}$), not covered in the original work. We combine this with a model for the joint number counts of flat and inverted sources, also at 5 GHz, from Tucci et al. (2011), described as a broken power law fitted to existing data. The abundance of steep sources is then given by the difference between these two combined populations:

The number counts of flat and inverted sources are then calculated by rescaling the flat+inverted distribution by the odds of a source being of either type given their spectral index. I.e.:

where:

The distribution of effective spectral indices $\alpha^5_1$ for the different populations is modelled, again following Tucci et al. (2011), as the sum of two truncated Gaussian distributions, one for the steep and one for the flat+inverted populations, with mean and variance dependent on the flux at 5 GHz $S_{\nu_\ast}$. (see Table 1 in Tucci et al. (2011) for the peak and sigma of these distributions in four flux ranges).

Freq. [GHz]	$\mu$K-arcmin baseline	$\mu$K-arcmin extended	$\ell_{\mathrm{knee}}$	$\alpha_{\mathrm{knee}}$	FWHM [′]
27	49.5	12.0	30	2.4	91
39	29.7	7.7	30	2.4	63
93	3.7	1.2	50	2.5	30
145	4.7	1.4	50	3.0	17
225	8.9	3.9	70	3.0	11
280	22.6	9.6	100	3.0	9

High-frequency spectra

To estimate source fluxes at the higher frequencies we use here for masking and analysis, we must extrapolate from the flux at the reference frequency $S_{\nu_\ast}$ (with $\nu_\ast= 5\,\mathrm{GHz}$) assuming a SED. The starting point for the three ERS populations considered here is the effective spectral index $\alpha^5_1$ in the low-frequency range (which, as we discussed above, is itself a stochastic spectral parameter). The SED of the different populations is then modelled as follows.

Steep-spectrum sources. The spectrum above 5 GHz is modelled as a power law with a high-frequency spectral index $\alpha_h=\alpha^5_1-\Delta\alpha_{\mathrm{steep}}$, where the steepening index $\Delta\alpha_{\mathrm{steep}}$ is itself a stochastic variable. Following Tucci et al. (2011), the distribution of $\Delta\alpha_{\mathrm{steep}}$ is Gaussian with mean 0.3 and standard deviation 0.2. The flux density at frequency $\nu$ is then

Flat-spectrum sources. We use the C0 model described in Tucci et al. (2011). As in the case of steep sources, the spectrum above $\nu_\ast=5\,\mathrm{GHz}$ is described as a power law with a steepened spectral index $\alpha_h=\alpha^5_1-\Delta\alpha_{\mathrm{flat}}$. The distribution of steepening indices is then modelled as a Gaussian with mean and standard deviations depending on the value of $\alpha^5_1$:

More complex models for the high-frequency SED of the flat-spectrum population were also proposed by Tucci et al. (2011). However, we adopt the simpler C0 as it represents the most conservative scenario, predicting the largest number of sources at high frequencies.

Inverted-spectrum sources. In this case, following Tucci et al. (2011), the high-frequency spectrum is modelled as a peaked spectrum. The spectrum rises up to a peak frequency $\nu_{\mathrm{peak}}$, and then falls with a post-peak spectral index:

with $\alpha_0=-0.1$. In this case, both $\nu_{\mathrm{peak}}$ and $\Delta\alpha_{\mathrm{inv}}$ are stochastic spectral parameters, with the distributions described in Tucci et al. (2011). $\nu_{\mathrm{peak}}$ follows a piecewise-constant distribution (see Table 3 in Tucci et al. (2011)), while the distribution of $\Delta\alpha_{\mathrm{inv}}$ is parametrised as:

$S_{\mathrm{lim}}$ [Jy]	$N_{\mathrm{PS}}$	$f_{\mathrm{sky}}$ [%]	$r$ ($\times10^3$) [true:0]	$A_{\mathrm{lens}}$ [true:0.3]	$\alpha_s$ [true:-0.8]	$\beta_s$ [true:-3]	$A_s$ [true:1]	$\beta_p$	$A_p$ ($\times10^{-4}$)	$\chi^2$
Baseline w/ nominal beams
ERS not modelled
0.01	25144	0.4	$-2.3 \pm11.4$	$0.282 \pm0.077$	$-0.02 \pm0.98$	$-3.03 \pm0.10$	$1.86 \pm1.57$	-	-	0.03
0.1	1699	7.8	$0.1 \pm2.9$	$0.303 \pm0.017$	$-0.60 \pm0.64$	$-3.03 \pm0.10$	$1.16 \pm0.57$	-	-	0.2
1	88	9.6	$-1.0 \pm3.0$	$0.352 \pm0.020$	$0.28 \pm0.50$	$-3.02 \pm0.10$	$1.91 \pm0.70$	-	-	4
ERS modelled as PL
0.01	25144	0.4	$-3.4 \pm11.7$	$0.354 \pm0.166$	$0.38 \pm1.07$	$-3.03 \pm0.10$	$3.40 \pm5.14$	$-2.38 \pm0.84$	$-6.4 \pm14.0$	0.09
0.1	1699	7.8	$0.0 \pm3.1$	$0.303 \pm0.042$	$-0.46 \pm0.80$	$-3.03 \pm0.10$	$1.27 \pm0.80$	$-2.28 \pm0.87$	$-0.3 \pm1.7$	0.1
1	88	9.6	$0.0 \pm2.9$	$0.310 \pm0.044$	$-0.45 \pm0.76$	$-3.02 \pm0.10$	$1.28 \pm0.69$	$-2.33 \pm0.78$	$1.8 \pm1.4$	0.16
Extended w/ nominal beams
ERS not modelled
0.01	21544	0.4	$0.5 \pm3.8$	$0.297 \pm0.022$	$-0.60 \pm0.79$	$-3.03 \pm0.10$	$0.99 \pm0.54$	-	-	0.1
0.1	1699	7.8	$-0.1 \pm0.9$	$0.306 \pm0.005$	$-0.62 \pm0.22$	$-3.01 \pm0.09$	$1.13 \pm0.26$	-	-	1
1	88	9.6	$0.0 \pm1.0$	$0.340 \pm0.006$	$1.02 \pm0.15$	$-2.91 \pm0.08$	$1.18 \pm0.27$	-	-	50
ERS modelled as PL
0.01	25144	0.4	$0.3 \pm3.9$	$0.304 \pm0.041$	$-0.27 \pm0.95$	$-3.02 \pm0.10$	$1.11 \pm0.68$	$-2.24 \pm0.85$	$-0.4 \pm1.4$	0.08
0.1	1699	7.8	$0.0 \pm0.9$	$0.299 \pm0.010$	$-0.78 \pm0.29$	$-3.01 \pm0.09$	$1.06 \pm0.27$	$-2.20 \pm0.83$	$0.2 \pm0.3$	0.2
1	88	9.6	$0.0 \pm0.9$	$0.301 \pm0.015$	$-0.78 \pm0.31$	$-3.02 \pm0.09$	$1.06 \pm0.28$	$-2.07 \pm0.25$	$2.3 \pm0.3$	0.4
Extended high resolution
ERS not modelled
0.01	21544	7.3	$-0.1 \pm0.5$	$0.301 \pm0.003$	$-0.72 \pm0.04$	$-3.00 \pm0.02$	$0.97 \pm0.04$	-	-	5
0.1	1699	9.5	$-0.2 \pm0.5$	$0.308 \pm0.003$	$-0.16 \pm0.04$	$-2.91 \pm0.02$	$0.70 \pm0.03$	-	-	237
1	88	9.7	$-1.0 \pm0.5$	$0.328 \pm0.003$	$1.41 \pm0.03$	$-2.36 \pm0.01$	$0.21 \pm0.01$	-	-	772
ERS modelled as PL
0.01	25144	7.3	$0.0 \pm0.6$	$0.301 \pm0.003$	$-0.79 \pm0.05$	$-3.00 \pm0.02$	$1.00 \pm0.04$	$-2.87 \pm0.60$	$0.0014 \pm0.0007$	0.4
0.1	1699	9.5	$-0.1 \pm0.5$	$0.303 \pm0.003$	$-0.79 \pm0.05$	$-3.00 \pm0.02$	$1.00 \pm0.03$	$-2.38 \pm0.06$	$0.0093 \pm0.0007$	1
1	88	9.7	$-0.1 \pm0.5$	$0.319 \pm0.004$	$-0.73 \pm0.07$	$-3.01 \pm0.02$	$0.99 \pm0.04$	$-2.22 \pm0.01$	$0.0649 \pm0.0013$	28

with $\sigma_\alpha=0.53$, and $\Delta\alpha_{\mathrm{inv}}$ constrained to the range $0\leq\Delta\alpha_{\mathrm{inv}}\leq2.1$. The distribution of $\alpha_h=\alpha_0-\Delta\alpha$ peaks at $\alpha_h\simeq-0.7$, and can be seen in Figure 9 of Tucci et al. (2011).

To summarise, the spectral parameters for the three populations studied here are:

As we have discussed, these spectral parameters are stochastic, and described by probability distributions which, moreover, are conditional on the population type ($t$) and, in some cases, on the flux at 5 GHz ($S_{\nu_\ast}$), $p(\theta_t\mid S_{\nu_\ast},t)$. The expectation values introduced in equations 3 and 6 are then given by:

where $Q$ is either $\Theta$ or $f_{\nu_1}S_{\nu_1}f_{\nu_2}S_{\nu_2}$. Note that these quantities are independent from the redshift, hence for ERSs equations 2 and 6 simplify by taking the expectation values outside the integral over redshift, which allows us to use the differential number counts already integrated over redshift, as presented in Sect. 2.3.1. We compute these expectation values by deterministic numerical quadrature on a multidimensional grid in the spectral-parameter space. For each spectral parameter $\theta^j$, we define a grid $\{\theta^j_i\}$ and the corresponding normalised quadrature weights

where $j$ denotes the spectral parameter type and $i$ the grid point. For parameters whose probability distribution does not depend on $S_{\nu_\ast}$, such as $\Delta\alpha_{\rm steep}$, the grid and its weights can be precomputed. For the others, such as $\alpha_1^5$, the corresponding quadrature weights are instead evaluated at each value of $S_{\nu_*}$. Finally, when the distribution of one spectral parameter depends on another spectral parameter (e.g. $\Delta\alpha_{\rm flat}$ and $\alpha_1^5$), its quadrature weights are evaluated conditionally at each grid point of the latter. Models with a higher-dimensional parameter space are therefore more computationally expensive to calculate.

Summary

To summarise, the number counts of masked ERSs and the polarisation cross-frequency power spectra of the residual point source contamination after masking are calculated using equations (2) and (6). The masking strategy, determining the unmasked fractions $f_\nu$ in Eq. (6), is described in § 2.2. The differential number counts for the three different source types considered here (steep, flat, and inverted) are described in § 2.4, while the source SEDs, and the probability distributions of their spectral parameters (required to estimate the ensemble averages in equations (3) and (6)) are described in § 2.5. We compute only the Poisson term of the ERS power spectrum as the clustering term is negligible compared to the former ((Toffolatti et al. 1998; Gonzalez-Nuevo et al. 2005)) unless aggressive flux cuts, $S_{\rm lim} \ll10$ mJy, are considered, which is not the case in our analysis.

The final ingredient required to predict the level of residual point source contamination from polarised ERSs is the squared-mean polarisation fraction $\langle\Pi^2\rangle$. Here we follow the approach of Lagache et al. (2020), assuming a constant polarisation fraction across all sources with $\langle\Pi^2\rangle^{1/2} = \langle\Pi\rangle= 2.8\%$, in agreement with Planck, SPT and ACT data ((Puglisi et al. 2018; Gupta et al. 2019; Datta et al. 2019)). In § 4.6.1 we will consider the impact of scatter in the distribution of polarisation fractions, and instead use a larger value $\langle\Pi^2\rangle^{1/2} = 4.94\%$ (see Puglisi et al. (2018)) for the flat-spectrum population, which dominates the residual point source contamination.

The software developed to compute the number counts of masked ERSs, the polarisation cross-frequency power spectra of the residual contamination after masking and the loss of sky area due to PS masking is made publicly available¹.

Model for polarised dusty star-forming galaxies sources

Dust grains in the interstellar medium absorb optical and ultra-violet radiation and emit radiatively in the far infrared and microwave spectrum with what is typically well described by a modified black body spectrum. Because of their elongated shape, dust grains may partially align with the interstellar magnetic field, and thermal dust emission in dusty star-forming galaxies is therefore expected to show a certain degree of polarisation. The signal we receive from one DSFG, however, is integrated over the whole galaxy, so its observed fractional polarisation is expected to be much lower than that coming from any one point in the galaxy. Few direct measurements exist: for individual objects, Greaves & Holland (2002) measured a global polarisation degree of only 0.4% at 353 GHz for M82, while Seiffert et al. (2007) set a 99% confidence upper limit of 1.54% for Arp 220 at the same wavelength. Stacking $\sim4700$ DSFGs in the Planck maps, Bonavera et al. (2017) reported mean fractional polarisations of $3.10 \pm0.75\%$ and $3.65 \pm0.66\%$ at 217 and 353 GHz, respectively. A re-analysis of the same data conducted by Trombetti et al. (2018) showed no significant detection and derived 90% confidence upper limits of $\langle\Pi\rangle< 2.2\%$ at 353 GHz and $< 3.9\%$ at 217 GHz. On the theoretical side, De Zotti et al. (2018) estimated a mean polarisation degree of 1.4% for spiral galaxies. We hence adopt this latter value and, as in the case of ERSs, we do not model polarised number counts directly, but work with the intensity counts and convert predictions for the residual cross-spectrum in intensity to polarisation through this polarisation fraction.

In intensity, the spectrum of a DSFG is dominated by a modified black body spectrum, peaking in the far infrared ($T \sim20$ K) and falling as $\nu^{2+\beta_{\rm dust}}$ towards the millimetre, with $\beta_{\rm dust}$ the dust emissivity spectral index. The frequency range of a typical CMB experiment (20 GHz $\lesssim\nu\lesssim$ 300 GHz) samples this steeply rising part of the spectrum, where the DSFG contribution grows rapidly with frequency and may become relevant only towards the top of the band range.

We follow the two-population model for the Cosmic Infrared Background (CIB) of Zagatti et al. (2024), which we summarise in here, and extend to polarised point sources. The model is built on the galaxy evolution framework of Granato et al. (2004) and Cai et al. (2013), a hybrid approach in which massive proto-spheroidal galaxies are treated with a physical, semi-analytic prescription, while late-type and starburst galaxies are described through an empirical approach. Within this framework, DSFGs are split in two samples which occupy dark matter haloes differently and emit at different epochs:

• the early-type (ET) galaxies, are the massive proto-spheroids, which dominate the emissivity at high redshifts, $z \gtrsim1.5$, and reside in the more massive haloes;

• the late-type (LT) galaxies, are modelled as the sum of the cold (normal late-type) and warm (starburst) components, dominate at smaller redshifts, $z \lesssim1.5$, and populate less massive haloes.

The two populations entering the analysis are described through different SEDs and are expected to live in different haloes, therefore they require separate halo occupation distributions (HODs). The differential number counts are computed at each frequency band from the luminosity functions and the SEDs of each population.

We note a difference with respect to the ERS case, for which the number counts at a reference frequency are the starting point. From an observational point of view, the millimetre counts of dusty galaxies are less determined than the radio counts at 5 GHz, while the infrared luminosity function is tightly constrained over a wide range of redshift by Herschel (Berta et al. 2011; Béthermin et al. 2012), Spitzer (Béthermin et al. 2010), and Planck (Negrello et al. 2013). Starting from the luminosity function allows for a straightforward computation of the number counts, without the need for a spectral extrapolation between a reference band and the observational bands. Therefore, for DSFG the reference frequency $\nu_*$ coincides with each measurement frequency $\nu$ (defined in Sect. 2.1), unlike in the ERS case.

Differential number counts from luminosity functions and SEDs

For a population $t$ with a bolometric infrared luminosity function $\Phi_t(\log_{10} L_{\rm IR},z)$, the surface density of sources per unit flux and redshift at the observed frequency $\nu$ follows from Eqs. (8) and (12) of Cai et al. (2013):

where $d^2V/dz\,d\Omega$ is the comoving volume element per unit solid angle, and the luminosity corresponding to a source of flux $S_\nu$ at redshift $z$ is

Here $L_{\mathrm{IR}}$ is the bolometric infrared luminosity, $D_L(z)$ is the luminosity distance, and $\varphi_t(\nu)$ is the SED template of population $t$ normalised to $\int\varphi_t(\nu)\,\mathrm{d}\nu= 1$. We evaluate Eq. (21) monochromatically at the nominal band centres and adopt best-fit Planck cosmological parameters to relate distances and redshifts, consistently with the halo model of Sect. 2.4.3.

In the notation of Sect. 2.1, for DSFG there is no spectral parameter dependence, but, contrary to the ERS case, the redshift dependence can not be neglected: the model assigns one deterministic SED to all sources of a given population at a given redshift $z$, so that the fluxes of a source at any two frequencies $\nu$ and $\nu_i$ entering the calculations are related by

The average $\langle\cdot\rangle_{\boldsymbol{\theta}_t}$ appearing in equations 3 and 6 therefore involves no distribution of spectral parameters, and reduces to a deterministic factor, that remains inside the integral over redshift weighted by the counts.

Luminosity functions and SEDs of the two populations

Both populations are described by the bolometric infrared luminosity function of Cai et al. (2013), as implemented by Zagatti et al. (2024), to which we refer for details. For the LT galaxies, these are modified log-Gaussian functions following Saunders et al. (1990), whose characteristic density and luminosity evolve with redshift up to $z \simeq1$ and are then smoothly cut off. For the ET galaxies the luminosity function comes from the physical model of Granato et al. (2004), and is truncated below $z \simeq1$ (see Figure 1).

Each sub-population carries a single SED template: for the proto-spheroids, that of a strongly lensed $z \simeq2.3$ submillimetre galaxy (Swinbank et al. 2010; Lapi et al. 2011); for the cold and warm late types, templates built from Dale & Helou (2002) and Smith et al. (2012) at $\log_{10}(L_{\mathrm{IR}}/L_\odot)=9.75$ and $11.25$ respectively. At lower rest frequencies the templates also include the radio emission associated with star formation, namely synchrotron radiation from supernova remnants and free--free emission from Hii regions. This falling component, together with the rising dust one, produces a minimum in $S_\nu$ at a rest frequency of about 68 GHz. We stress that this emission is physically distinct from that of the ERS of Sect. 2.3, which originates in the jets and lobes powered by an active nucleus, so that adding the two models does not double count it. Nevertheless, assuming that the radio components of the DSFGs should be scaled by the same polarised fraction $\Pi$ as their thermal dust component is likely incorrect (certainly so in the case of the free-free contribution). As we will see, the impact of polarised DSFGs is significantly smaller than that of ERSs and thus, the choice to scale all the DSFG components equally is a reasonable approximation.

Since each sub-population carries one template with no scatter in the SED, Eq. (23) is fixed for a given $(t,z)$. The cross-frequency correlation coefficients predicted by the model are therefore upper limits: introducing a scatter in the SEDs would reduce the coherence between bands.

Power spectra and clustering

Unlike ERSs, DSFGs trace the large-scale structure coherently and contribute a potentially non-negligible clustering term in addition to the Poisson contribution. We therefore write the total DSFG power spectrum as

where the first (Poisson) term is Eq. (6) evaluated with DSFG counts, and the second is computed with the two-population halo model of Zagatti et al. (2024).

The quantity that carries the frequency dependence in the model is the emissivity (Figure 1), $j_t(\nu,z)\equiv\mathrm{d}I_\nu/\mathrm{d}z$, the contribution to the specific intensity of the unresolved background per unit redshift, in Jy sr$^{-1}$ per unit $z$, so that $\int j_t(\nu,z)\,\mathrm{d}z$ is the total intensity from population $t$. We evaluate it separately for each population from the counts of Eq. (21):

with $t\in\{\mathrm{ET},\mathrm{LT}\}$. Here $f_\nu$ is the masking factor of Sect. 2.2. It is worth noting that, in this case, the emissivity is dominated by faint sources, and is therefore essentially unaffected by the masking procedure. Most of the dependence on the masking strategy resides in the Poisson term.

Projecting the emissivities with the Limber approximation gives the angular power spectrum of the clustering part,

with $k_\ell\equiv(\ell+ 1/2)/\chi$ in the Limber approximation (Limber 1953; LoVerde & Afshordi 2008), and $P_{\mathrm{gal}}^{tt'}(k,z)$ the halo model galaxy power spectrum between sources of type $t$ and $t'$, including both the one-halo and two-halo contributions, built from the HODs of the two populations. We assigned different minimum hosting halo masses, $M_{\min}=10^{12}\,M_\odot h^{-1}$ for ET and $10^{11}\,M_\odot h^{-1}$ for LT, reflecting the fact that ET galaxies form in more massive structures.

Converting the intensity clustering into polarisation requires more than the $\langle\Pi^2\rangle/2$ factor alone used for the Poisson term, since we cannot assume random polarisation angles of physically associated sources. Following Lagache et al. (2020) we write

where $\rho_{1h}$ and $\rho_{2h}$ quantify the degree of alignment of the polarisation angles within a halo and between distinct haloes, respectively. We set $\rho_{2h}=0$: galaxies in different haloes have no physical reason to share a preferred orientation, so the two-halo term does not contribute to the clustering polarisation power. We set $\rho_{1h}=1$, that is, we assume perfect alignment within a halo. This is an upper bound rather than a prediction.

Summary of quantities computed by the model

The model used for DSFGs takes as inputs

where $\mathrm{HOD}_{\mathrm{ET/LT}}$ denotes the HOD model used to estimate the clustering of ET or LT galaxies. We use these ingredients to calculate the differential number counts of DSFGs (Eq. (21)), and from that the two observables defined in Sect. 2.1, needed for this analysis:

• The angular number density of DSFGs whose flux at the masking frequency exceeds $S_{\mathrm{lim}}$. In the notation used in this section, this reads

with the total obtained as $N_{\nu_m}=\sum_t N_{\nu_m}^{t}$.

• The $BB$ cross-spectra, now with both a Poisson and a clustering contribution, as in Eq. (26). Both terms are built from the same counts, but enter with a different power of the flux. The Poisson term weights each source by the product of its fluxes in the two bands,

with $S_{\nu_2}$ and $S_{\nu_1}$ obtained from $S_\nu$ from Eq. (23). The clustering contribution is given by Eq. (24) scaled by $\langle\Pi^2\rangle/2$, with $\rho_{1h}=1$ and $\rho_{2h}=0$.

We use the polarisation fraction of De Zotti et al. (2018), used also by Lagache et al. (2020), $\langle\Pi_{\mathrm{IR}}\rangle=1.4\%$, independent of frequency and flux, and approximate $\langle\Pi^2\rangle$ by $\langle\Pi\rangle^2$. As for the ERS, the results are converted from $\mathrm{Jy\,sr}^{-1}$ to $\mu\mathrm{K}_{\mathrm{CMB}}^2$ with the monochromatic conversion factors evaluated at band centres.

Analysis framework

We propagate the impact of point source contamination by generating synthetic cross-frequency power spectrum data containing CMB and Galactic foreground contributions as well as the residual PS contribution calculated as described in § 2. We then derive constraints on cosmological and foreground parameters through a power spectrum-based component separation and likelihood framework. The impact of point sources is then reflected in the bias of the inferred parameters (caused by the residual PS contamination), and in their statistical uncertainties (caused largely by the loss of sky area due to masking).

Synthetic data

We build synthetic data vectors of cross-frequency polarisation power spectra. The models used to simulate the CMB and Galactic components are as follows:

(i) The CMB signal is generated with CAMB (Lewis et al. 2000) for a $\Lambda$CDM model with the best-fit cosmological parameters from Planck. We fix the tensor-to-scalar ratio to $r=0$ (no primordial $B$ modes), and rescale the $B$-mode power spectrum due to gravitational lensing by an amplitude $A_{\mathrm{lens}}=0.3$. This represents the level of delensing (70%) that could be achievable by an experiment such as SO (Hertig et al. 2024).

(ii) The dust and synchrotron components are modelled as having a power-law power spectrum of the form

with $\ell_0=80$. This is the power spectrum for component $x$ (with $x=d$ for dust and $x=s$ for synchrotron) at a pivot frequency $\nu_0^x$. We use $\nu_0^s=23\ \mathrm{GHz}$ and $\nu_0^d=353\ \mathrm{GHz}$. These power spectra are scaled to other frequencies assuming a modified black-body spectrum for dust and a power-law spectrum for synchrotron:

Note that both of these spectra are defined in Rayleigh-Jeans temperature units. We use the following values for the foreground parameters defining these contributions:

These values are representative of the amplitude, scale dependence, and frequency spectra of both components in the d0 and s0 models implemented in PySM (Thorne et al. 2017) within the angular footprint of an SO-like (i.e. Chile-based) experiment.

We note that the true Galactic sky is likely more complex than this, with each foreground component not exactly described by a power law $\mathcal{D}_{\ell}$ model and more complex spectral dependence, with foreground SEDs varying across the Galaxy, both across the sky and along any given line-of-sight (e.g. Krachmalnicoff et al. 2018b; Planck Collaboration [61]; Martire et al. 2022; Weiland et al. 2022). However, considering such a simple model allows us to explore an optimistic scenario recovering the tightest possible constraints on $r$ as well as all other model parameters, making the results more sensitive to the impact of residual point source contamination. The impact and performance of more complex foreground parametrisations on $B$-mode constraints have been studied in the literature (e.g. Azzoni et al. 2021; Vacher et al. 2022; Wolz et al. 2024; Liu et al. 2025).

Error model

The synthetic data described in the previous section are generated without any measurement noise component (neither due to cosmic/sample variance nor instrumental noise). This ensures that any deviation in the inferred model parameters from the input values used to generate the synthetic data represents a genuine bias caused by the presence of point source contamination. Nevertheless, in order to obtain representative statistical uncertainties on the model parameters, against which any such biases may be compared, we calculate the covariance matrix of the synthetic data accounting for all sources of measurement uncertainty, including the impact of point source masking.

The power spectrum covariance is calculated assuming all fields involved are Gaussian, as is standard in most CMB power spectrum analyses ((Hamimeche & Lewis 2008; Aghanim et al. 2016; Atkins et al. 2025)).² To account for the impact of mode-coupling due to the presence of a finite sky mask we use the improved Narrow Kernel Approximation ((García-García et al. 2019; Nicola et al. 2021)) as implemented in NaMaster³ ((Alonso et al. 2019)).

As a mask representative of the sky region targeted by current ground-based experiments, we use the publicly-available hits map of the SO Small Aperture Telescopes (SATs, The Simons Observatory Collaboration (2019)), shown in Figure 2. To account for the impact of PS masking in the covariance matrix, we re-scale the covariance matrix obtained in the absence of any PS mask by sky area remaining after masking:

Here $f_{\mathrm{sky,no\,PS}}$ is the effective area of the original mask ($f_{\mathrm{sky,no\,PS}} = 0.097$ in the case of the mask shown in Figure 2). The available sky fraction after point source masking is given by

where $N_{\nu_m}(>S_{\mathrm{lim}})$ is the angular number density of masked point sources (see Eq. (2)), and $\Omega_{\mathrm{mask}} = \pi R_{\mathrm{mask}}^2$ is the solid angle associated with the mask of a single point source (with $R_{\mathrm{mask}}$ given in Eq. (7)).

We include the contribution from instrumental noise in the covariance calculation. For this, we use the noise model presented in The Simons Observatory Collaboration (2019), with

Here $N_{\mathrm{white}}$ is the white noise amplitude, and $\ell_{\mathrm{knee}}$ and $\alpha_{\mathrm{knee}}$ characterise the non-white (“$1/\ell$”) component, largely associated with the atmosphere. These three parameters vary across frequency channels. We will consider an SO-like experiment, with 6 frequency bands listed in Table 1, together with the values of $N_{\mathrm{white}}$, $\ell_{\mathrm{knee}}$, and $\alpha_{\mathrm{knee}}$ considered in this work, as well as the beam FWHM apertures. To model the $1/\ell$ noise component we use the values of $\ell_{\mathrm{knee}}$ and $\alpha_{\mathrm{knee}}$ corresponding to the “pessimistic” model of The Simons Observatory Collaboration (2019). We investigate two white noise scenarios, the “baseline” noise level of The Simons Observatory Collaboration (2019), corresponding to a 5-year experiment, and the “extended” configuration of The Simons Observatory Collaboration (2026), with up to $\sim3$ times higher sensitivity. The resulting noise spectra, in both configurations, are shown in Figure 3, including the impact of deconvolving the instrumental beam.

When considering ERSs we will assume they are masked in the lowest frequency channel (i.e. $\nu_m = 27$ GHz), where the beam FWHM is $91^{\prime}$, as this is the channel with more and brighter sources. Masking radio sources may thus lead to a significant loss of area. A different choice could be made for the DSFGs, however we do not explore this as we find that no source masking is required for DSFGs in our analysis.

Component separation and likelihood

We use a power spectrum-based Bayesian framework to separate the CMB and foreground components and derive constraints on the tensor-to-scalar ratio and other model parameters. The methodology is described in detail in e.g. Wolz et al. (2024) (see Pipeline A), and we only outline its main elements here. Cross-frequency $B$-mode power spectra are modelled as a linear combination of contributions from different components:

The components considered include the CMB $B$-mode signal, Galactic thermal dust and synchrotron. The Galactic foregrounds are modelled as described in § 3.1 in terms of both scale and frequency dependence. The CMB contribution is described in terms of $r$ and $A_{\mathrm{lens}}$, which multiply tensor-only and lensing-only templates respectively, calculated fixing all other $\Lambda$CDM cosmological parameters to the Planck best-fit values. The free parameters of the fiducial model are then

In addition to this minimal model, we will also consider the possibility of absorbing any residual point source contamination in terms of an additional component, labelled $p$, with an effective SED modelled as a power law in Rayleigh-Jeans temperature units $f_p(\nu) = (\nu/\nu_0)^{\beta_p}$ and with a white angular power spectrum characterised by a free amplitude $A_p$ (i.e. $C_{\ell}^{p} = A_p\ \forall\ell$). Finally, we have also tested the interplay between residual extragalactic PS contamination and additional Galactic foreground complexity, by allowing for frequency decorrelation as modelled in BICEP2 Collaboration et al. (2018) and checking that such a term is not detected.

We approximate the power spectrum-level likelihood as a multivariate Gaussian distribution, which has been shown to be a good approximation on the scales explored here (Wolz et al. (2024)). We adopt broad, largely uninformative priors on all model parameters, although we enforce the positivity of all foreground amplitude parameters ($A_{\mathrm{d}}$, $A_{\mathrm{s}}$, $A_p$). The model parameters are then recovered via Metropolis-Hastings sampling using Monte Carlo Markov Chains (MCMC). This multi-frequency power spectrum approach has been extensively used in CMB analysis, e.g. Planck Collaboration [60]; Ade et al. (2021); Louis et al. (2025), as well as in forecast papers, e.g. Wolz et al. (2024); Azzoni et al. (2021). We use the specific implementation of this method available in BBPower⁴.

The synthetic data vector described in § 3.1 consists of all auto- and cross-correlations between the 6 frequency channels listed in Table 1, sampled in multipole bandpowers of width $\Delta\ell= 10$ in the range $32 < \ell\le600$.

Results

We first present results regarding the impact of residual contamination from polarised ERSs, as the impact of DSFG contamination on $B$ modes is significantly smaller. We present results for DSFGs in § 4.7.

ERS forecasts

Number counts

Using Eq. (2), we compute the ERS angular number density above the flux cut $S_{\mathrm{lim}}$, $N_{\nu_m}^{T}(> S_{\mathrm{lim}})$, multiplied by the SO SAT sky area, at the different SAT frequencies $\nu$. The result is shown in Figure 4. As expected, the number counts are largely dominated by the flat-spectrum sources. Moreover, for a fixed $S_{\mathrm{lim}}$ the number of sources increases towards lower frequency, motivating the choice of $\nu_m = 27$ GHz. We note that also the inverted spectrum sources follow this behaviour, since their number count is dominated by sources whose spectra peak before the SO lowest frequency of 27 GHz.

In our analysis, these number counts are relevant for determining the sources to mask. We choose to identify the ERSs at the instrument lower frequency of 27 GHz (common both to the SATs and the LAT instruments). In principle, ERSs could be identified at even lower frequencies leveraging external data. However, we do not consider this necessary, as extending the source identification to lower frequency would increase the uncertainty associated with extrapolating the source frequency scaling over broad frequency ranges. Moreover, the choice of $\nu_m$ is not determined solely by maximising the number counts, but also by choosing at what frequency to prioritize the ERS masking, as the brightest sources at a given frequency are possibly not the same at other frequencies. We hence choose to prioritize the masking at the lowest frequency band, as that is the one where we have higher ERS signal (but not necessary higher point source signal-to-noise ratio, as we show later). Finally, the highest-frequency bands are not well suited for identifying the ERS to be masked, since thermal dust emission becomes increasingly important at these frequencies and could dominate the observed flux.

Angular power spectra

Having chosen the masking frequency $\nu_m$, we can then compute the amplitude of the residual ERS power spectrum after masking with different flux cuts $S_{\mathrm{lim}}$. Figure 5 shows the amplitude of the ERS signal when masking with $\nu_m = 27$ GHz, and $f_{\mathrm{FWHM}} = 3$, which means each masked PS is removed from the map almost completely (less than 0.2% of the flux remains). This conservative approach may be accurate if we are only interested in suppressing the residual PS amplitudes, but may not be the best trade-off when also accounting for the noise and the fraction of sky lost due to masking (see § 4.3).

Residual ERS contamination and component separation

Per-band masking

We start by considering the flux limits and number counts studied in Lagache et al. (2020). This allows us to get a sense of what flux limits may be reasonably expected in future experiments, as well as a rough estimate of the parameter biases caused by different levels of PS contamination. Note that, although Lagache et al. (2020) assumed that each frequency channel would be masked according to the sources detected in it by the instrument being analysed, we will show results for a uniform mask applied across frequencies, built from the sources detected in a single channel for its corresponding flux limit. The flux limits estimated by Lagache et al. (2020) for the SO SAT and Large Aperture Telescope (LAT) bands are shown in Table 2, together with the associated expected number counts. In particular, at 27 GHz they find $S_{\mathrm{lim}}^{\mathrm{SAT}\ 27\ \mathrm{GHz}} = 1.983$ Jy, growing at 39 GHz due to the decaying nature of radio spectra, but decreasing afterwards due to the better sensitivity and smaller beam sizes.

$\nu$ [GHz]	SAT $S_{\mathrm{lim}}$ [Jy]	SAT $N_{\nu}(>S_{\mathrm{lim}})$	LAT $S_{\mathrm{lim}}$ [Jy]	LAT $N_{\nu}(>S_{\mathrm{lim}})$
27	1.983	27	0.01502	1.565 \times 10^{4}
39	4.177	6	0.00994	$2.004 \times10^{4}$
93	3.385	8	0.00296	5.028 \times 10^{4}
145	1.593	31	0.00261	$4.966 \times10^{4}$
225	0.7558	94	0.00568	$2.000 \times10^{4}$
280	0.4351	191	0.01011	$1.067 \times10^{4}$

Adopting the flux limit found for the SATs at 27 GHz, and propagating the impact of the residual contamination to the model parameters, we find the results shown by the black contours of Figure 6. We find that the synchrotron parameters, particularly the tilt $\alpha_s$ and amplitude $A_s$, as well as the CMB lensing amplitude $A_{\mathrm{lens}}$, are significantly biased. This is not surprising, as synchrotron is the Galactic foreground that most closely matches the SEDs of radio sources, while the scale dependence of the lensing $B$-mode component is largely flat on the scales ($\ell< 600$) studied here, similar to Poisson sources. Fortunately, the tensor-to-scalar ratio $r$ remains unbiased. Even at the more stringent flux limit found at 280 GHz (even though masking at this high frequency is not realistic given the associated emission from dust in the same sources, disregarded here), we find that these parameter biases do not disappear. This analysis is performed with $f_{\mathrm{FWHM}} = 1$, but because of the high $S_{\mathrm{lim}}^{27\ \mathrm{GHz}}$ the result does not change with $f_{\mathrm{FWHM}} = 3$.

We may also consider the case where the positions of sources to be masked come not from the same instrument being analysed, but from an external dataset, more sensitive to point source emission (e.g. because of better resolution). In particular, let us use the detection limit for the SO LAT at 27 GHz, more than 100 times smaller than the SAT flux limit at the same frequency. At this limit we are able to detect and mask a much larger number of sources, reducing the parameter bias but also the size of the sky area available for the analysis. We find that the parameter biases reported earlier disappear completely, albeit at the cost of a significant increase in the posterior parameter uncertainties (by a factor 1.8 in the case of $r$). While the discussion in this section has focused on the baseline configuration similar results hold with the higher sensitivity of the extended configuration.

For the remainder of this paper we will therefore explore a range of flux limits, spanning the fluxes available to SAT-like and LAT-like experiments, to study whether a suitable trade-off between these two effects (reduced bias and increased variance) can be achieved. Specifically, we will study three cases: $S_{\mathrm{lim}} = 0.01$ Jy, 0.1 Jy and 1 Jy.

Choice of the masking radius for each ERS

The choice of masking strategy depends not only on the flux limit used to detect and mask sources, but also on the aperture size used for masking. In our analysis we assume that sources are masked at the lowest frequency $\nu_m = 27$ GHz, where $\theta_{\mathrm{FWHM}}(\nu_m) = 91'$. The choice of masking fraction $f_{\mathrm{FWHM}}$ (which determines the masking radius – Eq. (7)) in this case is not trivial: large values of $f_{\mathrm{FWHM}}$ ensure that sources are fully masked at the masking frequency $\nu_m$, and certainly so at all other (higher) frequencies, where the beam is smaller, although this is at the cost of an increased loss of usable area. Choosing a smaller $f_{\mathrm{FWHM}}$ thus reduces this loss, increasing the residual contamination at $\nu_m$ but perhaps not so much at the higher frequencies (e.g. 90-150 GHz), where the beam is smaller and the overall sensitivity is typically higher. We explore different choices of $f_{\mathrm{FWHM}}$:

(i) $R_{\mathrm{mask}}(f_{\mathrm{FWHM}} = 0.5) = 22.75'$,

(ii) $R_{\mathrm{mask}}(f_{\mathrm{FWHM}} = 1) = 45.5'$,

(iii) $R_{\mathrm{mask}}(f_{\mathrm{FWHM}} = 2) = 91'$,

(iv) $R_{\mathrm{mask}}(f_{\mathrm{FWHM}} = 3) = 136.5'$

We study the trade-off between area loss and residual PS contamination of these choices in Figure 7. We find that $f_{\mathrm{FWHM}} = 1$ (corresponding to a circular mask the size of the beam FWHM at $\nu_m = 27$ GHz) achieves a somewhat reasonable balance: this choice allows us to mask up to 1000 sources across the sky without an excessive loss of sky area and, at the same time, although half of the masked PS flux remains at $\nu_m$, it vanishes almost completely at the more sensitive 90 and 150 GHz bands. Other choices lead to either a significantly larger sky loss or a non-negligible contamination at intermediate frequencies.

To quantify the impact of this choice more rigorously, we study the parameter constraints obtained with each of the options above. We assume the extended configuration in this case. The result is shown in Figure 8. For simplicity we adopt a flux cut of $S_{\mathrm{lim}} = 0.1$ Jy in all cases but $f_{\mathrm{FWHM}} = 0.5$, where a more stringent cut of $S_{\mathrm{lim}} = 0.01$ may be adopted with only a mild loss of area. We find that the fiducial case of $f_{\mathrm{FWHM}} = 1$ achieves the best balance of all explored: it achieves the smallest statistical uncertainties while incurring virtually no bias on the model parameters. The mild $\sim1\sigma$ bias on $A_{\mathrm{lens}}$ (which we will return to in Section 4.2.4) is still apparent in the more conservative cases $f_{\mathrm{FWHM}} = 2$ and 3, showing that it is driven by the residual contribution from unmasked sources, and not by the residual flux of the masked ones. Interestingly, even with the more stringent flux cut, the case $f_{\mathrm{FWHM}} = 0.5$ is significantly more biased, showing that, in this case, the dominant effect is the residual flux of the masked sources which, as we showed in Figure 7, extends to the high-sensitivity 90 GHz channel.

Masking: bias-uncertainty trade-off

We show in Figure 9 residual contribution to multi-frequency power spectra from the residual ERS contamination as a function of $S_{\mathrm{lim}}$ for $f_{\mathrm{FWHM}} = 1$. As noted in the previous section, the two lowest-frequency channels are the most affected by the residual flux of masked sources, manifested here in the plateau found for $S_{\mathrm{lim}} \lesssim0.1$ Jy in the case of the auto- and cross-correlations of both of these channels. To better understand the impact of PS contamination on the data, Figure 10 shows the angular auto-power spectra at $\nu= 27$ GHz and 93 GHz, including all contributions (CMB, Galactic foregrounds, and residual PS contamination) for three flux cuts: $S_{\mathrm{lim}} = 1$ Jy, 0.1 Jy, and 0.01 Jy. The error power spectra for both frequencies is also shown for comparison. Owing to the better resolution and white-noise amplitude, the 93 GHz achieves a much better sensitivity across all frequencies, and is comparatively more strongly affected by the presence of point sources, particularly on scales $\ell\gtrsim150$. This reinforces the point made in the previous section, when choosing a value of $f_{\mathrm{FWHM}}$ that ensures point sources are fully masked at the more sensitive central channels. The amplitude of the spectra are visibly enhanced above $\ell\simeq150$ (more evidently so for the higher flux cut), following a $D_\ell\propto\ell^2$ behaviour characteristic of both lensing $B$ modes (at these scales) and point sources. On larger scales the power spectrum is dominated by Galactic foregrounds with a red scale dependence.

Figure 11 shows the parameter constraints obtained for these flux limits, for the extended noise configuration. Comparing the cases $S_{\lim}=0.01\,\mathrm{Jy}$ and $0.1\,\mathrm{Jy}$ we find that the former leads to an excessive increase in statistical uncertainties, with the error on $r$ growing by $\sim280\%$. On the other hand, $S_{\lim}=1\,\mathrm{Jy}$ introduces very significant biases on the synchrotron and lensing parameters. In particular $\alpha_s$ is biased high by $12\sigma$ ($2\sigma$ for the baseline configuration), recovering a growing “synchrotron” spectrum (in fact corresponding to the PS contamination). Interestingly, $\beta_s$ also experiences a milder $\sim1\sigma$ upwards bias, shifting the synchrotron SED to the larger spectral indices found in extragalactic sources on average. Finally, $A_{\mathrm{lens}}$ is biased high by nearly $7\sigma$ ($3\sigma$ for baseline). These large biases disappear almost completely when adopting a more conservative $S_{\lim}=0.1\,\mathrm{Jy}$. $A_{\mathrm{lens}}$ and $\alpha_s$ are now only $\sim1\sigma$ high ($\sim0.3\sigma$ for baseline). In turn, the parameter errors grow only modestly: 20% for $A_{\mathrm{lens}}$, 40% for $\alpha_s$, and around 10% or less for all the other parameters. We summarise the relevant information for these three cases and both configurations in Table 3. The table also reports the $\chi^2$ values associated with the best-fit model in each case. Since the synthetic data was generated without statistic noise, these values should be understood as “excess” $\chi^2$ due to the presence of PS contamination. In all cases explored in this section, even if the $\chi^2$ increases for larger values of $S_{\lim}$, the degradation to the goodness of fit is still rather small, and one would not suspect the presence of ERS contamination purely on the basis of this metric.

Marginalising over residual PS contamination

Although a PS flux limit of $S_{\lim}=0.1\,\mathrm{Jy}$ provides a reasonable balance between bias and variance, a small $\sim1\sigma$ bias in $A_{\mathrm{lens}}$ and $\alpha_s$ remains, which could affect the interpretation of $B$-mode analyses with future datasets. In addition to this, current uncertainties in the modelling of polarised point sources at 20-200 GHz frequencies are significant, and could make this mild biases more significant. It is therefore interesting to explore whether the residual PS contamination may be effectively marginalised over, and what the penalty of doing so is in terms of posterior parameter constraints.

To explore this, we introduce a fourth sky component incorporating this contamination. We model the effective SED of this component as a power law in Rayleigh-Jeans unit, with a spectral index $\beta_p$ ($f_p(\nu)=(\nu/\nu_0^p)^{\beta_p}$), and with a flat/white power spectrum, characterised by an amplitude $A_p$. In detail, this can only be an approximate model for the residual contamination, not only because the combination of emission from different ERS samples is not described by an exact power law, but also because the distinct effective masking of sources at different frequencies will cause further frequency decorrelation. Nevertheless, this simplistic model may be enough to capture the bulk of the residual contamination, reducing parameter biases below their statistical uncertainties.

The result of adding this component on the inferred model parameters is shown in Table 3 for the different configurations explored here, and in Figure 12 in the case $S_{\lim}=1\,\mathrm{Jy}$ for the extended noise configuration. Overall we find that including the additional component resolves the parameter biases described in the previous section in almost all cases, even those that exhibit a significant residual ERS signal. For instance, where we found biases of up to $7\sigma$ in $A_s$, and $12\sigma$ in $\alpha_s$ for $S_{\lim}=1\,\mathrm{Jy}$, both parameters are recovered with no apparent bias when including the additional component. A penalty is paid in the form of a $\sim2$-fold increase in the posterior uncertainties in both parameters, although the constraints on $r$ remain largely unchanged. From the constraints on the amplitude $A_p$, we find that the PS residual is detected at the $\sim7.5\sigma$ level. In cases where the residual contamination is less significant, either due to using a more stringent flux cut, or to the larger noise in the baseline configuration, a mild bias in $\alpha_s$, below the $1\sigma$ level, remains. This is largely due to the three-way degeneracy between $A_p$, $\beta_p$, and $\alpha_s$, leading to volume/projection effects in the posterior constraints. All other parameters, including $A_{\mathrm{lens}}$, remain unbiased, however. An important point in this is the choice of pivot frequency $\nu_0^p$ for this component: when the residual PS contribution is close to the threshold of detection, fitting it with a non-linear two-parameter model $(A_p,\alpha_p)$ can lead to strong projection effects when $\nu_0^p$ is far from the channel with the highest SNR. Although these may be avoided by imposing a Jeffreys prior on the non-linear parameter $\alpha_p$, we find that these effects largely disappear when we choose $\nu_0^p=93\,\mathrm{GHz}$,

our most sensitive channel. With this choice, we show in Figure 13 that we are able to accurately reconstruct the true residual ERS signal amplitude at 93 GHz, at the price of measuring a slightly biased amplitude at the lower frequencies. In other words, setting the pivot frequency at 93 GHz we can correctly measure $A_p$, but we can not constrain $\alpha_p$ when there is insufficient signal at the other frequencies. In summary, we conclude that this simple effective model (a Poisson-like component with power-law SED) is sufficient to capture the residual point source contamination, even with minimal masking ($S_{\mathrm{lim}} = 1$ Jy), and recover unbiased constraints on cosmological and foreground parameters, without the need to fit for individual white noise amplitudes in each cross-frequency spectrum.

High angular resolution

The impact of point sources in the cases studied so far is mitigated by the fact that the ERS signal, peaking at low frequencies, dominates over DSFGs (see § 4.7). At these low frequencies, SAT-like instruments have poor sensitivity on small scales, where PS contamination may dominate over Galactic foregrounds. It is therefore interesting to quantify the impact of ERS contamination on potential future experiments with better angular resolution.

To do so we study an ideal case in which the instrument beam does not impact observations on the angular scales considered here, $32 < \ell\le600$. Specifically, we repeat our analysis assuming all frequency bands to be characterised by the same Gaussian beam with $\theta_{\mathrm{FWHM}} = 9'$. In this case we are able to mask much more of each PS without incurring a large penalty in terms of area loss, even when masking $\sim10^4$ sources. Therefore we use $f_{\mathrm{FWHM}} = 3$ (i.e. fully masking each source at all frequencies) for this analysis. Furthermore, the fact that, in this high-resolution setup, all the bands have the same beam, means that the astrophysical frequency scaling of the point sources is now preserved.

The results of this study are summarised in Table 3. We find that, unsurprisingly, the improved sensitivity on small scales increases all parameter biases due to residual contamination. In particular $A_{\mathrm{lens}}$ remains biased at the $\gtrsim3\sigma$ level even for $S_{\mathrm{lim}} = 0.1$ Jy, where all synchrotron parameters also exhibit large biases. However, a more stringent flux cut, at $S_{\mathrm{lim}} = 0.01$ Jy, suppresses these biases almost entirely, at no significant cost in the posterior parameter uncertainties. The only exception is the synchrotron tilt $\alpha_s$, which retains a $\sim2\sigma$ bias. $r$, as well as all dust parameters, remain unbiased in almost all cases, although the least stringent setup, with $S_{\lim}=1$ Jy recovers a $\sim2\sigma$ negative bias on $r$.

These results suggest that, although the use of a more stringent mask may be sufficient to reduce the most important sources of bias, we must inevitably marginalise over residual contamination to fully recover unbiased synchrotron constraints. We thus repeat our analysis including a fourth effective component, as described

in the previous section, characterised by a spectral index $\beta_p$, and a Poisson amplitude $A_p$. The results are summarised in Table 3. We find that including this additional component cancels the residual bias on $\alpha_s$ in the case with $S_{\lim}=0.01$ Jy (see also Figure 14), and suppresses almost all parameter biases for $S_{\lim}=0.1$ Jy and even 1 Jy, with the only exception of $A_{\mathrm{lens}}$ in the latter case, which retains a $\sim4.5\sigma$ bias. We did not find simple extensions to the effective point source models, such as including an effective frequency decorrelation parameter to characterise deviations from a perfect power law SED, to help reduce this bias. High-resolution experiments will therefore require a minimal level of polarised point source masking to achieve unbiased constraints.

Polarisation fraction scatter

The largest source of uncertainty in our model for the ERS contamination is the assumed polarisation fraction, which is not well constrained by current data. Throughout our results we have assumed a constant polarisation fraction across the sources of 2.8% as motivated in Lagache et al. (2020). We now explore a more conservative scenario in which the scatter in the distribution of polarisation fractions, as suggested by Puglisi et al. (2018), increases the effective fraction relevant to the residual ERS power spectrum to $\langle\Pi^2 \rangle^{1/2}=4.94\%$. This is the effective fraction, accounting for scatter, estimated for flat spectrum sources in Puglisi et al. (2018), which dominate the ERS emission at the frequencies studied here. This change corresponds to rescaling the ERS cross-spectra we used so far by a factor 3.1.

Unsurprisingly, this enhances the biases found in previous sections by a similar amount. For example, for $S_{\lim}=0.1$ Jy, the biases on $A_{\mathrm{lens}}$ and $\alpha_s$ grow from $\sim1\sigma$ to $\sim3\sigma$. The tensor-to-scalar ratio, as well as all dust parameters, remain unbiased. However, as before, all these biases disappear when including an effective Poisson-like component with a power-law SED to describe the residual ERS contamination. In this case, this effective component is detected at a higher (although still marginal) significance of $2.3\sigma$.

Dusty galaxies

We now consider the impact of polarised DSFGs, modelled as described in § 2.6. First, we compute the number counts for the DSFGs, see Figure 15. Contrary than what found for a fixed $S_{\mathrm{lim}}$ the number of sources decreases towards lower frequency, hence favouring masking at $\nu_m = 280$ GHz. And the number count differences among the different frequencies is more important for DSFGs than ERSs. Also, comparing with Figure 4, we find that there are more faint DSFGs than ERSs, but less bright DSFGs than ERSs. In particular, in the SO SAT sky area, at $S_{\mathrm{lim}} \simeq1$ Jy there are only $\sim6$ DSFGs at $\nu_m = 280$ GHz, while $\sim88$ DSFGs at $\nu_m = 27$ GHz. Such a low number of DSFGs above $S_{\mathrm{lim}} \simeq1$ Jy hints that we may avoid masking DSFGs at all.

Then, it is instructive to compare the ERS and DSFGs signals, both with each other and with Galactic foregrounds. Figure 16 shows, in the left panel, the residual ERS power spectrum for a flux limit $S_{\mathrm{lim}} = 0.1$ Jy together with the Galactic synchrotron signal in the SAT patch. In turn, the right panel shows the estimated polarised DSFG signal without masking any sources, together with the Galactic dust component. All signals are shown for the auto-correlation of the six frequency channels studied here. The figure also shows the lensing signal for $A_{\mathrm{lens}} = 0.3$, and the primordial $B$-mode signal for $r$ in the range $0.001 \le r \le0.01$ for comparison.

Already at the power spectrum level we can see that DSFGs are subdominant with respect to ERSs. At every frequency band, except 225 GHz and 280 GHz, the $BB$ power spectra of the DSFGs is lower than the ERS signal, even when we mask out ERSs with $S_{\mathrm{lim}} = 0.1$ Jy and do not mask any DSFGs. In addition to this, the Galactic dust signal is higher than Galactic synchrotron, and dominates over the DSFG contribution at almost all scales relevant for this analysis. Finally, the DSFG signal is lower than the $B$-mode lensing contribution at all relevant scales.

It is worth noting that the DSFG signal, as modelled here, reaches a minimum at around 68 GHz, growing towards higher frequencies with a dust-like modified black-body spectrum, but also towards lower frequencies, due to the synchrotron contribution from star-forming galaxies. Thus, the signal at 27 and 39 GHz is comparable to that in the higher-frequency channels. The larger beam limits the sensitivity to this contribution at the lower frequencies, however.

Finally, it is worth recalling that the DSFG signal discussed here contains both a shot-noise and clustered 1-halo contribution, although the former dominates over the latter here.

We propagate the DSFG residual shown in Figure 16 to parameter constraints as done before in the case of ERSs. We find that, even without masking any DSFGs, as in the spectra shown in the figure, their contribution is negligible for both the baseline and extended configurations, with no impact on the recovery of either the Galactic foreground or the cosmological parameters. No significant bias is found on any of the model parameters, including $\alpha_s$, $\alpha_d$ or $A_{\mathrm{lens}}$. Even in the high-resolution configurations explored in § 4.6, only a minor $0.3\sigma$ bias develops for $\alpha_s$, sourced by the radio component of the SED template used to describe DSFGs.

Conclusions

In this work we have studied the impact from polarised extragalactic point sources in ground-based $B$-mode analyses. Using the models of Tucci et al. (2011); Toffolatti et al. (1998); Lagache et al. (2020) for ERSs and Zagatti et al. (2024) for DSFGs, we have propagated the residual point source contamination through component separation into constraints on $r$, $A_{\mathrm{lens}}$, and Galactic foreground parameters. We have done so using a power-spectrum-based component separation pipeline, assuming a self-consistent masking strategy that uses the same point source mask at all frequencies.

We find that, while residual DSFG contamination does not lead to significant biases, ERS contamination must be mitigated and taken into account in near-future experiments, such as the Simons Observatory. Specifically focusing on small-aperture experiments, such as the SO SATs, we have shown that an effective strategy requires masking the bulk of the source contamination, without necessarily removing it entirely when the angular resolution varies substantially across frequency channels. Specifically, a flux cut of approximately $S_{\mathrm{lim}} = 0.1$ Jy at 27 GHz provides a suitable compromise between minimising bias and variance in the model parameters (although the exact value of $S_{\mathrm{lim}}$ may be optimised for individual instrumental араион configurations). More concretely, we find that, while the tensor-to-scalar ratio $r$ is largely unaffected by point source contamination, the lensing amplitude $A_{\mathrm{lens}}$ and the spectral tilt for Galactic synchrotron $\alpha_s$ may be significantly biased. We also show that any remaining parameter bias may be effectively suppressed by modelling the residual point source contamination as an additional component, with a power-law SED and a flat, Poisson-like, angular power spectrum, after marginalising over the free parameters of this model. Furthermore this is achieved at no cost in the statistical uncertainties on $r$, although the errors on $A_{\mathrm{lens}}$ and $\alpha_s$ approximately double.

We additionally investigate an hypothetical experiment with higher angular resolution, where point source masking results in a much milder loss of survey area. In this case, a more stringent flux cut, $S_{\mathrm{lim}} = 0.01$ Jy, achievable by high-resolution telescopes such as the SO LAT, provides the most effective mitigation strategy. This must be accompanied by the modelling and marginalisation over the residual contamination as described above.

Our analysis has focused only on the impact of point source contamination at the power spectrum level, and we have not explicitly considered the practical implications of introducing a large number of source holes in the analysis mask, beyond the increase in power spectrum uncertainties due to the loss of survey area. In practice, utmost care must be taken in the presence of complex masks to minimise the level of $E$-to-$B$ leakage, which can severely degrade the precision of the $BB$ power spectrum. Such effects will, however, need to be carefully accounted for in applications to real data, where extensive point source masking may be unavoidable. A natural extension of the approach proposed here would therefore be to combine point source masking with inpainting techniques, in order to mitigate these effects. Finally, future, more sensitive experiments, or those covering larger sky areas (such as LiteBIRD) may be more affected by point source contamination than the configurations explored here, and may require the use of more strict masking criteria, or the implementation of more sophisticated model to describe the residual point source contributions than the relatively simple effective power law behaviour proposed here.

Acknowledgements

AR thanks Gilles Weymann-Despres for stimulating scientific discussions and acknowledges support from the Dennis Sciama Junior Research Fellowship at Wolfson College, Oxford. GZ is funded by the European Union (ERC, ReLICS, project number 101116027). GZ acknowledges the financial support from the COSMOS network (www.cosmosnet.it) through the ASI (Italian Space Agency) Grants 2016-24-H.0 and 2016-24-H.1-2018.

Note that this is not an official Simons Observatory collaboration publication.

Data availability

No new data were generated or analysed in support of this research. This paper has been typeset from a TeX/LaTeX file prepared by the author.