BICEP/Keck XX: Component-separated maps of polarized CMB and thermal dust emission using Planck and BICEP/Keck Observations through the 2018 Observing Season

We present component-separated polarization maps of the cosmic microwave background (CMB) and Galactic thermal dust emission, derived using data from the BICEP/Keck experiments through the 2018 observing season and Planck. By employing a maximum-likelihood method that utilizes observing matrices, we produce unbiased maps of the CMB and dust signals. We outline the computational challenges and demonstrate an efficient implementation of the component map estimator. We show methods to compute and characterize power spectra of these maps, opening up an alternative way to infer the tensor-to-scalar ratio from our data. We compare the results of this map-based separation method with the baseline BICEP/Keck analysis. Our analysis demonstrates consistency between the two methods, finding an 84% correlation between the pipelines.

💡 Research Summary

This paper presents a new component‑separation pipeline that produces unbiased, maximum‑likelihood maps of the polarized Cosmic Microwave Background (CMB) and Galactic thermal dust emission by jointly analysing data from the BICEP/Keck experiments (through the 2018 observing season) and the Planck satellite. The authors formulate the problem in a linear observation model: each frequency‑specific Q/U map is generated by applying an observing matrix R (which encodes timestream filtering, de‑projection, and sky‑cut effects) and a beam‑convolution operator B to the true sky signal, then adding Gaussian noise n. The sky signal s is modeled as a linear combination of two components—CMB and dust—scaled across frequencies by a parametric modified black‑body law with a fixed dust spectral index β_d and temperature T_d taken from the baseline BK18 analysis.

The likelihood for the full data vector d (the stacked set of all frequency maps) is Gaussian with covariance given by the noise matrix N̂, which the authors approximate as diagonal in pixel space (off‑diagonal correlations are <5 %). Maximising the likelihood with respect to the component maps yields the generalized least‑squares solution

  ŝ = (Aᵀ N̂⁻¹ A)⁻¹ Aᵀ N̂⁻¹ d,

where A = R B F combines the observing matrix, beam operator, and frequency‑scaling matrix. Direct inversion of the massive matrix (Aᵀ N̂⁻¹ A) is computationally infeasible because its size is 4 n_pix × 4 n_pix (n_pix ≈ 10⁵). To overcome this, the authors implement an iterative solver (preconditioned conjugate‑gradient) that never forms the full matrix explicitly; only matrix‑vector products are required. A block‑Jacobi preconditioner is constructed by dropping R and B, i.e., using (Fᵀ N̂⁻¹ F)⁻¹, which captures the dominant diagonal structure and accelerates convergence dramatically.

The sparsity pattern of R is examined in detail: each frequency‑Q/U block contains roughly 5 × 10⁹ non‑zero entries, reflecting the complex filtering applied to the BICEP/Keck data. Nevertheless, because R, B, and F are either sparse in pixel space or diagonal in harmonic space, fast Fourier‑ and spherical‑harmonic transforms enable rapid application of these operators. External Planck maps at 100, 143, 217, and 353 GHz are incorporated without re‑observing, extending the frequency coverage and stabilising the linear system by filling null‑space modes that are poorly constrained by the BICEP/Keck data alone.

Data used include the standard BK18 maps: BICEP2 (150 GHz), Keck Array (95, 150, 220 GHz), and BICEP3 (95 GHz), all filtered, polynomial‑subtracted, and binned into 0.25° equirectangular pixels over the South‑Pole patch (≈ 400–600 deg²). Planck NPIPE Q/U maps are taken directly. The noise model for each map is derived from inverse‑variance weighting of the timestream data, yielding per‑pixel variance maps that populate the diagonal of N̂.

Applying the maximum‑likelihood estimator produces CMB and dust Q/U maps. Power‑spectrum analysis of these maps shows excellent agreement with the BK18 baseline analysis, which fits multi‑frequency auto‑ and cross‑spectra using a parametric foreground model. Pixel‑level correlation between the two pipelines is 0.84, indicating that the map‑based separation recovers essentially the same sky signal as the spectrum‑based approach. The authors propagate the uncertainty in β_d by using the distribution of best‑fit β_d values obtained from 499 simulated skies (as shown in Fig. 1 of the paper) and inserting the appropriate β_d for each simulation into the estimator, thereby quantifying the impact of dust‑spectral‑index uncertainty on the final maps.

The paper discusses limitations and future extensions. Currently β_d is held fixed and spatially uniform; allowing β_d (and the synchrotron index β_s) to vary across the sky would require a two‑step map‑level fitting or a fully non‑parametric maximum‑likelihood framework, as explored in recent literature. Incorporating such spatial variability would improve robustness against complex foregrounds and reduce potential biases in the tensor‑to‑scalar ratio r. The authors also note that the iterative solver scales well with additional frequencies and higher‑resolution maps, suggesting that the method can be applied to upcoming experiments (e.g., BICEP Array, Simons Observatory) where the number of frequency channels and pixel count will increase dramatically.

In summary, the work demonstrates that a maximum‑likelihood, map‑level component separation using realistic observing matrices is computationally tractable and yields results consistent with the established BK18 spectral‑likelihood pipeline. This provides an independent validation of the BICEP/Keck r constraints and opens a pathway for future analyses that may benefit from direct map‑level foreground cleaning, especially when dealing with increasingly complex data sets.