© The Authors 2026

1. Introduction

The cosmic-microwave-background (CMB) photons, which decoupled from the baryon-photon plasma approximately 400 000 years after the Big Bang, encode information about the primordial fluctuations that seeded the large-scale structure of the Universe. In the previous decade, the Planck experiment successfully mapped the temperature anisotropies of the CMB with high precision, reaching sensitivity close to the cosmic-variance limit over a wide range of multipoles (Planck Collaboration I 2020; Planck Collaboration V 2020). In the new era of precision cosmology, next-generation experiments such as the Simons Observatory (SO; Lee et al. 2019), CMB-S4 (Abazajian et al. 2019), and LiteBIRD (Hazumi et al. 2019) are targeting the CMB polarisation field with unprecedented sensitivity up to a resolution scale of a few arcminutes for the first two experiments and around 20 arcminutes for the latter.

Inflationary models predict the production of the primordial gravitational waves (GWs) in the early Universe that imprint a large-scale divergence-free pattern – so-called B-modes – in the CMB polarisation field (Starobinskii 1979; Guth & Pi 1982; Linde 1983; Fabbri & Pollock 1983; Abbott & Wise 1984; Polnarev 1985; Kamionkowski et al. 1997; Seljak & Zaldarriaga 1997). Detecting these tensor B-mode signals is one of the primary science goals of current and upcoming CMB experiments. The detection of such signal would provide indirect evidence of inflation. Because no detection has been made far, only an upper limit on the amplitude of the primordial GWs can be constrained. Conventionally, this amplitude is expressed in terms of the so-called tensor-to-scalar ratio, r, which is defined as a ratio of the tensor (primordial GWs) and scalar perturbation power spectrum amplitudes for a characteristic pivot scale. Combined measurements from the BICEP/Keck array, Planck, and WMAP experiments give the current upper bound (at 2σ level) on the tensor-to-scalar ratio as r < 0.036 at 0.05 Mpc⁻¹ pivot scale (BICEP/Keck XIII 2021).

The detection of tensor B-modes is challenged by secondary effects such as extragalactic and Galactic foreground contaminations and weak gravitational lensing of the CMB. The intervening matter distribution deflects the CMB photons along its path, remapping the CMB anisotropies. Thus, weakly lensed CMB anisotropies trace the integrated matter density distribution of the Universe along the line of sight to the last scattering surface. Standard deviation of the deflection angle caused by the lensing effect is ∼2 arcminutes when the lensing potential wells along the path are not correlated (Lewis & Challinor 2006). The line-of-sight integrated deflection field can be reconstructed using quadratic estimators (QEs) – a method that is described well in the literature exploits the lensing-induced correlations present in observed CMB maps (Hu & Okamoto 2002; Okamoto & Hu 2003; Planck Collaboration XVIII 2014; Carron et al. 2017, 2022).

The lensing effect can introduce an additional B-mode polarisation pattern by re-mapping and distorting the primordial rotation-free polarisation patterns – so-called E-modes (Zaldarriaga & Seljak 1998; Kesden et al. 2002; Sayre et al. 2020). For the current upper bound on the tensor-to-scalar ratio and high signal-to-noise observations of polarisation fields, subtracting the lensing effect from B-mode polarisation becomes essential for the detection of the tensor B modes (Seljak & Hirata 2004; Smith et al. 2012; Simard et al. 2015; Sherwin & Schmittfull 2015; Baleato Lizancos et al. 2022). The subtraction of the lensing effect from observations, often called de-lensing, can be performed using a template of lensing B modes, which is constructed from an observed E mode and a high signal-to-noise estimate of the deflection field on small and intermediate scales, respectively (Manzotti et al. 2017; Namikawa 2017; BICEP/Keck XII 2021; Planck Collaboration VIII 2020).

With increasing signal-to-noise ratio of current and forthcoming CMB polarisation experiments, the polarisation-based lensing of potential power spectrum estimators, which uses the higher order statistics of the CMB polarisation field, has become more important (Namikawa et al. 2022; Belkner et al. 2024). However, these estimators are sensitive to secondary sources of polarised emission such as radio galaxies and dusty star-forming galaxies, which can bias the lensing reconstruction and therefore the delensed B-mode signal (Ferraro & Hill 2018; Schaan & Ferraro 2019; Sailer et al. 2021, 2023). A diffused Galactic foreground, and particularly polarised emission from thermal dust and synchrotron emission, pose an additional hindrance to delensing of CMB B-modes. These foregrounds display spatial variations and strong frequency dependence, and they are inherently non-Gaussian on both large and small scales (Gold et al. 2011; Fauvet et al. 2012; Planck Collaboration X 2016; Remazeilles et al. 2015; Planck Collaboration int. XXII 2015; Planck Collaboration IV 2020; Planck Collaboration XI 2020). The foreground contamination can introduce spurious correlations, which will bias both the lensing reconstruction and the estimation of the tensor-to-scalar ratio. Multi-frequency component-separation techniques can mitigate some of these effects, although residual foreground contamination can still propagate into delensed CMB maps (Fantaye et al. 2012; Beck et al. 2020; Belkner et al. 2024; Hertig et al. 2024; Abril-Cabezas et al. 2025; Bianchini et al. 2025).

In this work, we investigated the impact of Galactic foreground on the reconstruction of CMB-lensing power spectra, the de-lensing of B-mode maps, and the resulting constraints on the tensor-to-scalar ratio in detail. Particular attention was paid to the effects of non-Gaussian small-scale Galactic emissions. We quantified the biases introduced in the lensing reconstruction noise and the estimated lensing power spectra for three foreground models with varying degrees of small-scale non-Gaussianity. In this work, we did not include extragalactic foregrounds, as their impact on CMB-lensing reconstruction and de-lensing has already been widely studied in the literature, and in some cases mitigation methods have already been developed.

Using realistic simulations for a next-generation CMB-S4–like experiment, we explored how increasing foreground complexity affects the constraints on primordial gravity amplitudes. The paper is organized as follows. Section 2 presents the theoretical framework for component separation, lensing reconstruction, de-lensing, and the overall analysis methodology. The simulations employed are described in detail in Sect. 3. Section 4 discusses the impact of residual Galactic foreground on de-lensing efficiency and the resulting constraints on the tensor-to-scalar ratio. In Sect. 5, we conclude our findings, discuss the limitations of our work, and provide future prospects of an extension of our work.

2. Theoretical background and methodology

The CMB photons are deflected by the intervening matter distribution of the Universe along its journey from the last scattering surface to us. Weak lensing breaks the statistical isotropy of primordial CMB signal introducing coupling between different spherical harmonic modes. As a consequence of lensing, the CMB temperature and polarisation fluctuations observed in the direction $\widehat{\mathit{n}}$ in the sky are coming from primordial fluctuation fields in the $\widehat{\mathit{n}}+\mathrm{\nabla }\varphi (\widehat{\mathit{n}})$ direction, where $\mathrm{\nabla }\varphi (\widehat{\mathit{n}})$ is a deflection angle. Thus, the lensing causes the re-mapping of the primordial fluctuation fields as given by (Lewis & Challinor 2006)

$$\begin{array}{cc}\hfill T(\widehat{\mathit{n}})& =\stackrel{\sim }{T}(\widehat{\mathit{n}}+\mathrm{\nabla }\varphi (\widehat{\mathit{n}}))\hfill \\ \hfill [Q+iU](\widehat{\mathit{n}})& =[\stackrel{\sim }{Q}+i\stackrel{\sim }{U}](\widehat{\mathit{n}}+\mathrm{\nabla }\varphi (\widehat{\mathit{n}})).\hfill \end{array}$$(1)

Here, T denotes the lensed CMB temperature anisotropies and Q, U are the Stokes parameters of lensed CMB polarisation fields. Meanwhile, $\stackrel{\sim }{T},\stackrel{\sim }{Q},\stackrel{\sim }{U}$ represents the unlensed CMB fields. The lensing-potential field, $\varphi (\widehat{\mathit{n}}),$ is defined by (Lewis et al. 2000; Hu & Okamoto 2002; Okamoto & Hu 2003)

$$\begin{array}{c}\hfill \varphi (\widehat{\mathit{n}})=-2{\displaystyle {\int }_0^{{\chi }^{\ast }}}d\chi \frac{D_{\mathrm{A}}({\chi }^{\ast }-\chi )}{D_{\mathrm{A}}(\chi )D_{\mathrm{A}}({\chi }^{\ast })}\mathrm{\Psi }(\chi (\widehat{\mathit{n}}),\eta ),\end{array}$$(2)

which is the line-of-sight projection of the Weyl gravitational potential, $\mathrm{\Psi }(\chi (\widehat{n}))$, integrated along the comoving distance χ, from the observer to the last scattering surface at the comoving distance χ^* (≈14 000 Mpc). The quantity η denotes the conformal time when the photon was at distance of χ in the direction $\widehat{\mathbf{n}}$. D_A(χ) is the angular diameter distance that depends on the curvature of the Universe. For a flat Universe, the angular diameter distance is D_A(χ) = χ.

The unlensed temperature and polarisation fields and the lensing-potential field are assumed to be Gaussian and statistically isotropic. Hence, the statistical properties can be characterized by using the angular power spectra defined as (Zaldarriaga & Seljak 1997)

$$\begin{array}{cc}\hfill \langle {\stackrel{\sim }{X}}_{\ell m}^{\ast }{\stackrel{\sim }{Y}}_{{\ell }^{\prime }{m}^{\prime }}\rangle & ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}{\stackrel{\sim }{C}}_{\ell }^{\mathit{XY}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \langle {\varphi }_{\ell m}^{\ast }{\varphi }_{{\ell }^{\prime }{m}^{\prime }}\rangle & ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}{C}_{\ell }^{\varphi \varphi },\hfill \end{array}$$(4)

where ${\stackrel{\sim }{X}}_{\ell }^m$ and ϕ_ℓ^m are spherical harmonic coefficients of unlensed $\stackrel{\sim }{X}\in \{\stackrel{\sim }{T},\stackrel{\sim }{E},\stackrel{\sim }{B}\}$ fields and the lensing-potential field ϕ, respectively. The brackets ⟨.⟩ denote the ensemble average over multiple realisations, and δ_ij is the Kronecker delta symbol. The observed angular power spectra of T, E, and B also contain instrumental noise characterised by power spectra deconvolved for a beam:

$$\begin{array}{cc}\hfill N_{\ell }^{\mathit{TT}}& ={\sigma }_{\mathrm{T}}^{2}exp\frac{\ell (\ell +1){\theta }_{\text{fwhm}}^2}{8log2}=\frac{{\sigma }_{\mathrm{T}}^2}{B_{\ell }^2},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill N_{\ell }^{\mathit{XX}}& ={\sigma }_{\mathrm{P}}^{2}exp\frac{\ell (\ell +1){\theta }_{\text{fwhm}}^2}{8log2}=\frac{{\sigma }_{\mathrm{P}}^2}{B_{\ell }^2},\hfill \end{array}$$(6)

where XX ∈ {EE, BB }, B_ℓ is the beam-transfer function for a Gaussian profile, θ_fwhm is the angular resolution of the instrument expressed as the full width at half maximum (FWHM) of the beam, and σ_T and σ_P are the sensitivities corresponding to the intensity and polarisation detectors, respectively. The specification of noise spectra varies for different wavelength observations.

2.1. Component separation

Diffuse Galactic emission is one of the main sources of contamination in precise CMB measurements. The polarised Galactic foreground emission consists of two components that includes thermal dust emission and synchrotron radiation. Due to the alignment of dust grains with respect to the Galactic magnetic field (GMF; Stein 1966; Hildebrand 1988), the emission is partially linearly polarised (Planck Collaboration XI 2014; Planck Collaboration int. XXII 2015; Planck Collaboration XI 2020). The synchrotron emission is generated by relativistic electrons in the interstellar medium that spiral in the GMF. It is strongly polarised and dominates at low microwave frequencies (Haslam et al. 1982; Bennett et al. 2003; Gold et al. 2011; Planck Collaboration XXV 2016; Planck Collaboration IV 2020; Martire et al. 2022).

Multi-frequency observations and accurate component separation methods are essential for cleaning CMB signal from the Galactic foreground. A widely used method known as internal linear combination (ILC) recovers the CMB signal as a linear combination of multi-frequency maps with weights that minimise contributions from components of the maps other than the CMB (Bennett et al. 2003). The sum of the weights is constrained to be equal to one, which ensures that the clean CMB output map has unit response. This method of component separation is blind in the sense that it does not assume any prior knowledge about the foreground emissions. Bennett et al. (2003) showed the implementation of this method for the Wilkinson Microwave Anisotropy Probe (WMAP) temperature maps. Eriksen et al. (2004) extended the analysis to polarisation maps through the use of Lagrange multipliers. Similarly to pixel-based ILC, the harmonic ILC (HILC) method uses the spherical-harmonic coefficients of CMB anisotropies in a linear combination of the coefficients of different frequency maps (Tegmark et al. 2003; Kim et al. 2009). An advantage of the HILC method is that the Harmonic ILC weights depend on angular scale, whereas the pixel-based ILC method provides scale-independent ILC weights.

2.2. Quadratic estimator

Following the seminal work by Hu & Okamoto (2002) and its generalisation for the full sky analysis by Okamoto & Hu (2003), we implemented the quadratic estimator (QE) to reconstruct the lensing-potential field. The lensing of CMB photons breaks the statistical isotropy of CMB fluctuations introducing coupling between spherical harmonic modes. This gives rise to non-zero off-diagonal covariance between two lensed CMB fields, which up to first order in the lensing potential, ϕ, is given by

$$\begin{array}{c}\hfill {\langle X_{\ell m}{Y}_{{\ell }^{\prime }m\prime }\rangle }_{\mathrm{CMB}}={\displaystyle \sum _{\mathit{LM}}}{(-1)}^M\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& L\\ m& m\prime & -M\end{array}\right)f_{\ell L\ell \prime }^{\mathit{XY}}{\varphi }_{\mathit{LM}}.\end{array}$$(7)

The ensemble average 〈.〉_CMB is given over the lensed CMB fields assuming, contrary to the average introduced in Eq. (4), a fixed lensing-potential field. The lensing response functions, f_{ℓℓ′ L}^XY, for different pairs of fields are shown in Table 1. The coefficients, F_ℓLℓ′^(s), are defined as the mode-coupling ones due to the lensing (Okamoto & Hu 2003). In Eq. (7), we use the Wigner-3j symbol representation (Okamoto & Hu 2003; Smith et al. 2012).

An estimate of the lensing-potential field can be obtained by taking a weighted sum of quadratic combination of filtered observed fields as (Hu & Okamoto 2002; Okamoto & Hu 2003)

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{\mathit{LM}}^{\mathit{XY}}=\frac{A_L^{\mathit{XY}}}{L(L+1)}{\displaystyle \sum _{\ell m}}{\displaystyle \sum _{{\ell }^{\prime }{m}^{\prime }}}{(-1)}^M\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& L\\ m& m^{\prime }& -M\end{array}\right){\mathcal{W}}_{\ell L\ell \prime }^{\mathit{XY}}{\overline{X}}_{\ell m}{\overline{Y}}_{{\ell }^{\prime }{m}^{\prime }.}\end{array}$$(8)

Here, ${\overline{X}}_{\ell m}$ and ${\overline{Y}}_{\ell \prime m\prime }$ are inverse variance-filtered observed CMB fluctuation fields. The inverse-variance filtering of an observed field, ${\widehat{X}}_{\ell m}$, in harmonic space is given by

$$\begin{array}{c}\hfill {\overline{X}}_{\ell m}={(C_{\ell }^{\mathit{XX}}+N_{\ell }^{\mathit{XX}})}^{-1}{\widehat{X}}_{\ell m},\end{array}$$(9)

under the assumption that the noise covariance matrix is diagonal. In presence of foreground residual, the inverse weight becomes (C_ℓ^XX + N_ℓ^XX + F_ℓ^res)⁻¹, including the residual foreground power spectrum, F_ℓ^res.

The weights, 𝒲_ℓLℓ′^XY, in Eq. (8) are obtained by minimizing the variance $\langle |{\widehat{\varphi }}_{\mathit{LM}}{|}^2\rangle$ assuming a Gaussian distribution of CMB anisotropy. The normalisation factor is defined so that the expectation of the estimator is unbiased:

$$\begin{array}{c}\hfill \langle {\widehat{\varphi }}_{\mathit{LM}}\rangle ={\varphi }_{\mathit{LM}}.\end{array}$$(10)

The normalisation factor, $A_L^{\mathit{XY}}$, and the minimum variance weights, 𝒲_{ℓℓ′ L}^XY, are well described in Okamoto & Hu (2003). The angular power spectrum of the reconstructed lensing field is given by

$$\begin{array}{c}\hfill {\widehat{C}}_L^{\varphi \varphi ,XYX\prime Y\prime }=\frac{1}{2L+1}{\displaystyle \sum _M}{\widehat{\varphi }}_{\mathit{LM}}^{XY\ast }{\widehat{\varphi }}_{\mathit{LM}}^{X\prime Y\prime },\end{array}$$(11)

which, after averaging over ensemble, can be expressed as

$$\begin{array}{c}\hfill \langle {\widehat{C}}_L^{\varphi \varphi ,XY{X}^{\prime }{Y}^{\prime }}\rangle =C_L^{\varphi \varphi }+N_L^{XY{X}^{\prime }{Y}^{\prime }}.\end{array}$$(12)

Here, N_L^{XXX′ Y′} denotes the reconstruction noise power spectra associated with the quadratic estimator. The reconstruction noise depends on the instrumental noise levels and the resolution of the experiment. The noise term contains different bias contributions with different orders of dependency on $C_L^{\varphi \varphi }$:

$$\begin{array}{c}\hfill N_L^{XYX\prime Y\prime }=N_L^{(0)}+N_L^{(1)}+N_L^{(2)}+\cdots .\end{array}$$(13)

In the QE approach, the lensing-power spectrum is essentially estimated from the four-point function of the CMB fields. Following Wick’s theorem, the contribution from disconnected four-point contractions gives rise to the zeroth-order noise term, $N_L^{(0)}$, which have no dependency on $C_L^{\varphi \varphi }$. The $N_L^{(0)}$ noise would be present even in the case of purely Gaussian fields without any lensing effect present. The first-order $N_L^{(1)}$ noise is from connected lensed trispectrum contributions due to secondary lensing effects (Kesden et al. 2003). The $N_L^{(1)}$ bias has first-order dependency on $C_L^{\varphi \varphi }$ and becomes significant at small angular scales (L ≳ 1000). Following Hanson et al. (2011), we used lensed theoretical total power spectra in the calculation of $A_L^{\mathit{XY}}$ and 𝒲_ℓLℓ′^XY in Eq. (8) to mitigate higher-order bias contributions in ${\widehat{C}}_L^{\varphi \varphi }$.

The disconnected bias contribution term, $N_L^{(0)}$, is expressed analytically as

$$\begin{array}{cc}\hfill N_L^{(0)XYX\prime Y\prime }\left[C_{\ell }\right]& =\frac{A_L^{XY\ast }{A}_L^{X\prime Y\prime }}{2L+1}{\displaystyle \sum _{\ell {\ell }^{\prime }}}{\mathcal{W}}_{\ell L{\ell }^{\prime }}^{XY\ast }[{\mathcal{W}}_{\ell L{\ell }^{\prime }}^{X\prime Y\prime }{C}_{\ell }^{XX\prime }{C}_{\ell }^{YY\prime }\hfill \\ \hfill & +{(-1)}^{\ell +{\ell }^{\prime }+L}{\mathcal{W}}_{\ell L{\ell }^{\prime }}^{X\prime Y\prime }{C}_{\ell }^{XY\prime }{C}_{\ell }^{X\prime Y}],\hfill \end{array}$$(14)

where $C_{\ell }=\langle {\widehat{C}}_{\ell }\rangle$ is the expectation of estimated total power spectra with noise and foreground contributions.

To simplify analysis hereafter, we considered identical pairs of fields in both legs of the lensing-power-spectrum estimator; that is, we assumed that X = X′ and Y = Y′. Then, the $N_L^{(0)}$ noise follows:

$$\begin{array}{c}\hfill N_L^{(0)XY}[C_{\ell }]=A_L^{\mathit{XY}}.\end{array}$$(15)

In Fig. 1, we show the $N_L^{(0)}$ bias contribution for different QEs for CMB-S4-like experiments. Following Okamoto & Hu (2003), we also show a minimum-variance (MV) estimator combined from all different estimators that has the highest signal-to-noise ratio. Nevertheless, in our analysis we use only polarisation based EB estimator because it is the least sensitive to the mean-field effect arising from survey inhomogeneities.

Fig. 1. Analytical zeroth-order lensing reconstruction bias, $N_L^{(0)}$, for different estimators given by Eq. (14). The instrumental specification was considered for a CMB-S4-like experiment with a beam resolution of 2.5 arcminutes and white-noise temperature level of ΔT = 2 μK arcmin and polarisation of ΔP = 2.8 μK arcmin.

2.3. Mean-field effect

The reconstructed lensing-potential field using QE is sensitive not only to the anisotropy generated by the lensing effect, but also to other sources of statistical anisotropy. One of the most significant is anisotropy related to the incompleteness of full sky coverage for the ground-based observations and the removal of the parts of the sky most contaminated by Galactic foreground (Planck Collaboration XVII 2014; Hanson & Lewis 2009; Benoit-Lévy et al. 2013). For the CMB-S4-like survey, only the southern half of the sky will be observed. The sky coverage is shown in Fig. 2. Besides anisotropy from masking part of the sky, there are also anisotropies induced by inhomogeneous instrumental noise (Hanson et al. 2009) and asymmetric beams (Hanson et al. 2010). The latter is the least significant, so we neglected it in our analysis. The remaining anisotropies introduce systematic bias in the reconstructed lensing field, that is, the so-called mean field. We estimated the effect of masking and inhomogeneous noise in the mean-field signal, ${\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{MF}}$, by taking an average over different realisations of the lensing field reconstructed using the QE, $\langle {\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{QE}}\rangle$. The ensemble average of a random Gaussian lensing field, ${\widehat{\varphi }}_{\mathit{LM}}$, should go to zero, but the presence of anisotropies in the CMB fields used in the estimator leads to non-zero mean-field residuals. To obtain unbiased estimates of the lensing field, we subtracted the mean-field from the quadratic estimate of the lensing potential (Planck Collaboration XVIII 2014):

$$\begin{array}{cc}\hfill {\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{MF}}& =\langle {\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{QE}}\rangle ,\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill {\widehat{\varphi }}_{\mathit{LM}}& ={\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{QE}}-{\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{MF}}.\hfill \end{array}$$(17)

Fig. 2. Masks used in our analysis to mimic CMB-S4 sky coverage, shown in Ecliptic coordinates with Galactic emission at 145 GHz in the background. (a) The left panel shows the mask used for component separation for the wide sky survey corresponding to the LAT configuration (fsky = 0.64). (b) The middle panel shows the mask used for lensing reconstruction for the wide sky survey (fsky = 0.37). (c) The right panel shows the sky area considered for the deep survey corresponding to the SAT configuration, which was used for the delensing study (fsky = 0.025).

2.4. Foreground effects

The presence of Galactic foreground affects lensing reconstruction in two ways. Firstly, it modifies the reconstruction noise power spectrum associated with the QE due to additional power of the foreground signal. In particular, assuming Gaussianity of the foreground, the reconstruction noise defined in Eq. (14) becomes a function of the total power spectra that includes the foreground angular power spectrum, F_ℓ, given as (Beck et al. 2020)

$$\begin{array}{c}\hfill N_L^{(0),XY{X}^{\prime }{Y}^{\prime }}\left[{\widehat{C}}_{\ell }\right]\to N_L^{(0),XY{X}^{\prime }{Y}^{\prime }}[{\widehat{C}}_{\ell }+F_{\ell }].\end{array}$$(18)

On the other hand, non-isotropic distribution of the foreground emission can contribute to the lensing-potential estimation biasing the estimator. The power spectrum of this term, $F_L^{\mathrm{syst}}$, is also a systematic bias for the lensing-potential power spectrum (Beck et al. 2020). As it depends on the connected four-point correlation function of the foreground emission, the bias is very sensitive to the level of non-Gaussianity of the foreground. The $F_L^{\mathrm{syst}}$ term can be computed by implementing the quadratic estimator on the Galactic-foreground emission maps. We also considered systematic bias arising from residual foreground present in the HILC cleaned CMB maps and denoted it as $F_L^{\mathrm{syst},\mathrm{res}}$. Similarly to the $F_L^{\mathrm{syst}}$ term, $F_L^{\mathrm{syst},\mathrm{res}}$ was estimated from the residual foreground maps.

2.5. B-mode template de-lensing

In the Taylor expansion of the lensed polarisation fields, the lensed B modes at first order in the lensing potential field become

$$\begin{array}{c}\hfill B_{\ell m}^{\mathrm{len}}={\displaystyle \sum _{l_1{m}_{1}LM}}{g}_{\ell l_1{l}_2}^{\mathit{EB}}\left(\begin{array}{ccc}\ell & l_1& L\\ m& m_1& M\end{array}\right)E_{l_1{m}_1}^{\ast }{\varphi }_{\mathit{LM}}^{\ast },\end{array}$$(19)

which corresponds to convolution between unlensed E modes, E_ℓm, and lensing potential, ϕ. The coupling coefficients, g_ℓl1L^EB, are defined as (Smith et al. 2012)

$$\begin{array}{c}\hfill g_{\ell l_{1}L}^{\mathit{EB}}=\frac{F_{\ell l_{1}L}^{(-2)}-F_{\ell l_{1}L}^{(2)}}{2i}.\end{array}$$(20)

The mode coupling due to the lensing is captured by F_ℓℓ′L^( ± 2) (Smith et al. 2012; Okamoto & Hu 2003). The first-order approximation is sufficient for de-lensing studies at large angular scales relevant for the constraints on the tensor-to-scalar ratio (Smith et al. 2012; Challinor & Lewis 2005).

To de-lens observed B modes, a template of lensed B modes can be obtained using the Wiener-filtered observed E modes, $E_{\ell m}^{\mathrm{obs}}$, and the Wiener-filtered reconstructed lensing potential, $\widehat{\varphi }$. Thus, the lensed B-mode template at gradient order is given by (Smith et al. 2012)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill B_{\ell m}^{\mathrm{temp}}& ={\displaystyle \sum _{\ell \prime m\prime LM}}{g}_{\ell \ell \prime L}^{\mathit{EB}}\left(\begin{array}{ccc}\ell & \ell \prime & L\\ m& m\prime & M\end{array}\right)\left(\frac{C_{{\ell }^{\prime }}^{\mathit{EE}}}{C_{{\ell }^{\prime }}^{\mathit{EE}}+N_{{\ell }^{\prime }}^{\mathit{EE}}}\right)E_{{\ell }^{\prime }{m}^{\prime }}^{obs\ast }\hfill \\ \hfill & \times \left(\frac{C_L^{\varphi \varphi }}{C_L^{\varphi \varphi }+N_L^{\varphi \varphi }}\right){\widehat{\varphi }}_{\mathit{LM}}^{\ast }.\hfill \end{array}\end{array}$$(21)

In the presence of foreground residuals, we used residual foreground and noise power spectra in the Wiener filtering of E modes to suppress the contribution of highly contaminated modes. Then, the de-lensed B-mode map is given by,

$$\begin{array}{c}\hfill B_{\ell m}^{\mathrm{del}}=B_{\ell m}^{\mathrm{obs}}-B_{\ell m}^{\mathrm{temp}}.\end{array}$$(22)

Baleato Lizancos et al. (2021b) showed that using lensed E modes in Eq. (21) leads to a cancellation of higher order terms in the de-lensed power spectra, thereby reducing the residual lensing floor. Other widely used de-lensing methods that construct B-mode templates from non-perturbative remapping of observed E-modes have been shown to leave one-order-higher residual lensing compared to gradient-order templates (Baleato Lizancos et al. 2021b). Anti-lensing approaches, in which the observed E modes are first de-lensed by inverse lensing methods and then used to generate a non-perturbatively re-mapped lensed B-mode template (Han et al. 2021), achieve similar de-lensing efficiency to gradient-order templates, but they are computationally more expensive. The gradient-order template method adopted in this paper is therefore well suited for forecasting studies due to its analytical simplicity and computational efficiency.

The angular power spectra for de-lensed B-modes can be written as

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{del}}=C_{\ell }^{\mathrm{tens}}+C_{\ell }^{\mathrm{res}}+N_{\ell }^{\mathrm{del}}+N_{\ell }^{\mathrm{res}}+F_{\ell }^{\mathrm{res}},\end{array}$$(23)

where C_ℓ^tens is the primordial tensor component power spectrum, C_ℓ^res is the lensing residual spectra, N_ℓ^del is the de-lensing bias, N_ℓ^res is the residual instrumental noise spectra, and F_ℓ^res is the residual foreground spectra after component separation. The de-lensing bias includes contributions from six-point and higher order correlation functions of CMB fields present in the convolution of observed E-mode and internally reconstructed $\widehat{\varphi }$ fields using pairs of E- and B-mode maps (Baleato Lizancos et al. 2021a).

For our analysis, we estimated the residual lensing spectra, C_ℓ^res, from Monte Carlo (MC) simulations of CMB maps by taking the difference between lensing only the B-mode map and lensing B-mode templates over a set of 400 simulations. A power spectrum of the difference also contains the de-lensing-bias contribution, N_ℓ^del, by definition, so this term in Eq. (23) can be omitted. The residual instrumental noise and foreground power spectra were computed by taking an average over a set of 100 simulations for each of the three foreground models described in Sect. 3.1. All the power spectra used in our analysis were computed on a cut sky corresponding to the SAT configuration (sky fraction, f_sky = 0.025) using the Namaster¹ package, which implements the MASTER algorithm (Hivon et al. 2002; Alonso 2025).

2.6. Likelihood analysis

We used the power spectrum of cleaned and de-lensed B-mode map to constrain the tensor-to-scalar ratio, r. In our analysis, we adopted the Hamimeche–Lewis (H–L) likelihood approximation, which closely reproduces the exact full-sky likelihood (Hamimeche & Lewis 2008). Unlike a naive Gaussian approximation, the H-L likelihood captures the non-Gaussian statistics introduced by cut-sky analyses and delivers improved accuracy. For a single B-mode polarisation field, the log-likelihood function is written as

$$\begin{array}{c}\hfill -2ln\mathcal{L}\left(r|{\widehat{C}}_{\ell }\right)={\displaystyle \sum _{l{l}^{\prime }}}\left[g({\widehat{C}}_{\ell }/C_{\ell })C_{\mathit{fl}}\right]{\left[{\mathbf{M}}_f^{-1}\right]}_{l{l}^{\prime }}\left[g({\widehat{C}}_{{\ell }^{\prime }}/C_{{\ell }^{\prime }})C_{f{\ell }^{\prime }}\right],\end{array}$$(24)

where

$$\begin{array}{c}\hfill g(x)=\mathrm{sign}(x-1)\sqrt{2(x-ln(x)-1)}.\end{array}$$

[M_f] is the covariance matrix estimated from a set of 400 simulated CMB realisations that include a purely Gaussian Galactic foreground. These Gaussian-foreground simulation sets are described in details in Sect. 3.1. The covariance matrix computed from a limited number of simulations can be biased (Hartlap et al. 2007). In our analysis, we used the Hartlap factor to obtain an unbiased estimation of the inverse covariance matrix given as

$$\begin{array}{c}\hfill {\mathbf{M}}_f^{-1}=\frac{N-d-2}{N-1}{\mathbf{M}}_{\ast }^{-1},\end{array}$$(25)

where N is the number of simulations and d is the dimension of the covariance matrix estimator. Here, M_* denotes the biased covariance matrix estimator given as

$$\begin{array}{c}\hfill {[{\mathbf{M}}_{\ast }]}_{\ell {\ell }^{\prime }}=\langle ({\widehat{C}}_{\ell }-C_{\mathit{fl}})({\widehat{C}}_{{\ell }^{\prime }}-C_{f{\ell }^{\prime }})\rangle .\end{array}$$(26)

The model angular power spectrum is given by

$$\begin{array}{c}\hfill C_{\ell }=r{C}_{\ell }^{\mathrm{tens},r=1}+C_{\ell }^{\mathrm{res}}+N_{\ell }^{\mathrm{res}}+F_{\ell }^{\mathrm{res}},\end{array}$$(27)

where C_ℓ^{tens, r = 1} is the theoretical power spectrum of primordial GWs with r = 1 and spectral index n_t = 0. We define a fiducial model C_fl = C_ℓ by setting r = 3 × 10⁻³, which is the same as our input value of r for the simulated observations. We considered the covariance matrix to be effectively model independent, as it exhibits no significant variation for changes in r around the fiducial value, and the primordial tensor component contributes negligibly to the uncertainty in the likelihood analysis compared to other components, as shown in Fig. 9.

3. Simulations

To simulate realistic sky observations, we considered lensed CMB maps, instrumental noise, and foreground emission. We calculated theoretical unlensed CMB spectra and lensing potential power spectra with the publicly available CAMB² package (Lewis et al. 2000). We used a standard Λ-cold dark matter (Λ-CDM) cosmology with the cosmological parameters from the Planck 2018 release (Planck Collaboration VI 2020) : Ω_ch² = 0.120, Ω_bh² = 0.022, n_s = 0.965, τ = 0.054, H₀ = 67.4. We considered sets of simulations for two tensor-to-scalar ratios: r = 0 and r = 3 × 10⁻³. Using the healpy package³ (Górski et al. 2005), we generated full sky maps of the Stokes parameters Q and U from the obtained theoretical CMB power spectra and Gaussian random realisations of lensing-potential field from the lensing power spectrum. The CMB maps were then lensed using publicly available lenspyx⁴ package that provides a curved-sky lensing algorithm for the full sky (Reinecke et al. 2023). We simulated all maps at a resolution parameter of N_side = 2048. The lensed full-sky Stokes parameter maps were converted to full-sky E- and B-mode maps and were then convolved with Gaussian beams corresponding to the seven frequency bands of a CMB-S4–like reference experiment, as listed in Table 2. Subsequently, we added beam-smoothed, simulated E- and B-mode thermal dust and synchrotron polarisation maps to the respective frequency channels. In the Monte Carlo simulations, we kept the Galactic-foreground emission templates fixed and identical across all realisations. We provide the details of the foreground emission maps and their frequency scaling in subsequent section. We then added Gaussian E- and B-mode noise maps and then applied the appropriate survey masks for different telescope configurations.

3.1. Galactic foreground

Linearly polarised thermal dust emission dominates the CMB signal at frequencies higher than 90 GHz, where upcoming CMB surveys will operate, while the synchrotron emission dominates the polarisation signal at lower frequencies.

The intensity of the thermal dust emission follows the scaling law of a modified-black-body (MBB) emission (Planck Collaboration int. XXII 2015):

$$\begin{array}{cc}\hfill I_{\nu }^{\mathrm{dust}}& =I_{{\nu }_0}^{\mathrm{dust}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_d}\frac{B_{\nu }(T_d)}{B_{{\nu }_0}(T_d)}.\hfill \end{array}$$(28)

Here, I_ν0^dust is a template of dust emission at a reference frequency, ν₀. The spectral index, β_d, and dust temperature, T_d, have typical values of 1.5 and 20 K, respectively, for both intensity and polarisation, but a spatial variation of the order of 10% in high Galactic latitude is observed (Planck Collaboration XI 2020; Planck Collaboration X 2016). The Stokes parameter Q and U maps also follow the same MBB law with a frequency-dependent polarisation fraction that determines the level of polarisation degree. The polarisation fraction also shows spatial variations and is taken into account in the modelling on small scales using a polarisation-fraction tensor formalism (Pan-Experiment Galactic Science Group 2025).

Synchrotron emission arises from relativistic electrons accelerated along the Galactic magnetic field. Synchrotron emission follows a power law (Gold et al. 2011):

$$\begin{array}{c}\hfill I_{\nu }^{\mathrm{sync}}=I_{{\nu }_0}^{\mathrm{sync}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_s}.\end{array}$$(29)

Similar scaling is applicable to Q and U polarisation fields. Here, the power-law index, β_s, varies across the sky with values ranging from β_s = −3.0 near the Galactic plane to β_s = −3.3 off the plane.

We considered three different pairs of models of Galactic emission, each with different levels of complexity. We used the PySM3⁵ package to implement the scaling law for the foreground maps (Thorne et al. 2017; Zonca et al. 2021; Pan-Experiment Galactic Science Group 2025). The foreground models are listed and described in detail below.

3.1.1. d1s1:

This is the base model where the spectral indices used in the frequency scaling of dust and synchrotron emission have spatial variations. The reference dust templates used are the 545 GHz intensity maps and 353 GHz polarisation maps from the Planck 2015 results (Planck Collaboration X 2016). The dust template maps were smoothed to 2.6 degrees. The synchrotron emission templates for temperature and polarisation were taken from the Haslem maps at 408 MHz (Remazeilles et al. 2015) scaled to 23 GHz and the WMAP nine-year data at 23 GHz (Bennett et al. 2013). Both templates were smoothed to five-degree resolution. Small scales with Gaussian distributions were then added to both dust and synchrotron maps by fitting a power law to the template maps.

3.1.2. d10s5:

The thermal dust model d10 and synchrotron model s5 are of medium complexity, with amplitude-modulated small scales chosen to imitate non-Gaussian behaviour. The thermal dust-emission template was taken from the component-separated, thermal dust-emission maps using the Generalised Needlet ILC (GNILC) method from the Planck survey (Planck Collaboration IV 2020). These template maps were obtained at 353 GHz with a resolution of 21.8′ in terms of temperature and a variable resolution of 21.8′−80′ for polarisation maps. The s5 model is similar to s1, except for the re-scaling of the spectral index map based on the S-band polarisation All Sky Survey (S-PASS; Kamionkowski et al. 1997). Small scales were added to both emission models following the Logarithm of the polarisation-fraction tensor (logpoltens) formalism, which is well described in the PySM documentation (Pan-Experiment Galactic Science Group 2025). The small scales added have a modulation adopted according to the intensity maps, but we excluded the parts of the sky with the strongest Galactic emission to avoid excessive power in the high-intensity regions of the Galactic plane. The logpoltens formalism and the intensity modulation introduce some non-Gaussianity to the added small angular scales, which is important for our analysis. Additionally, the spectral indices of dust and synchrotron emission and dust temperature maps also contain random fluctuations added to small scales.

3.1.3. d12s7:

The d12 model is a three-dimensional line-of-sight-integrated model consisting of six layers, with each layer having a different realisation of dust emission at ν₀ = 353 GHz, different dust temperature maps, and spatially varying spectral indices for frequency scaling (Martínez-Solaeche et al. 2018). On large scales, the model matches observed dust emission, and on smaller scales it is extrapolated from polarised dust-emission power spectra. The small-scale realisations for each layer were modulated on a large-scale polarised intensity level to obtain non-Gaussian statistics. The synchrotron model s7 is a more complex version of the s5 model, with the introduction of a curvature term to the spectral index in Eq. (29):

$$\begin{array}{c}\hfill I_{\nu }^{\mathrm{sync}}=I_{{\nu }_0}^{\mathrm{sync}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_s+Cln(\nu /{\nu }_0)}.\end{array}$$(30)

The curvature term, C, was defined for the reference frequency ν₀ = 23 GHz, which is based on the measurement from the balloon-borne experiment of the Absolute Radiometer for Cosmology, Astrophysics, and Diffuse Emission (ARCADE-2; Fixsen et al. 2011; Kogut 2012). The dust model in d12s7 is of the highest complexity, as all individual layers contain intensity-modulated small scales, deviating from Gaussian statistics. The frequency dependency also contributes to non-Gaussianity due to different realisations of spatially varying spectral indices in each layer. The synchrotron emission in d12s7 contains an additional curvature dependence in frequency scaling, which increases its complexity compared to the simple power law in the d1s1 and d10s5 cases.

To estimate the covariance matrix used in the likelihood analysis (Eq. (26)), in the case of each considered foreground model we used 400 simulations of purely Gaussian foreground maps. We fitted a power law to the power spectra of thermal dust and synchrotron emission at each frequency, and then used these fitted spectra to produce random Gaussian realisations of foreground maps. These maps were added to the CMB realisations along with their corresponding noise realisations, which were then processed through the HILC and de-lensing pipelines to obtain the de-lensed power spectra. Finally, using Eq. (26), we computed the covariance matrices for each model, capturing the off-diagonal contributions arising from mask-induced anisotropies. Additionally, we generated one separate set of simulations with CMB-S4–like noise properties that are free of foreground contamination, but processed in a consistent way through the HILC pipeline; we used these to validate our pipeline.

3.2. Sky coverage and instrumental noise

In our simulated maps, we considered observation of the sky area corresponding to the CMB-S4 experiment. It covers around 60% of the sky and overlaps with footprints of large-scale structure surveys such as the Vera Rubin Observatory (VRO) Legacy Survey of Space and Time (LSST; LSST Science Collaboration 2020) and the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016). Significant overlap is especially important for future cross-correlation studies between the CMB lensing potential and the large-scale structure tracers. We considered the two sky cuts shown in Fig. 2 for two separate setups for telescopes used for observations: small-aperture telescopes (SATs) and large-aperture telescopes (LATs). The SATs are considered to perform ultra-deep observations targeting large-scale B-mode signals where primordial B modes are expected to have the highest power. Meanwhile, the LATs observe the small-scale polarisation field in a wide field, which is important to lensing science. The weak lensing information from LATs is relevant to de-lensing analyses and large-scale structure studies.

The HILC method was implemented on the cut-sky for both LAT and SAT configurations. We applied the corresponding masks to the scalar E- and B-mode fields directly to mitigate the E − B mixing problem arising from the cut-sky transformation of the Stokes parameters Q and U (Lewis et al. 2001; Bunn et al. 2003). In actual observations, one has to be careful regarding the E − B mixing problem when working with cut-sky data. For the LAT configuration, we used the mask shown in the left panel of Fig. 2, which retains 65% of the sky to exclude the highly contaminated Galactic-plane region. In the HILC implementation, we used cut-sky harmonic coefficients, which led to higher residual contamination around the mask edges due to the strong signal gradient at the boundary, which is commonly referred to as the ringing effect (Grain et al. 2009). We found that this spurious HILC residue near the mask edge introduces additional biases in the lensing reconstruction. To mitigate this, we employed an extended mask by removing a narrow band around the edges of the component-separation mask. We used this final mask in our lensing reconstruction, which covers 37% of the sky, as shown the middle panel of Fig. 2. For the SAT configuration, we considered a 3% sky-coverage region located in a low-Galactic-contamination area. To further suppress residual Galactic contamination, we adopted a narrower mask with 2.5% sky coverage excluding a band around the edges. This mask, shown in the right panel of Fig. 2, was used for the de-lensing and for the likelihood analysis.

The instrumental noise contributions are assumed to be Gaussian. The instrumental properties for different frequencies are listed in Table 2. The temperature and polarisation noise power spectra follow Eqs. (5) and (6), with an additional modification at low multipoles due to atmospheric noise given as

$$\begin{array}{c}\hfill N_{\ell }^X=\frac{{\sigma }_X^2}{B_{\ell }^2}(1+{\left(\frac{\ell }{{\ell }_{\mathit{knee}}}\right)}^{\alpha }),X\in \{T,E,B\}.\end{array}$$(31)

Here, B_ℓ is the beam transfer function, and σ_T and ${\sigma }_{\mathrm{E},\mathrm{B}}=\sqrt{2}{\sigma }_{\mathrm{T}}$ are the temperature and polarisation noise sensitivity, respectively. We took ℓ_knee and α parameter values from the CMB-S4 DRAFT⁶ tool. The parameter, ℓ_knee, and power-law index, α, for SATs and LATs in our polarisation measurements are listed in Table 3. In Fig. 3, we show the beam-de-convolved noise spectra for different frequencies corresponding to LAT specifications. We used the noise properties and sky coverage of the LATs for the lensing reconstruction, employing the extended LAT mask after additional edge cuts. For the de-lensing study, we used mock observations with the noise properties of the SAT configuration, while the lensing potential used for de-lensing was estimated from the LAT configuration.

Fig. 3. Beam-de-convolved noise power spectra for temperature (solid) and polarisation (dashed) in LAT configuration.

4. Results

4.1. Harmonic ILC products

The Harmonic ILC method outputs a single minimum-variance combination of observations from all frequency channels that have unit responses to the CMB signal by definition. To combine the maps from different frequency channels in spherical harmonic space, the maps were deconvolved to a common resolution, usually to the highest resolution available. However, we used the common resolution of 2.5 arcminutes to avoid any artefacts arising from de-convolving in harmonic space, especially for the 20, 30, and 270 GHz channels as their input resolutions show a huge difference compared to the highest resolution. As mentioned in Sect. 3.2, for the maps corresponding to the LAT configuration we used a 5% masking in the Galactic plane to avoid the most heavily contaminated region. We considered the HILC weights to be the same in the unmasked part of the sky. Dividing the sky into different regions to obtain the spatially varying HILC weights can further improve the component separation.

For the LAT configuration, the E modes are signal-dominated up to multipoles of L = 3000, while the lensing B modes are signal-dominated in the range of 200 < L < 1500. The residual foreground only dominates the B modes beyond L = 3000. In our lensing reconstruction analysis, we therefore used modes up to L = 3000. To suppress foreground contamination at high L, we fitted a simple noise model given by Eq. (31) to the residual foreground power spectra of the E and B modes, and we used the fitted theoretical model for inverse weighting of the input maps in Eq. (9). A detailed description of HILC products and the HILC weights for LAT configuration is provided in Appendix A.

The tensor-to-scalar ratio constraints were derived from the deep observations of the SATs. The SAT maps at each frequency were weighted by a normalised hit-count map. In Fig. 4, we compare the average level of foreground contamination in the 95 GHz and 145 GHz channels, which is approximately four orders of magnitude higher than the residual foreground level in the component-separated maps. The HILC-cleaned B-mode maps and residual foreground for all three foreground models are shown in Fig. A.4. The foreground residues are about one order of magnitude smaller than the signal.

Fig. 4. Angular power spectra for CMB B-mode-cleaned maps for the SAT configuration and sky coverage. The solid grey zigzag lines are BB power spectra of 100 simulations, and the average is shown as red stars. The theoretical input power spectra are shown as solid black lines. The residual HILC noise and foreground levels are shown as dotted and dashed lines, respectively. The solid magenta line shows the average de-lensed B-mode spectra over 100 simulations. The dashed blue line shows the average power spectra of Galactic-foreground B modes at the frequency channels 95 and 145 GHz. The dashed grey line shows the tensor B-mode level for r = 0.003.

4.2. Mean-field estimation

We implemented the EB estimator on the HILC-cleaned E-mode and B-mode maps for the LAT configuration and sky coverage. We computed the mean-field effect on lensing reconstruction using a set of 100 simulations as described in Sect. 3. We considered homogeneous instrumental noise and an azimuthally symmetric beam, so the dominant contribution to the mean-field effect arises from masking due to incomplete sky coverage. In addition, the Galactic foreground can be regarded as an anisotropic contamination to the CMB signal, which can further enhance the mean-field bias. To quantify the impact of these foregrounds, we compared two cases. In the first case, we used co-added, foreground-contaminated, noise-free CMB maps of 95 + 145 GHz before component separation. These co-added maps include the average foreground contribution of the 95 GHz and 145 GHz channels. We then added residual HILC noise realisations to these maps and performed lensing reconstruction. We refer to these foreground contaminated mock observations as the medium-frequency (MF) maps in this paper. The resulting power spectrum of the mean-field effect for MF maps before component separation is shown in red in Fig. 5 for the d10s5 model. In the second case, we evaluated the mean-field effect after component separation for both the d10s5 and d12s7 foreground models. The comparison demonstrates that component separation reduces the overall mean-field contribution by a factor of two. However, the mean-field level is comparable to the $N_L^{(1)}$ reconstruction noise at low multipoles, and it is non-negligible even after component separation. The mean-field contamination does not depend on the foreground models because power spectra for the two models overlap with each other. This is expected as the mean-field effect is mostly dominated by the anisotropies introduced due to masking. The power spectrum of the mean field in all cases is comparable to the $N_L^{(1)}$ bias at low ℓ ≲ 100, underlining the necessity of carefully modeling and subtracting it to achieve an unbiased lensing reconstruction. In Fig. 5, we show the mean-field effect for the EE quadratic estimator of component-separated E-mode maps corresponding to d10s5 model. We can notice that the effect dominates the lensing signal at low ℓ. One of the main reasons why we did not consider a minimum-variance (MV) combination of all quadratic estimators is discussed in Sect. 2.2. The mean-field effect contributions from same-pair estimators such as the TT and EE estimators makes the reconstruction with an MV estimator noisy at low multipoles, and a careful subtraction is needed in that case.

Fig. 5. Angular power spectra of reconstructed lensing field shown as a solid blue line. The mean field effect on the EB estimator for MF maps for d10s5 and HILC-cleaned maps for d10s5 and d12s7 are shown as solid red, green, and grey lines, respectively. The mean-field effect on the EE estimator for HILC-cleaned maps of d10s5 is shown as a solid cyan line. Theoretical lensing power spectra are shown as solid black lines. The dashed magenta line is the sum of theory power spectra and lensing reconstruction noise terms, N(0) and N(1), shown as dashed (orange) and dash-dotted (purple) lines, respectively.

4.3. Foreground bias

Contamination from Galactic foreground increases the reconstruction noise in the lensing power spectrum coming from Gaussian contribution and also introduces a systematic bias to the reconstructed lensing power spectrum, denoted as $F_L^{\mathrm{syst}}$, which arises from non-Gaussian small scales. We performed lensing reconstruction both before and after component separation to understand these biases in detail.

We defined the fractional change in the zeroth-order reconstruction noise as

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{N}_L^{(0)}}{N_L^{(0)}}=\frac{N_L^{(0),\mathrm{FG}}-N_L^{(0)}}{N_L^{(0)}},\end{array}$$(32)

where, $N_L^{(0),\mathrm{FG}}$ is a noise power spectrum with foreground (or residual foreground) contamination, and $N_L^{(0)}$ is a noise spectrum without the presence of foreground contamination. The $N_L^{(0)}$ was computed analytically using Eq. (14), and $N_L^{(0),FG}$ was computed following Eq. (18) using foreground (or residual foreground) power spectra. In Fig. 6, we compare the biases due to the foreground for MF maps before component separation and for the HILC products after component separation, which contains a residual foreground. We used the same noise realisations in both cases, before and after component separation, while generating simulated maps to only compare the impact of foreground.

Fig. 6. Relative bias in semi-analytic $N_L^{(0)}$ noise due to foreground contamination in MF channels (solid lines) before component separation and residual foreground in HILC maps (dashed lines) after component separation for different models.

The d1s1 model induces the largest increase in the reconstruction noise prior to component separation, with unity-order changes at low multipoles. For the more realistic d10s5 and d12s7 models with non-Gaussian small scales, the impact on reconstruction noise is less severe, but still non-negligible, leading to a 10–20% increase in the reconstruction noise across multipoles. In the case of the HILC-cleaned CMB maps, the bias due to residual foreground significantly drops to 2–3%. Although unmitigated foreground can substantially bias the reconstruction noise, component separation methods are effective in suppressing these contributions.

In Fig. 7, we show the systematic bias, $F_L^{\mathrm{syst}}$, for the MF and HILC-cleaned CMB maps. The results show that the systematic bias for the MF map before component separation is below the cosmic variance limit of the total reconstructed lensing power spectrum, including the lensing potential spectrum ($C_L^{\varphi \varphi }$) and reconstruction noise ($N_L^{(0)}+N_L^{(1)}$). The component separation further reduces the systematic bias to three orders of magnitude lower than $N_L^{(0)}$, making it negligible for our analysis. For highly non-Gaussian models such as the d12s7 model, we see the largest $F_L^{\mathrm{syst}}$ bias. A comparison of Figs. 6 and 7 shows that the $\mathrm{\Delta }{N}_L^{(0)}$ bias in reconstruction noise is more pronounced for foreground models dominated by the small-scale Gaussian structure, while models with stronger non-Gaussianity contribute significantly to the systematic bias in the reconstructed power spectra. In Fig. 7, we also present the bias coming from the Gaussian contribution of $\mathrm{\Delta }{N}_L^{(0)}$ for the d12s7 model after component separation, which is almost two orders of magnitude larger than corresponding systematic bias arising from non-Gaussian foreground residues, and it is comparable to the cosmic-variance limit.

Fig. 7. The $F_L^{\mathrm{syst}}$ bias is shown for MF maps (dark colours) and HILC-cleaned CMB maps (light colours) for three different models. The black curve shows the fiducial lensing power spectrum, and the dashed red curve is the sum of reconstruction noises, N(0) + N(1). The dash-dotted magenta line is the bias in reconstruction noise, $\mathrm{\Delta }{N}_L^{(0)}$, due to residual Galactic foreground after component separation for the d12s7 model. The shaded grey region represents the cosmic variance limit on the reconstructed lensing power spectrum.

4.4. De-lensing and r constraint

Implementing Eq. (21), we obtained lensing B-mode templates using Wiener-filtered, foreground-cleaned (HILC) E-mode maps and the reconstructed lensing potential field over the LAT configuration sky coverage (f_sky = 0.37). We included residual foreground and noise power spectra in the Wiener filtering. We used this template to de-lens the B-mode map using Eq. (22) for the deep-field observations from the SAT configuration. Subtracting the B-mode template does not fully remove the lensing-induced modes, leaving some residual lensing signal as shown in Fig. C.1. The de-lensed B-mode maps also include noise and residual foreground components, which contribute to the measured power spectra. In Fig. 4, we show the average power spectrum of de-lensed CMB maps for a set of 100 simulations. The delensed power-spectra level is one order of magnitude larger than the residual noise and foreground-spectra level, which shows that the de-lensed maps are dominated by residual lensing signal.

To estimate the tensor-to-scalar ratio, we performed a likelihood analysis of the de-lensed B-mode power spectrum using Eq. (24). The model power spectrum is given in Eq. (27). The full covariance matrix was estimated from 400 purely Gaussian realisations of the foreground, allowing us to propagate the impact of residual foreground into the likelihood. The residual foreground power spectra, F_ℓ^res, in Eq. (27) was obtained by taking an average over 400 purely Gaussian realisations of Galactic-foreground maps. The steps followed to simulate the Gaussian foreground are explained in Sect. 3.1. To show the off-diagonal contributions of the covariance matrix, we present the corresponding correlation matrices for our three models in Fig. B.1. The model d10s5 shows the largest off-diagonal contribution in the correlation matrix.

We implemented the full pipeline on four sets of simulated maps: one set with foreground-free (FG-free) maps (see Sect. 3.1) and three sets each including foreground contamination based on the three models. All the maps in each set were passed through the HILC pipeline, lensing reconstruction, template building, and de-lensing pipeline before we performed the likelihood analysis. We present the mean estimate of the tensor-to-scalar ratio, r, obtained from a set of 100 simulations each for all four cases. The standard deviation of the mean shown in Fig. 8 represents the spread of the mean of the posterior of each simulation. The error on the mean, r, varies slightly between the models, with d12s7 showing the largest uncertainty due to its higher level of complexity. No significant deviation from the fiducial value is observed, demonstrating that the residual foreground biases in lensing reconstruction after component separation and correcting for the biases in the analysis do not have a significant impact on the estimation of r. The presence of residual foreground in SAT maps increases the uncertainty of the r measurement compared to the foreground-free case.

Fig. 8. Mean value of the tensor-to-scalar ratio, r, and its standard deviation, σ(r), estimated from 100 simulations for the three considered foreground models and FG-free maps. The vertical dashed line corresponds to the fiducial value of r.

4.5. Error budget on r

The additional B-mode components from residual instrumental noise and residual foreground contributes to the uncertainties in parameter estimation. To estimate contributions from different components to the error on r estimation, we used a Fisher matrix formalism. We considered the full-sky log-likelihood scaled by a factor of f_sky for the cut-sky power spectrum estimator with a given model power spectrum of Eq. (27). Then, we obtained the uncertainty on the r parameter from Fisher information matrix,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }^2(r)& ={\langle -\frac{{\partial }^{2}ln\mathcal{L}}{\partial r^2}\rangle }^{-1}\hfill \\ \hfill & =\frac{1}{f_{\mathrm{sky}}}{\displaystyle \sum _{{\ell }_b}}\frac{2}{(2{\ell }_{\mathrm{b}}+1)\mathrm{\Delta }\ell }\left(\frac{C_{{\ell }_{\mathrm{b}}}^{\mathrm{tens},r=1}}{r{C}_{{\ell }_{\mathrm{b}}}^{\mathrm{tens},r=1}+C_{{\ell }_{\mathrm{b}}}^{\mathrm{res}}+N_{{\ell }_{\mathrm{b}}}^{\mathrm{res}}+F_{{\ell }_{\mathrm{b}}}^{\mathrm{res}}}\right).\hfill \end{array}\end{array}$$(33)

Here, we used binned power spectra with a bin width of Δℓ. The summation is done over the multipoles at the centre, ℓ_b, of each bin. We considered the bins within 50 ≤ ℓ_b ≤ 300 with a width of Δℓ = 30.

In Fig. 9, we show the individual contributions to the total uncertainty in the tensor-to-scalar ratio, r, arising from different components present in the de-lensed B-mode power spectrum. The total uncertainty, σ(r), represents the expected statistical error due to the cosmic-variance limit on the measured de-lensed power spectrum for a single realisation. We obtained the uncertainty for each component of the model power spectrum in Eq. (33) by considering contribution from that particular component while setting all other components to zero. This approach allowed us to isolate the effect of each source of uncertainty.

Fig. 9. Contribution to σ(r) from different components of observed B modes in the presence of foreground (d10s5).

The tensor (r = 0.003) term in Fig. 9 reflects the fundamental limit imposed by the primordial B-mode signal itself, corresponding to the cosmic variance of the tensor component for a fiducial value of r = 0.003 considered in this paper. To quantify the impact of residual foreground in the template map used for de-lensing, we computed the difference in the variance of r, obtained using component-separated CMB maps with residual foregrounds and using analogous maps with subtracted residual foreground (foreground-free HILC CMB maps). This uncertainty level from foreground residue in template building – shown in Fig. 9 as foreground res. temp. – gives a quantitative measure of how residual foreground contamination in the B-mode template propagates into the constraint on r. The foreground res. and the noise res. terms arise from residual foreground and instrumental noise, respectively, which is present in the component-separated B-mode SAT maps. The lensing res. term captures the remaining B-mode power due the incomplete removal of lensing-induced B-modes. Finally, the combined uncertainty shows the total uncertainty when all contributions are included simultaneously, illustrating how each component cumulatively increases σ(r). Note that the combined uncertainty is not simply the sum of each individual contribution from different components.

5. Conclusions

In this work, we studied the impact of Galactic foreground and its residuals on CMB gravitational-lensing reconstruction, the de-lensing of CMB maps, and primordial gravitational-wave searches. Using simulations of Galactic foreground based on three different models, we quantified the mean-field effect, the bias in the reconstruction noise, and the non-Gaussian systematic bias. We demonstrated that the mean-field bias arises primarily from masking anisotropies, and that its amplitude for the EB quadratic estimator is comparable to the $N_L^{(1)}$ bias. We further showed that foreground affects the reconstruction of a lensing power spectrum in two different ways. The Gaussian component of the foreground primarily increases the reconstruction noise at intermediate and high multipoles, while non-Gaussian small-scale features induce a systematic bias at low multipoles.

We demonstrated that component separation is an essential step for both lensing power-spectrum reconstruction and tensor-to-scalar ratio estimation. In lensing reconstruction, component separation reduces the bias in $N_L^{(0)}$ reconstruction noise, as well as systematic biases from non-Gaussian small-scale features. As shown in Fig. 4, for the estimation of the tensor-to-scalar ratio, the de-lensed spectra are initially dominated by foreground contamination, where the cleanest frequency channels before component separation exhibit strong foreground emission, which is reduced by nearly two orders of magnitude after component separation. We see that the increase in reconstruction noise, $\mathrm{\Delta }{N}_L^{(0)}$, for the most complex model, d12s7, is comparable to a cosmic-variance limit at high multipoles. Therefore, we included the bias in the Wiener filtering when generating the template B-mode map in order to obtain optimal de-lensing efficiency and improved constraints on the tensor-to-scalar ratio. For the three considered foreground models, the estimated values of r from the likelihood analysis are consistent with the fiducial input value of r = 0.003. No significant bias in the mean estimated value of r is observed due to residual Galactic foreground. We observe only an ∼10% deviation in the case of the highest complexity model, d12s7, which could be due to the fact that the covariance matrix computed from purely Gaussian foreground does not fully capture the contribution of residual foreground in the delensed spectra.

The average uncertainty on the tensor-to-scalar ratio across the three models is about 70 − 80% of the fiducial value in our analysis, with the d12s7 model showing the largest uncertainty. We compare our results with those from Beck et al. (2020) and find approximately twice the level of uncertainty on the value of r. This difference can be attributed to the fact that we used the full Hamimeche and Lewis likelihood without the Gaussian approximation, thereby including off-diagonal correlations that tend to broaden the σ(r). In addition, we considered more realistic and updated models of the Galactic foreground, such as d10s5 and d12s7, which include higher levels of complexity and non-Gaussian structures on small angular scales. Furthermore, we employed a non-parametric component-separation technique, the Harmonic ILC method, which allowed us to avoid relying on prior assumptions about foreground.

In a similar work, Bianchini et al. (2025) investigated the impact of Galactic-foreground complexity on the estimation of the tensor-to-scalar ratio by comparing multiple component-separation pipelines applied to simulated CMB-S4–like data. Their analysis focused on biases and uncertainties in r arising from foreground cleaning and pipeline choices. In contrast, our work focused on the impact of non-Gaussian Galactic-foreground residues on CMB lensing reconstruction and the resulting constraints on r using a gradient-order, template-based de-lensing method. Bianchini et al. (2025) employed an iterative de-lensing procedure that removes nearly 90% of the lensing signal, yielding tighter constraints on r. Another key difference in the methodologies is that, rather than marginalising over residual foregrounds as in Bianchini et al. (2025), we included the residual foreground contribution directly in the covariance matrix and in the model power spectra. Together, these studies provide complementary perspectives on how foreground complexity propagates through different stages of CMB polarisation analyses.

As we show in Fig. 9, a 3σ detection of primordial gravitational waves with r ∼ 10⁻³ using a simple quadratic EB estimator in the de-lensing of the B-mode map is not achievable for CMB-S4–like experiments. The dominant contribution to the total uncertainty arises from the lensing residuals present in the de-lensed B-modes. Using a polarisation-based EB estimator, we were able to remove approximately 65% of the lensing-induced B-mode power. Finding optimal methods for de-lensing was beyond the of scope of our analysis; however, there are approaches enabling the improvement of the de-lensing efficiency. One such approach, especially relevant for small angular scales, is to combine the internally reconstructed lensing potential with other large-scale structure (LSS) tracers to generate the B-mode template used for de-lensing (Manzotti 2018; Hazumi et al. 2019; Ade et al. 2021; Hertig et al. 2024). One can also use iterative de-lensing methods (Carron & Lewis 2017; Belkner et al. 2024), removing up to about 90% of the lensing contamination. Such improvements would significantly reduce the contribution of lensing residuals to σ(r), potentially enabling a 3σ detection of tensor modes down to r ∼ 10⁻³. In such cases, residuals of the Galactic foreground, together with instrumental noise, provide a dominant contribution to the errors on the tensor-to-scalar ratio. For this reason, robust component separation and accurate correcting for the foreground residues for CMB lensing reconstruction analysed in our work will be highly relevant for the extraction of unbiased lensing potential and the detection of the primordial B-mode signal.

Acknowledgments

The authors acknowledge the use of the public software packages : CAMB, healpy, lenspyx, cmblensplus, plancklens, NaMaster, emcee, PySM3, Numpy, Scipy and matplotlib. We thank Radek Stompor, Dominic Beck, Carlo Baccigalupi and Anto I. Lonappan for helpful comments and suggestions on the results of this work. We thank Giuseppe Puglisi and Tuhin Ghosh for their suggestions on Galactic foreground modeling. We also thank Jacques Delabrouille, Julien Carron, Shamik Ghosh and Chandra Shekhar Saraf for their useful comments on delensing methods and on likelihood analysis techniques. KD acknowledges partial support from STER programme BPI/STE/2021/1/00033/U/00001 funded by NAWA. The authors are thankful for the computational resources provided by the CIŚ NCBJ cluster. This work used the NASA Astrophysics Data System Bibliographic Services.

References

Appendix A: HILC method implementation

Assuming multi-frequency foreground contaminated CMB observations in the N_c frequency bands and frequency maps convolved with an instrumental beam function, the observed signal for the i^th frequency band, where i ∈ 1, 2, …, N_c, in terms of harmonic coefficients of a given mode-index (ℓ,m), can be written as

$$\begin{array}{c}\hfill a_{\ell m}^{f,i}=B_{\ell }^i(a_{\ell m}^{c,i}+a_{\ell m}^s)+a_{\ell m}^{n,i},\end{array}$$(A.1)

where the coefficients $a_{\ell m}^{f,i}$, $a_{\ell m}^{c,i}$, a_ℓm^s and $a_{\ell m}^{n,i}$ denotes the total observed signal, net foreground emission, CMB signal and detector noise contribution, respectively. B_ℓⁱ is the beam transfer function of the instrument for i^th frequency band. Before applying HILC method to these coefficients, all the maps are first deconvolved with a beam of frequency-specific resolution and then convolved with common beam, generally to the highest available resolution amongst observed frequencies. Then, the cleaned CMB signal can be obtained as a linear combination of all N_c observed frequency maps in harmonic space,

$$\begin{array}{c}\hfill a_{\ell m}^{\mathit{Clean}}={\displaystyle \sum _{i=1}^{N_c}}{w}_{\ell }^i\frac{B_{\ell }^o}{B_{\ell }^i}{a}_{\ell m}^{f,i},\end{array}$$(A.2)

where B_ℓ^o is the common output beam function and w_ℓⁱ is the corresponding weight function of i^th band for a given multipole ℓ. The weights w_ℓⁱ are obtained by minimising the variance of the clean CMB signal under the constraint

$$\begin{array}{c}\hfill {\displaystyle \sum _i}{w}_i=1.\end{array}$$(A.3)

The constraint on the weights in Eq. (A.3) makes sure that the CMB signal is preserved. The minimum variance weight can be written as (Tegmark et al. 2003; Kim et al. 2009)

$$\begin{array}{c}\hfill w_{\ell }^i=\frac{{\sum }_j{[\widehat{}\}{C}_{\ell }^{-1}]}^{\mathit{ij}}}{{\sum }_{\mathit{ij}}{[\widehat{}\}{C}_{\ell }^{-1}]}^{\mathit{ij}}}.\end{array}$$(A.4)

Here, ${[\widehat{\mathbf{C}}}_{\ell }^{-1}{]}^{\mathit{ij}}$ denotes an element of the inverse of the covariance matrix having the i^th row index and j^th column index. Corresponding element of the covariance matrix is defined as

$$\begin{array}{c}\hfill {[{\widehat{C}}_{\ell }]}^{\mathit{ij}}=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{\ell }}{a}_{\ell m}^{f,i}{a}_{\ell m}^{f,j\ast }.\end{array}$$(A.5)

The HILC-cleaned maps still contains HILC noise and residual Galactic foreground signal. The residual noise and residual foreground maps are defined as the HILC weighted sum of noise and Galactic foreground contributions of each frequency channel. We estimate the residual noise and residual foreground and their power spectra as

$$\begin{array}{c}\hfill n_{\ell m}^{\mathit{res}}={\displaystyle \sum _i}{w}_{\ell }^i{a}_{\ell m}^{n,i},f_{\ell m}^{\mathit{res}}={\displaystyle \sum _i}{w}_{\ell }^i\left(\frac{B_{\ell }^o}{B_{\ell }^i}\right)a_{\ell m}^{f,i},\end{array}$$(A.6)

$$\begin{array}{c}\hfill N_{\ell }^{\mathit{res}}=1/(2\ell +1){\displaystyle \sum _m}|n_{\ell m}^{\mathit{res}}{|}^2,F_{\ell }^{\mathit{res}}=1/(2\ell +1){\displaystyle \sum _m}{|f_{\ell m}^{\mathit{res}}|}^2.\end{array}$$(A.7)

Here, $n_{\ell m}^{\mathit{res}}$ and $f_{\ell m}^{\mathit{res}}$ denote residual noise and foreground spherical harmonic coefficients, respectively, and N_ℓ^res and F_ℓ^res are corresponding estimators of the power spectra.

In Fig. A.1, we present the mean HILC weights and their standard deviation for the frequency channels in LAT configuration in the case of d10s5 Galactic foreground model computed from a set of 100 realisations. We see that, the cleanest channels 95 and 145 GHz has the highest contributions in the final HILC product. The weights for rest of the channels are closer to zero and do not contribute at high multipoles. The weights are estimated in a band of multipoles which is seen as step-like pattern in the plot.

Fig. A.1. Mean and standard deviation of HILC-derived weights over 100 realisation for d10s5 model in LAT configuration (fsky = 0.37).

The E- and B-mode HILC-cleaned LAT maps contain residual noise and residual foreground, which can be computed using Eqs. (A.6) and (A.7). In Fig. A.2 we compare the residual foreground for E and B modes from different foreground models in LAT configuration. Figure A.3 shows a comparison of the residual noise and residual Galactic foreground power spectra in the HILC B-mode products with the theoretical E- and B-mode input power spectra. The residual noise and residual foreground in the E-mode power spectra are of the same order of magnitude.

Fig. A.2. Residual foreground in HILC-cleaned E and B modes maps for one realisation in LAT configuration (fsky = 0.37).

Fig. A.3. Angular power spectra for CMB E and B-mode maps cleaned using the HILC method for the LAT configuration and sky coverage. The orange and blue solid zigzag lines are average power spectra over 100 simulations for the E and B-mode maps, respectively. The theoretical input power spectra are shown as black solid line. The map residual B-mode noise and foreground spectra for B modes are shown as dotted and dash-dotted lines, respectively.

Fig. A.4. HILC-cleaned CMB maps (upper row) and foreground residual maps (lower row) for SAT configuration and sky coverage for the Galactic foreground models d1s1 (left), (b) d10s5 (middle) and d12s7 (right). These cleaned maps (upper row) includes residual HILC noise and foreground defined by Eq. (A.6).

Appendix B: Covariance matrix for likelihood analysis

The covariance matrices defined by Eq. 26 captures non-zero off-diagonal contributions for the three foreground models considered. We computed correlation matrix for the covariance matrix shown in Fig. B.1 and we see that off-diagonal terms are non-negligible. Therefore, we considered the inverse of full covariance matrix in the likelihood analysis to constraint the tensor-to-scalar ratio.

Fig. B.1. Correlation matrices for the covariance matrix in Eq. (26) for the three foreground models. The full covariance matrices have off-diagonal contributions.

Appendix C: Delensed B-mode maps

Delensed B-mode maps includes different components of B-mode signal as mentioned in the power spectrum model in Eq. 27. In Fig. C.1, we show each different component separately to compare their level of contribution to the delensed B-mode map.

Fig. C.1. Different components of the CMB B-mode maps for SAT configuration and sky coverage. All the maps only includes multipoles in the range 30 ≤ ℓ ≤ 300. Observed B-mode map in top-left panel includes residue after component separation noise and foreground for the d10s5 model. The remaining maps are noise-free and foreground-free.

All Tables

All Figures

	Fig. 1. Analytical zeroth-order lensing reconstruction bias, $N_L^{(0)}$, for different estimators given by Eq. (14). The instrumental specification was considered for a CMB-S4-like experiment with a beam resolution of 2.5 arcminutes and white-noise temperature level of ΔT = 2 μK arcmin and polarisation of ΔP = 2.8 μK arcmin.
In the text

Fig. 2. Masks used in our analysis to mimic CMB-S4 sky coverage, shown in Ecliptic coordinates with Galactic emission at 145 GHz in the background. (a) The left panel shows the mask used for component separation for the wide sky survey corresponding to the LAT configuration (fsky = 0.64). (b) The middle panel shows the mask used for lensing reconstruction for the wide sky survey (fsky = 0.37). (c) The right panel shows the sky area considered for the deep survey corresponding to the SAT configuration, which was used for the delensing study (fsky = 0.025).

In the text

	Fig. 3. Beam-de-convolved noise power spectra for temperature (solid) and polarisation (dashed) in LAT configuration.
In the text

Fig. 4. Angular power spectra for CMB B-mode-cleaned maps for the SAT configuration and sky coverage. The solid grey zigzag lines are BB power spectra of 100 simulations, and the average is shown as red stars. The theoretical input power spectra are shown as solid black lines. The residual HILC noise and foreground levels are shown as dotted and dashed lines, respectively. The solid magenta line shows the average de-lensed B-mode spectra over 100 simulations. The dashed blue line shows the average power spectra of Galactic-foreground B modes at the frequency channels 95 and 145 GHz. The dashed grey line shows the tensor B-mode level for r = 0.003.

In the text

Fig. 5. Angular power spectra of reconstructed lensing field shown as a solid blue line. The mean field effect on the EB estimator for MF maps for d10s5 and HILC-cleaned maps for d10s5 and d12s7 are shown as solid red, green, and grey lines, respectively. The mean-field effect on the EE estimator for HILC-cleaned maps of d10s5 is shown as a solid cyan line. Theoretical lensing power spectra are shown as solid black lines. The dashed magenta line is the sum of theory power spectra and lensing reconstruction noise terms, N(0) and N(1), shown as dashed (orange) and dash-dotted (purple) lines, respectively.

In the text

	Fig. 6. Relative bias in semi-analytic $N_L^{(0)}$ noise due to foreground contamination in MF channels (solid lines) before component separation and residual foreground in HILC maps (dashed lines) after component separation for different models.
In the text

Fig. 7. The $F_L^{\mathrm{syst}}$ bias is shown for MF maps (dark colours) and HILC-cleaned CMB maps (light colours) for three different models. The black curve shows the fiducial lensing power spectrum, and the dashed red curve is the sum of reconstruction noises, N(0) + N(1). The dash-dotted magenta line is the bias in reconstruction noise, $\mathrm{\Delta }{N}_L^{(0)}$, due to residual Galactic foreground after component separation for the d12s7 model. The shaded grey region represents the cosmic variance limit on the reconstructed lensing power spectrum.

In the text

	Fig. 8. Mean value of the tensor-to-scalar ratio, r, and its standard deviation, σ(r), estimated from 100 simulations for the three considered foreground models and FG-free maps. The vertical dashed line corresponds to the fiducial value of r.
In the text

	Fig. 9. Contribution to σ(r) from different components of observed B modes in the presence of foreground (d10s5).
In the text

	Fig. A.1. Mean and standard deviation of HILC-derived weights over 100 realisation for d10s5 model in LAT configuration (fsky = 0.37).
In the text

	Fig. A.2. Residual foreground in HILC-cleaned E and B modes maps for one realisation in LAT configuration (fsky = 0.37).
In the text

Fig. A.3. Angular power spectra for CMB E and B-mode maps cleaned using the HILC method for the LAT configuration and sky coverage. The orange and blue solid zigzag lines are average power spectra over 100 simulations for the E and B-mode maps, respectively. The theoretical input power spectra are shown as black solid line. The map residual B-mode noise and foreground spectra for B modes are shown as dotted and dash-dotted lines, respectively.

In the text

	Fig. A.4. HILC-cleaned CMB maps (upper row) and foreground residual maps (lower row) for SAT configuration and sky coverage for the Galactic foreground models d1s1 (left), (b) d10s5 (middle) and d12s7 (right). These cleaned maps (upper row) includes residual HILC noise and foreground defined by Eq. (A.6).
In the text

	Fig. B.1. Correlation matrices for the covariance matrix in Eq. (26) for the three foreground models. The full covariance matrices have off-diagonal contributions.
In the text

	Fig. C.1. Different components of the CMB B-mode maps for SAT configuration and sky coverage. All the maps only includes multipoles in the range 30 ≤ ℓ ≤ 300. Observed B-mode map in top-left panel includes residue after component separation noise and foreground for the d10s5 model. The remaining maps are noise-free and foreground-free.
In the text