# Likelihood methods for CMB experiments

Martina Gerbino, Massimiliano Lattanzi, Marina Migliaccio, Luca Pagano, Laura Salvati, Loris Colombo, Alessandro Gruppuso, Paolo Natoli, Gianluca Polenta

► To cite this version:

Martina Gerbino, Massimiliano Lattanzi, Marina Migliaccio, Luca Pagano, Laura Salvati, et al. Likelihood methods for CMB experiments. Frontiers in Physics, 2020, 8, pp.15. 10.3389/fphy.2020.00015. hal-02372452

HAL Id: hal-02372452

https://hal.science/hal-02372452

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

Check for
updates

# Likelihood Methods for CMB Experiments

Martina Gerbino<sup>1,2*</sup>, Massimiliano Lattanzi<sup>1*</sup>, Marina Migliaccio<sup>3,4*</sup>, Luca Pagano<sup>1,5*</sup>, Laura Salvati<sup>6,7,8*</sup>, Loris Colombo<sup>9</sup>, Alessandro Gruppuso<sup>10,11</sup>, Paolo Natoli<sup>1,5</sup> and Gianluca Polenta<sup>12</sup>

<sup>1</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Ferrara, Italy, <sup>2</sup> High Energy Physics (HEP) Division, Argonne National Laboratory, Lemont, IL, United States, <sup>3</sup> Dipartimento di Fisica, Università di Roma Tor Vergata, Rome, Italy, <sup>4</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Roma 2, Rome, Italy, <sup>5</sup> Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Ferrara, Italy, <sup>6</sup> Institute for Fundamental Physics of the Universe, Trieste, Italy, <sup>7</sup> INAF - Osservatorio Astronomico di Trieste, Trieste, Italy, <sup>8</sup> Institut d'Astrophysique Spatiale, CNRS (UMR 8617) Université Paris-Sud, Orsay, France, <sup>9</sup> Dipartimento di Fisica, Università degli Studi di Milano, Milan, Italy, <sup>10</sup> INAF-OAS Bologna, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Istituto Nazionale di Astrofisica, Bologna, Italy, <sup>11</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Bologna, Italy, <sup>12</sup> Space Science Data Center - Agenzia Spaziale Italiana, Rome, Italy

**OPEN ACCESS**

**Edited by:**

Michele Liguori,
University of Padova, Italy

**Reviewed by:**

Sayantan Choudhury,
Max-Planck-Institut für
Gravitationsphysik, Germany

Andrew H. Jaffe,
Imperial College London,
United Kingdom

*Correspondence:

Martina Gerbino
gerbinom@fe.infn.it

Massimiliano Lattanzi
lattanzi@fe.infn.it

Marina Migliaccio
marina.migliaccio@roma2.infn.it

Luca Pagano
pgnlcu@unife.it

Laura Salvati
laura.salvati@inaf.it

**Specialty section:**

This article was submitted to
Cosmology,
a section of the journal
Frontiers in Physics

**Received:** 27 May 2019

**Accepted:** 16 January 2020

**Published:** 26 February 2020

**Citation:**

Gerbino M, Lattanzi M, Migliaccio M, Pagano L, Salvati L, Colombo L, Gruppuso A, Natoli P and Polenta G (2020) Likelihood Methods for CMB Experiments. Front. Phys. 8:15. doi: 10.3389/fphy.2020.00015

A great deal of experimental effort is currently being devoted to the precise measurements of the cosmic microwave background (CMB) sky in temperature and polarization. Satellites, balloon-borne, and ground-based experiments scrutinize the CMB sky at multiple scales, and therefore enable to investigate not only the evolution of the early Universe, but also its late-time physics with unprecedented accuracy. The pipeline leading from time ordered data as collected by the instrument to the final product is highly structured. Moreover, it has also to provide accurate estimates of statistical and systematic uncertainties connected to the specific experiment. In this paper, we review likelihood approaches targeted to the analysis of the CMB signal at different scales, and to the estimation of key cosmological parameters. We consider methods that analyze the data in the spatial (i.e., pixel-based) or harmonic domain. We highlight the most relevant aspects of each approach and compare their performance.

**Keywords:** CMB, likelihood function, data analysis - methods, parameter estimation, CMB statistics

## 1. INTRODUCTION

After 71 years from its first predictions, and after 55 years from its first observational evidences, the cosmic microwave background (CMB) is nowadays one of the most important probes in cosmology. During the past decades, theoretical efforts have elucidated the physics leading to the pattern of anisotropies in temperature and polarization (see e.g., [1] for a review on CMB physics and [2] for an exhaustive review on early Universe physics). Quantum fluctuations in the early universe generate metric perturbations. Scalar perturbations are converted into matter perturbations and radiation anisotropies that evolve in the expanding universe according to a set of coupled Einstein, Boltzmann and fluid equations. Matter perturbations eventually grow into galaxies and galaxy clusters. Primary CMB anisotropies are frozen at the time of matter-radiation decoupling, and subsequently modified during the propagation through evolving structures from the last scattering surface to the observer. Scattering between free electrons and CMB photons in two distinct epochs (recombination and reionization) further enriches the CMB structure with the addition of a polarization “curl-free” (E-mode) pattern in the CMB radiation. Gravitational lensing of the CMB due to the propagation of CMB photons throughout large-scale structures (LSS)
<figure><img src="image_2.png" /></figure>

Check for
updates

# Likelihood Methods for CMB Experiments

Martina Gerbino<sup>1,2*</sup>, Massimiliano Lattanzi<sup>1*</sup>, Marina Migliaccio<sup>3,4*</sup>, Luca Pagano<sup>1,5*</sup>, Laura Salvati<sup>6,7,8*</sup>, Loris Colombo<sup>9</sup>, Alessandro Gruppuso<sup>10,11</sup>, Paolo Natoli<sup>1,5</sup> and Gianluca Polenta<sup>12</sup>

<sup>1</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Ferrara, Italy, <sup>2</sup> High Energy Physics (HEP) Division, Argonne National Laboratory, Lemont, IL, United States, <sup>3</sup> Dipartimento di Fisica, Università di Roma Tor Vergata, Rome, Italy, <sup>4</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Roma 2, Rome, Italy, <sup>5</sup> Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Ferrara, Italy, <sup>6</sup> Institute for Fundamental Physics of the Universe, Trieste, Italy, <sup>7</sup> INAF - Osservatorio Astronomico di Trieste, Trieste, Italy, <sup>8</sup> Institut d'Astrophysique Spatiale, CNRS (UMR 8617) Université Paris-Sud, Orsay, France, <sup>9</sup> Dipartimento di Fisica, Università degli Studi di Milano, Milan, Italy, <sup>10</sup> INAF-OAS Bologna, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Istituto Nazionale di Astrofisica, Bologna, Italy, <sup>11</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Bologna, Italy, <sup>12</sup> Space Science Data Center - Agenzia Spaziale Italiana, Rome, Italy

## OPEN ACCESS

**Edited by:**
Michele Liguori,
University of Padova, Italy

**Reviewed by:**
Sayantan Choudhury,
Max-Planck-Institut für
Gravitationsphysik, Germany
Andrew H. Jaffe,
Imperial College London,
United Kingdom

**\*Correspondence:**
Martina Gerbino
gerbinom@fe.infn.it
Massimiliano Lattanzi
lattanzi@fe.infn.it
Marina Migliaccio
marina.migliaccio@roma2.infn.it
Luca Pagano
pgnlcu@unife.it
Laura Salvati
laura.salvati@inaf.it

**Specialty section:**
This article was submitted to
Cosmology,
a section of the journal
Frontiers in Physics

**Received:** 27 May 2019
**Accepted:** 16 January 2020
**Published:** 26 February 2020

**Citation:**
Gerbino M, Lattanzi M, Migliaccio M, Pagano L, Salvati L, Colombo L, Gruppuso A, Natoli P and Polenta G (2020) Likelihood Methods for CMB Experiments. Front. Phys. 8:15. doi: 10.3389/fphy.2020.00015

A great deal of experimental effort is currently being devoted to the precise measurements of the cosmic microwave background (CMB) sky in temperature and polarization. Satellites, balloon-borne, and ground-based experiments scrutinize the CMB sky at multiple scales, and therefore enable to investigate not only the evolution of the early Universe, but also its late-time physics with unprecedented accuracy. The pipeline leading from time ordered data as collected by the instrument to the final product is highly structured. Moreover, it has also to provide accurate estimates of statistical and systematic uncertainties connected to the specific experiment. In this paper, we review likelihood approaches targeted to the analysis of the CMB signal at different scales, and to the estimation of key cosmological parameters. We consider methods that analyze the data in the spatial (i.e., pixel-based) or harmonic domain. We highlight the most relevant aspects of each approach and compare their performance.

**Keywords:** CMB, likelihood function, data analysis - methods, parameter estimation, CMB statistics

## 1. INTRODUCTION

After 71 years from its first predictions, and after 55 years from its first observational evidences, the cosmic microwave background (CMB) is nowadays one of the most important probes in cosmology. During the past decades, theoretical efforts have elucidated the physics leading to the pattern of anisotropies in temperature and polarization (see e.g., [1] for a review on CMB physics and [2] for an exhaustive review on early Universe physics). Quantum fluctuations in the early universe generate metric perturbations. Scalar perturbations are converted into matter perturbations and radiation anisotropies that evolve in the expanding universe according to a set of coupled Einstein, Boltzmann and fluid equations. Matter perturbations eventually grow into galaxies and galaxy clusters. Primary CMB anisotropies are frozen at the time of matter-radiation decoupling, and subsequently modified during the propagation through evolving structures from the last scattering surface to the observer. Scattering between free electrons and CMB photons in two distinct epochs (recombination and reionization) further enriches the CMB structure with the addition of a polarization “curl-free” (E-mode) pattern in the CMB radiation. Gravitational lensing of the CMB due to the propagation of CMB photons throughout large-scale structures (LSS)
generate a polarized “divergence-free” (B-mode) patter from the distortion of the CMB E-modes. Perhaps more elusive, though of paramount importance, is the primordial B-mode signal sourced by tensor perturbations to the metric (gravitational waves).

On the other hand, observational efforts have progressively lead the field to the current stage of precision cosmology. Observations of the CMB sky from space missions [3–5], balloon-borne experiments [6–8], and from ground-based telescopes [9–13] provided measurements of CMB anisotropies in temperature and polarization over a wide range of angular scales. While we are writing this review, the Planck collaboration [14] is preparing the final public release of the Planck legacy products, which will likely represent the state of the art of CMB measurements from a single experiment for the next decade and more. Current observations [3] are in agreement with the standard cosmological model of a homogeneous and isotropic Universe at large scales, based on General Relativity and on the standard model (SM) of particle physics, complemented with a mechanism for the generation of primordial perturbations, i.e., the inflationary paradigm. When interpreted in this ΛCDM framework, cosmological data point to a spatially flat Universe composed by baryons ($\Omega_b h^2 = 0.02237 \pm 0.00015$, ~ 5% of the total density), dark matter ($\Omega_c h^2 = 0.1200 \pm 0.0012$, ~ 25%), and dark energy ($\Omega_\Lambda = 0.6847 \pm 0.0073$, ~ 70%), a component that behaves like a cosmological constant, and is responsible for the present accelerated expansion, plus photons (a few parts in $10^5$) and light neutrinos. Further advances in CMB observations are still to come. Planned upgrades of existing ground-based and balloon experiments are ongoing [15–18]. The next generation of CMB observatories is under construction and is paving the way to the “stage IV” experiments targeting the ultimate measurements of the CMB polarization field [19–23].

The long run that lead from the pivotal observations of Penzias and Wilson [24] to the Planck legacy release has seen the dramatic improvement of the sensitivity to key cosmological parameters. Planck 2018 data provides sub-percent constraints on the base-ΛCDM parameters<fnref n="1" /> [3]. Moreover, advances in experimental cosmology over the past decades made cosmology itself a new avenue to the investigation of fundamental physics properties complementary to laboratory searches. A clear example is given by the possibility to constrain neutrino properties, such as their number $N_{\text{eff}}$ and the sum of their masses $\sum m_\nu$. Indeed, the combination of Planck 2018 data and LSS information (in the form of measurements of baryon acoustic oscillations, BAO) can exclude at 95% c.l. the presence of light thermal relics decoupling after the QCD phase transition ($T < 100$ MeV) and provides a bound on the sum of the neutrino masses of $\sum m_\nu < 0.12$ eV at 95% c.l.<fnref n="2" /> [3].

In this context, a key ingredient is the suitable choice of the likelihood function to compare observed data with theoretical

predictions in order to constrain the model parameters. In the standard cosmological model of the early universe, primordial perturbations are Gaussian distributed, and so are CMB fluctuations. Therefore, all relevant physical information in the CMB field are contained in the variance of the distribution. This is the reason why the full-sky power spectra of CMB fluctuations are a sufficient statistics. The power spectra of observed data also provide an unbiased estimator of the ensemble averaged variance of the CMB fluctuations. In the simple case of full-sky observations and isotropic and mode-uncorrelated experimental noise, the likelihood function can be derived analytically. In particular, for correlated temperature and polarization field, the probability of the data given the theoretical model (i.e., the likelihood $\mathcal{L}$) is given by a Wishart distribution.

However, this simple case does not capture the properties of realistic observations. Depending on the experimental platform (satellite, balloon, ground), each telescope has access to fractions of the sky $f_{\text{sky}}$ of different size. As an example, compare the almost full-sky observations of the Planck satellite [14] with the $f_{\text{sky}} \sim 1\%$ sky coverage of the ground-based BICEP-2 experiment [11]. Even in the case of full-sky observations, only a certain fraction of the sky can be retained for cosmological analysis. Foreground emissions from astrophysical and galactic sources should be masked if particularly bright contaminants. In addition to limited access to the sky coverage, a particular choice of the observational (or scanning) strategies of the sky can break the assumption of isotropic noise, due to repeated visits to the same part of the sky. As an example, consider the Planck scanning strategy featuring a longer integration time in the proximity of the Ecliptic poles (i.e., at lower Galactic latitudes, where galactic foreground contaminations are smaller).

In general, complications to the simple case of full-sky and isotropic noise arising from realistic experimental conditions require a different likelihood analysis. First of all, specific estimators of the power spectra should be defined in the partial-sky regime, which take into account spurious correlations between fields induced by the incomplete sky coverage. Secondly, the use of a Wishart distribution as a likelihood function is no longer possible. Either the new estimators are no longer distributed according to a Wishart, and therefore this choice is not exact anymore. Or, the use of the exact likelihood is unfeasible as one moves to the analysis of smaller scales (larger multipoles) and higher-resolution maps, due to the huge computational cost of inverting large covariance matrices. At large scales and for low-enough angular resolutions, the exact likelihood in pixel space can still be adopted.

In all the above situations, approximate forms of the likelihood functions have been developed. At small scales, the central limit theorem allows to approximate the Wishart distribution as a Gaussian in the power spectra. In general, quadratic forms in some functions of the CMB spectra have been adopted as approximate likelihood functions, with various choices of the covariance matrix. To conclude this long introduction, it has to be stressed that the choice of the likelihood strongly depends on the characteristics of the experiment at large, i.e., on the observational strategy, on the range of scales probed, on the noise properties, etc.

<footnote n="1">¹The base ΛCDM parameters are: the angular size of the sound horizon at recombination $\theta_s$, the amplitude $A_s$ and tilt $n_s$ of the spectrum of primordial scalar perturbations, the reionization optical depth $\tau$, the energy density in baryonic matter $\Omega_b h^2$ and in cold dark matter $\Omega_c h^2$.</footnote>

<footnote n="2">²Constraints derived in the context of minimal extensions of the standard ΛCDM model.</footnote>
generate a polarized “divergence-free” (B-mode) patter from the distortion of the CMB E-modes. Perhaps more elusive, though of paramount importance, is the primordial B-mode signal sourced by tensor perturbations to the metric (gravitational waves).

On the other hand, observational efforts have progressively lead the field to the current stage of precision cosmology. Observations of the CMB sky from space missions [3–5], balloon-borne experiments [6–8], and from ground-based telescopes [9–13] provided measurements of CMB anisotropies in temperature and polarization over a wide range of angular scales. While we are writing this review, the Planck collaboration [14] is preparing the final public release of the Planck legacy products, which will likely represent the state of the art of CMB measurements from a single experiment for the next decade and more. Current observations [3] are in agreement with the standard cosmological model of a homogeneous and isotropic Universe at large scales, based on General Relativity and on the standard model (SM) of particle physics, complemented with a mechanism for the generation of primordial perturbations, i.e., the inflationary paradigm. When interpreted in this ΛCDM framework, cosmological data point to a spatially flat Universe composed by baryons ($\Omega_b h^2 = 0.02237 \pm 0.00015$, ~ 5% of the total density), dark matter ($\Omega_c h^2 = 0.1200 \pm 0.0012$, ~ 25%), and dark energy ($\Omega_\Lambda = 0.6847 \pm 0.0073$, ~ 70%), a component that behaves like a cosmological constant, and is responsible for the present accelerated expansion, plus photons (a few parts in $10^5$) and light neutrinos. Further advances in CMB observations are still to come. Planned upgrades of existing ground-based and balloon experiments are ongoing [15–18]. The next generation of CMB observatories is under construction and is paving the way to the “stage IV” experiments targeting the ultimate measurements of the CMB polarization field [19–23].

The long run that lead from the pivotal observations of Penzias and Wilson [24] to the Planck legacy release has seen the dramatic improvement of the sensitivity to key cosmological parameters. Planck 2018 data provides sub-percent constraints on the base-ΛCDM parameters<fnref n="1" /> [3]. Moreover, advances in experimental cosmology over the past decades made cosmology itself a new avenue to the investigation of fundamental physics properties complementary to laboratory searches. A clear example is given by the possibility to constrain neutrino properties, such as their number $N_{\text{eff}}$ and the sum of their masses $\sum m_\nu$. Indeed, the combination of Planck 2018 data and LSS information (in the form of measurements of baryon acoustic oscillations, BAO) can exclude at 95% c.l. the presence of light thermal relics decoupling after the QCD phase transition ($T < 100$ MeV) and provides a bound on the sum of the neutrino masses of $\sum m_\nu < 0.12$ eV at 95% c.l.<fnref n="2" /> [3].

In this context, a key ingredient is the suitable choice of the likelihood function to compare observed data with theoretical

predictions in order to constrain the model parameters. In the standard cosmological model of the early universe, primordial perturbations are Gaussian distributed, and so are CMB fluctuations. Therefore, all relevant physical information in the CMB field are contained in the variance of the distribution. This is the reason why the full-sky power spectra of CMB fluctuations are a sufficient statistics. The power spectra of observed data also provide an unbiased estimator of the ensemble averaged variance of the CMB fluctuations. In the simple case of full-sky observations and isotropic and mode-uncorrelated experimental noise, the likelihood function can be derived analytically. In particular, for correlated temperature and polarization field, the probability of the data given the theoretical model (i.e., the likelihood $\mathcal{L}$) is given by a Wishart distribution.

However, this simple case does not capture the properties of realistic observations. Depending on the experimental platform (satellite, balloon, ground), each telescope has access to fractions of the sky $f_{\text{sky}}$ of different size. As an example, compare the almost full-sky observations of the Planck satellite [14] with the $f_{\text{sky}} \sim 1\%$ sky coverage of the ground-based BICEP-2 experiment [11]. Even in the case of full-sky observations, only a certain fraction of the sky can be retained for cosmological analysis. Foreground emissions from astrophysical and galactic sources should be masked if particularly bright contaminants. In addition to limited access to the sky coverage, a particular choice of the observational (or scanning) strategies of the sky can break the assumption of isotropic noise, due to repeated visits to the same part of the sky. As an example, consider the Planck scanning strategy featuring a longer integration time in the proximity of the Ecliptic poles (i.e., at lower Galactic latitudes, where galactic foreground contaminations are smaller).

In general, complications to the simple case of full-sky and isotropic noise arising from realistic experimental conditions require a different likelihood analysis. First of all, specific estimators of the power spectra should be defined in the partial-sky regime, which take into account spurious correlations between fields induced by the incomplete sky coverage. Secondly, the use of a Wishart distribution as a likelihood function is no longer possible. Either the new estimators are no longer distributed according to a Wishart, and therefore this choice is not exact anymore. Or, the use of the exact likelihood is unfeasible as one moves to the analysis of smaller scales (larger multipoles) and higher-resolution maps, due to the huge computational cost of inverting large covariance matrices. At large scales and for low-enough angular resolutions, the exact likelihood in pixel space can still be adopted.

In all the above situations, approximate forms of the likelihood functions have been developed. At small scales, the central limit theorem allows to approximate the Wishart distribution as a Gaussian in the power spectra. In general, quadratic forms in some functions of the CMB spectra have been adopted as approximate likelihood functions, with various choices of the covariance matrix. To conclude this long introduction, it has to be stressed that the choice of the likelihood strongly depends on the characteristics of the experiment at large, i.e., on the observational strategy, on the range of scales probed, on the noise properties, etc.

<footnote n="1">¹The base ΛCDM parameters are: the angular size of the sound horizon at recombination $\theta_*$, the amplitude $A_S$ and tilt $n_S$ of the spectrum of primordial scalar perturbations, the reionization optical depth $\tau$, the energy density in baryonic matter $\Omega_b h^2$ and in cold dark matter $\Omega_c h^2$.</footnote>

<footnote n="2">²Constraints derived in the context of minimal extensions of the standard ΛCDM model.</footnote>
The aim of this manuscript is to review the basics of likelihood analysis in CMB experiments. This goal is motivated by the fact the we are at a crucial point in the history of observational CMB cosmology. The level of maturity and complexity reached by current CMB experiments boosted the theoretical efforts in finding smart solutions to the issue of identifying a suitable likelihood choice. A rich literature has been produced in this sense, although an exhaustive overview of the topic is not available, to the best of our knowledge. This review would fill the gap. Such a review could also serve as a good starting point for those who are approaching the field of CMB data analysis today or in the next future, and would be ideally contributing to the advances of CMB science in the next decades.

The structure of the manuscript is as follows. Section 2 is devoted to the statistics of the CMB fields. The approach to the topic is pedagogical, in a sense that we begin with a discussion in the single-field, temperature-only regime and introduce the basic statistical properties of CMB fluctuations. Then, we move to the more general case of correlated T, E, B fields. The discussion is carried over in the full-sky regime, with no distinction made between applications to large- and small-scales. We conclude section 2 with the introduction of the exact likelihood in full sky. Specific approximations to the exact likelihood are presented in section 3 (applications to the small-scale regime) and in section 5 (applications to the large-scale regime). In both sections, attention is devoted to complications due to partial-sky coverage and noise contamination. The inclusion of physical late-time Universe effects on the CMB photons in terms of gravitational lensing is detailed in section 4, whereas the important issues related to the presence of foreground emissions are described in section 6. The discussion of the various likelihood approaches in terms of computational cost (where applicable) and robustness with respect to the ability to provide unbiased estimations of cosmological parameters is detailed in section 7. Our conclusions are summarized in section 8. Some useful tools that will be mentioned throughout the main text are further discussed in Appendix A. In particular, in Appendix A, we review the basic notions of statistics needed to develop the formalism of CMB statistics. In Appendix B, we discuss the construction of power spectrum estimators, including pseudo-C_l, the “pure” formalism, and quadratic maximum likelihood (QML) estimator.

## 2. STATISTICS OF THE COSMIC MICROWAVE BACKGROUND

We now introduce some basic aspects of the statistics of the CMB. The basic object that we are interested in is the likelihood function L, i.e., the probability of the observed data d given a model, regarded as a function of the model itself. If the model is defined in terms of a vector of parameters θ, we thus have:

$$ \mathcal{L}(\theta) = p(d|\theta). \qquad (1) $$

The notation used throughout this review is presented in Appendix A, were we also recall some basic notion of probability and statistics.

We first derive the exact likelihood function for the CMB fields in harmonic space. The main statistical concepts are introduced in the limit of single field (section 2.1), i.e., temperature only, for the sake of simplicity. We then generalize these main findings in the case of joint temperature and polarization analysis (section 2.2). The exact likelihood in real space are derived in section 2.3.

We assume an ideal scenario of full-sky observations with infinite angular resolution and absence of noise and foreground contaminations. Obviously, this scenario is highly idealized. Nevertheless, it allows to easily derive the basic concepts of CMB statistics. Modifications to this picture arising from realistic observational issues (limited sky coverage, masked sky, experimental noise, finite angular resolution) are introduced in section 3.3. Foregrounds are briefly discussed in section 6. We also assume that the temperature and polarization fluctuations are Gaussian, thus neglecting any non-Gaussianity, either of primordial origin (which are anyway bounded to be small [25]), or coming from unresolved systematics (e.g., foreground residuals).

A final remark concerns the physical, late-time-Universe effects on the CMB fields due to the propagation of CMB photons from the last-scattering surface to the observer throughout the evolving large-scale structures. Weak gravitational lensing due to the gravitational potential of cosmological structures deflects CMB photons and modifies the observed statistics of CMB anisotropies with respect to the pattern arising at decoupling. In what follows, we will implicitly consider unensched CMB fields, i.e., we will ignore the effects of gravitational lensing for the sake of simplicity. The non-trivial modifications induced by the gravitational potential will be discussed later in section 4.

### 2.1. Statistics of CMB Temperature Field—Exact Likelihood in Harmonic Space

The CMB temperature field $T(\bar{x}, \hat{n}, \tau) = \bar{T}(\tau)[1 + \Theta(\bar{x}, \hat{n}, \tau)]$ observed<fnref n="3" /> in a given direction $\hat{n}$ is defined at every point $(\bar{x}, \tau)$ in space and time. The field has been decomposed in an isotropic background value $\bar{T}(\tau)$ and a small perturbation, the anisotropy field $\Theta(\bar{x}, \hat{n}, \tau) = (T - \bar{T})/\bar{T}$. Anisotropies are assumed to be the result of a Gaussian random process originated from quantum fluctuations in the early Universe. The observed temperature field generated by scalar fluctuations is a linear operator acting on three-dimensional perturbation fields:

$$ \Theta(\hat{n}) = \int d^3\bar{k} \xi(\bar{k}) \int_0^{\tau_0} d\tau e^{i(\bar{k} \cdot \hat{n})(\tau_0 - \tau)} S_T^{(s)}(k, \tau) \qquad (2) $$

where $\xi(\bar{k})$ is the primordial curvature perturbation, and the source function $S_T^{(s)}$ for scalar temperature fluctuations is a linear combination of the cosmological perturbation fields (see [26] for the explicit expression). A similar expression holds for temperature fluctuations generated by tensor perturbations [26].

<footnote n="3">³In this section, “observed” stands for what would be observed in the idealized situation considered here, i.e., negligible instrumental noise and absence of contaminants of any kind. A more proper term would be “realized,” but in this context we avoid to use it as we reserve it for simulated CMB sky signals.</footnote>
The aim of this manuscript is to review the basics of likelihood analysis in CMB experiments. This goal is motivated by the fact the we are at a crucial point in the history of observational CMB cosmology. The level of maturity and complexity reached by current CMB experiments boosted the theoretical efforts in finding smart solutions to the issue of identifying a suitable likelihood choice. A rich literature has been produced in this sense, although an exhaustive overview of the topic is not available, to the best of our knowledge. This review would fill the gap. Such a review could also serve as a good starting point for those who are approaching the field of CMB data analysis today or in the next future, and would be ideally contributing to the advances of CMB science in the next decades.

The structure of the manuscript is as follows. Section 2 is devoted to the statistics of the CMB fields. The approach to the topic is pedagogical, in a sense that we begin with a discussion in the single-field, temperature-only regime and introduce the basic statistical properties of CMB fluctuations. Then, we move to the more general case of correlated T, E, B fields. The discussion is carried over in the full-sky regime, with no distinction made between applications to large- and small-scales. We conclude section 2 with the introduction of the exact likelihood in full sky. Specific approximations to the exact likelihood are presented in section 3 (applications to the small-scale regime) and in section 5 (applications to the large-scale regime). In both sections, attention is devoted to complications due to partial-sky coverage and noise contamination. The inclusion of physical late-time Universe effects on the CMB photons in terms of gravitational lensing is detailed in section 4, whereas the important issues related to the presence of foreground emissions are described in section 6. The discussion of the various likelihood approaches in terms of computational cost (where applicable) and robustness with respect to the ability to provide unbiased estimations of cosmological parameters is detailed in section 7. Our conclusions are summarized in section 8. Some useful tools that will be mentioned throughout the main text are further discussed in Appendix A. In particular, in Appendix A, we review the basic notions of statistics needed to develop the formalism of CMB statistics. In Appendix B, we discuss the construction of power spectrum estimators, including pseudo-C_l, the “pure” formalism, and quadratic maximum likelihood (QML) estimator.

## 2. STATISTICS OF THE COSMIC MICROWAVE BACKGROUND

We now introduce some basic aspects of the statistics of the CMB. The basic object that we are interested in is the likelihood function L, i.e., the probability of the observed data d given a model, regarded as a function of the model itself. If the model is defined in terms of a vector of parameters θ, we thus have:

$$ \mathcal{L}(\theta) = p(d|\theta). \qquad (1) $$

The notation used throughout this review is presented in Appendix A, were we also recall some basic notion of probability and statistics.

We first derive the exact likelihood function for the CMB fields in harmonic space. The main statistical concepts are introduced in the limit of single field (section 2.1), i.e., temperature only, for the sake of simplicity. We then generalize these main findings in the case of joint temperature and polarization analysis (section 2.2). The exact likelihood in real space are derived in section 2.3.

We assume an ideal scenario of full-sky observations with infinite angular resolution and absence of noise and foreground contaminations. Obviously, this scenario is highly idealized. Nevertheless, it allows to easily derive the basic concepts of CMB statistics. Modifications to this picture arising from realistic observational issues (limited sky coverage, masked sky, experimental noise, finite angular resolution) are introduced in section 3.3. Foregrounds are briefly discussed in section 6. We also assume that the temperature and polarization fluctuations are Gaussian, thus neglecting any non-Gaussianity, either of primordial origin (which are anyway bounded to be small [25]), or coming from unresolved systematics (e.g., foreground residuals).

A final remark concerns the physical, late-time-Universe effects on the CMB fields due to the propagation of CMB photons from the last-scattering surface to the observer throughout the evolving large-scale structures. Weak gravitational lensing due to the gravitational potential of cosmological structures deflects CMB photons and modifies the observed statistics of CMB anisotropies with respect to the pattern arising at decoupling. In what follows, we will implicitly consider unensched CMB fields, i.e., we will ignore the effects of gravitational lensing for the sake of simplicity. The non-trivial modifications induced by the gravitational potential will be discussed later in section 4.

### 2.1. Statistics of CMB Temperature Field—Exact Likelihood in Harmonic Space

The CMB temperature field $T(\bar{x}, \hat{n}, \tau) = \bar{T}(\tau)[1 + \Theta(\bar{x}, \hat{n}, \tau)]$ observed<fnref n="3" /> in a given direction $\hat{n}$ is defined at every point $(\bar{x}, \tau)$ in space and time. The field has been decomposed in an isotropic background value $\bar{T}(\tau)$ and a small perturbation, the anisotropy field $\Theta(\bar{x}, \hat{n}, \tau) = (T - \bar{T})/\bar{T}$. Anisotropies are assumed to be the result of a Gaussian random process originated from quantum fluctuations in the early Universe. The observed temperature field generated by scalar fluctuations is a linear operator acting on three-dimensional perturbation fields:

$$ \Theta(\hat{n}) = \int d^3\bar{k} \xi(\bar{k}) \int_0^{\tau_0} d\tau e^{i(\bar{k} \cdot \hat{n})(\tau_0 - \tau)} S_T^{(s)}(k, \tau) \qquad (2) $$

where $\xi(\bar{k})$ is the primordial curvature perturbation, and the source function $S_T^{(s)}$ for scalar temperature fluctuations is a linear combination of the cosmological perturbation fields (see [26] for the explicit expression). A similar expression holds for temperature fluctuations generated by tensor perturbations [26].

<footnote n="3">³In this section, “observed” stands for what would be observed in the idealized situation considered here, i.e., negligible instrumental noise and absence of contaminants of any kind. A more proper term would be “realized,” but in this context we avoid to use it as we reserve it for simulated CMB sky signals.</footnote>
In Equation (2), we have suppressed the $\tau$ and $\bar{x}$ dependence in $\Theta$ as we are implicitly assuming that the temperature field is observed in a given position at a fixed point in time.

Assuming a given cosmological model, we cannot directly predict the particular realization of the temperature field. Instead, we shall infer statistical properties of the observed perturbation field. It is useful, to this purpose, to decompose the angular dependence of the temperature anisotropy field in spherical harmonics $Y_{\ell m}(\hat{p})$

$$ \Theta(\vec{x}, \hat{p}, \tau) = \sum_{\ell=1}^{\infty} \sum_{m=-\ell}^{\ell} a_{\ell m}(\vec{x}, \tau) Y_{\ell m}(\hat{p}) \quad (3) $$

where the harmonic $Y_{\ell m}$ corresponds to an angular scale $\theta \sim \pi/\ell$ with $(2\ell + 1)$ m-modes for each multipole $\ell$. Low multipoles (low-$\ell$) in the expansion correspond to large angular scales in the sky, whereas high multipoles (high-$\ell$) correspond to small scales. Since $\Theta$ is real, the decomposition coefficients $a_{\ell m}$ have to satisfy the reality condition

$$ a_{\ell m}^* = a_{\ell -m} \quad (4) $$

All the information about the $\bar{x}$ and $\tau$ dependence of $\Theta$ is now encoded in the $a_{\ell m}$'s.

We are interested in extracting information about the statistical properties of the $a_{\ell m}$'s from the observations. In the standard cosmological model, $a_{\ell m}$s follow a Gaussian distribution, with vanishing average ($\langle a_{\ell m} \rangle = 0$, since the $a_{\ell m}$ are expansion coefficients of the anisotropy field, whose mean vanishes), and covariance

$$ \langle a_{\ell m} a_{\ell' m'}^* \rangle = \delta_{\ell \ell'} \delta_{mm'} C_{\ell} \quad (5) $$

where the constraints imposed by the two Dirac delta functions follow from the $a_{\ell m}$ being independent random variables (diagonal covariance). Moreover, statistical isotropy ensures that the variance does not depend on $m$ (rotational invariance of $C_{\ell}$). The $C_{\ell}$'s are the angular power spectrum of the CMB temperature field. The power spectrum is related to the two-point correlation function of the field $C(\theta) = \langle \Theta(\hat{n}_1)\Theta(\hat{n}_2) \rangle$ observed at two directions $\hat{n}_1$ and $\hat{n}_2$ in the sky such that $\hat{n}_1 \cdot \hat{n}_2 = \cos \theta$:

$$ C(\theta) = \sum_{\ell} \frac{2\ell + 1}{4\pi} C_{\ell} P_{\ell}(\cos \theta), \quad (6) $$

where $P_{\ell}$ is the Legendre polynomial of order $\ell$.

If a random variable is Gaussian distributed, all the statistical properties are encoded in its mean and variance, which are the only momenta of the distribution we need to know. In fact, for a Gaussian distribution, odd momenta vanish and even momenta beyond the second can be recast as a function of the variance (Wick's theorem). Thus, the power spectrum $C_{\ell}$, or equivalently the two-point correlation function $C(\theta)$, completely characterizes the statistical properties of the anisotropy field.

Since the $a_{\ell m}$'s follow a Gaussian distribution with zero mean and variance $C_{\ell}$, we can readily write the probability density function $p(a_{\ell m}|C_{\ell})$ of the $a_{\ell m}$'s conditioned by the $C_{\ell}$'s:

$$ p(a_{\ell m} | C_{\ell}) = \frac{1}{\sqrt{2\pi} C_{\ell}} \exp \left( -\frac{|a_{\ell m}|^2}{2C_{\ell}} \right). \quad (7) $$

Given the observed temperature field and the corresponding $a_{\ell m}$'s, this expression already provides the likelihood function for the theoretical (model) $C_{\ell}$'s. However the information contained in the $a_{\ell m}$'s can be further compressed, as we shall see in the following.

Statistical isotropy of the $C_{\ell}$'s allows us to rewrite Equation (5) as:

$$ C_{\ell} = \frac{1}{2\ell + 1} \sum_{m=-\ell}^{\ell} \langle |a_{\ell m}|^2 \rangle \quad (8) $$

Some considerations are in order at this point. The average operation defined with the symbol $\langle...\rangle$ in Equations (5) and (8) is an ensemble average. As noted above, the CMB field is a realization of a random process and statistical information about the outcome of such a process should be obtained by averaging over all possible realizations. In practice, however, we can only observe a single realization of the CMB field. A way out is provided by the statistical omogeneity and isotropy of the CMB fluctuations, that in principle allows to substitute the ensemble average in Equation (5) with an average over different positions and directions. According to this *ergodic hypothesis*, different regions that are widely separated in the sky are statistically independent from each other and can be considered as different statistical realizations of the same stochastic process. Since we only have access to the CMB field observed at $\bar{x}_0$ and $\tau_0$, i.e., the CMB field here and now, what we are really left is the average over different directions, or equivalently over different values of $m$. In other words, for a given $\ell$, all the $a_{\ell m}$ are drawn from the same distribution, which can be therefore sampled by measuring all the $2\ell + 1$ coefficients. We are thus led to define an estimator of the observed power spectrum

$$ \hat{C}_{\ell} = \frac{1}{2\ell + 1} \sum_{m} |a_{\ell m}|^2, \quad (9) $$

with the property $\langle \hat{C}_{\ell} \rangle = C_{\ell}$. Note that in Equation (9) the ensemble average does not appear: we are forced to measure $C_{\ell}$ only with a limited number of values. This induces an intrinsic source of inaccuracy due to replacing the true variance $C_{\ell}$ with the observed power $\hat{C}_{\ell}$ (i.e., by replacing the ensemble average with the average over directions). This effect is known as *cosmic variance*:

$$ \left\langle \left( \frac{\hat{C}_{\ell} - C_{\ell}}{C_{\ell}} \right)^2 \right\rangle = -1 + \frac{1}{(2\ell + 1)^2 C_{\ell}^2} \sum_{m \text{ mm}'} \langle a_{\ell m} a_{\ell' m}^* a_{\ell m'} a_{\ell' m'}^* \rangle \\ = -1 + \frac{1}{(2\ell + 1)^2 C_{\ell}^2} \left[ \sum_m \langle a_{\ell m} a_{\ell' m}^* a_{\ell m} a_{\ell' m}^* \rangle \right] $$
In Equation (2), we have suppressed the $\tau$ and $\bar{x}$ dependence in $\Theta$ as we are implicitly assuming that the temperature field is observed in a given position at a fixed point in time.

Assuming a given cosmological model, we cannot directly predict the particular realization of the temperature field. Instead, we shall infer statistical properties of the observed perturbation field. It is useful, to this purpose, to decompose the angular dependence of the temperature anisotropy field in spherical harmonics $Y_{\ell m}(\hat{p})$

$$ \Theta(\vec{x}, \hat{p}, \tau) = \sum_{\ell=1}^{\infty} \sum_{m=-\ell}^{\ell} a_{\ell m}(\vec{x}, \tau) Y_{\ell m}(\hat{p}) \quad (3) $$

where the harmonic $Y_{\ell m}$ corresponds to an angular scale $\theta \sim \pi/\ell$ with $(2\ell + 1)$ m-modes for each multipole $\ell$. Low multipoles (low-$\ell$) in the expansion correspond to large angular scales in the sky, whereas high multipoles (high-$\ell$) correspond to small scales. Since $\Theta$ is real, the decomposition coefficients $a_{\ell m}$ have to satisfy the reality condition

$$ a_{\ell m}^* = a_{\ell -m} \quad (4) $$

All the information about the $\bar{x}$ and $\tau$ dependence of $\Theta$ is now encoded in the $a_{\ell m}$'s.

We are interested in extracting information about the statistical properties of the $a_{\ell m}$'s from the observations. In the standard cosmological model, $a_{\ell m}$s follow a Gaussian distribution, with vanishing average $\langle a_{\ell m} \rangle = 0$, since the $a_{\ell m}$ are expansion coefficients of the anisotropy field, whose mean vanishes), and covariance

$$ \langle a_{\ell m} a_{\ell' m'}^* \rangle = \delta_{\ell \ell'} \delta_{mm'} C_{\ell} \quad (5) $$

where the constraints imposed by the two Dirac delta functions follow from the $a_{\ell m}$ being independent random variables (diagonal covariance). Moreover, statistical isotropy ensures that the variance does not depend on $m$ (rotational invariance of $C_{\ell}$). The $C_{\ell}$'s are the angular power spectrum of the CMB temperature field. The power spectrum is related to the two-point correlation function of the field $C(\theta) = \langle \Theta(\hat{n}_1)\Theta(\hat{n}_2) \rangle$ observed at two directions $\hat{n}_1$ and $\hat{n}_2$ in the sky such that $\hat{n}_1 \cdot \hat{n}_2 = \cos \theta$:

$$ C(\theta) = \sum_{\ell} \frac{2\ell + 1}{4\pi} C_{\ell} P_{\ell}(\cos \theta), \quad (6) $$

where $P_{\ell}$ is the Legendre polynomial of order $\ell$.

If a random variable is Gaussian distributed, all the statistical properties are encoded in its mean and variance, which are the only momenta of the distribution we need to know. In fact, for a Gaussian distribution, odd momenta vanish and even momenta beyond the second can be recast as a function of the variance (Wick's theorem). Thus, the power spectrum $C_{\ell}$, or equivalently the two-point correlation function $C(\theta)$, completely characterizes the statistical properties of the anisotropy field.

Since the $a_{\ell m}$'s follow a Gaussian distribution with zero mean and variance $C_{\ell}$, we can readily write the probability density function $p(a_{\ell m}|C_{\ell})$ of the $a_{\ell m}$'s conditioned by the $C_{\ell}$'s:

$$ p(a_{\ell m} | C_{\ell}) = \frac{1}{\sqrt{2\pi} C_{\ell}} \exp \left( -\frac{|a_{\ell m}|^2}{2C_{\ell}} \right). \quad (7) $$

Given the observed temperature field and the corresponding $a_{\ell m}$'s, this expression already provides the likelihood function for the theoretical (model) $C_{\ell}$'s. However the information contained in the $a_{\ell m}$'s can be further compressed, as we shall see in the following.

Statistical isotropy of the $C_{\ell}$'s allows us to rewrite Equation (5) as:

$$ C_{\ell} = \frac{1}{2\ell + 1} \sum_{m=-\ell}^{\ell} \langle |a_{\ell m}|^2 \rangle \quad (8) $$

Some considerations are in order at this point. The average operation defined with the symbol $\langle...\rangle$ in Equations (5) and (8) is an ensemble average. As noted above, the CMB field is a realization of a random process and statistical information about the outcome of such a process should be obtained by averaging over all possible realizations. In practice, however, we can only observe a single realization of the CMB field. A way out is provided by the statistical omogeneity and isotropy of the CMB fluctuations, that in principle allows to substitute the ensemble average in Equation (5) with an average over different positions and directions. According to this *ergodic hypothesis*, different regions that are widely separated in the sky are statistically independent from each other and can be considered as different statistical realizations of the same stochastic process. Since we only have access to the CMB field observed at $\bar{x}_0$ and $\tau_0$, i.e., the CMB field here and now, what we are really left is the average over different directions, or equivalently over different values of $m$. In other words, for a given $\ell$, all the $a_{\ell m}$ are drawn from the same distribution, which can be therefore sampled by measuring all the $2\ell + 1$ coefficients. We are thus led to define an estimator of the observed power spectrum

$$ \hat{C}_{\ell} = \frac{1}{2\ell + 1} \sum_{m} |a_{\ell m}|^2, \quad (9) $$

with the property $\langle \hat{C}_{\ell} \rangle = C_{\ell}$. Note that in Equation (9) the ensemble average does not appear: we are forced to measure $C_{\ell}$ only with a limited number of values. This induces an intrinsic source of inaccuracy due to replacing the true variance $C_{\ell}$ with the observed power $\hat{C}_{\ell}$ (i.e., by replacing the ensemble average with the average over directions). This effect is known as *cosmic variance*:

$$ \begin{aligned} \left\langle \left( \frac{\hat{C}_{\ell} - C_{\ell}}{C_{\ell}} \right)^2 \right\rangle &= -1 + \frac{1}{(2\ell + 1)^2 C_{\ell}^2} \sum_{m \text{ mm}'} \langle a_{\ell m} a_{\ell' m}^* a_{\ell m'} a_{\ell' m'}^* \rangle \\ &= -1 + \frac{1}{(2\ell + 1)^2 C_{\ell}^2} \left[ \sum_m \langle a_{\ell m} a_{\ell' m}^* a_{\ell m} a_{\ell' m}^* \rangle \right] \end{aligned} $$
$$ 
\begin{aligned}
& + \sum_{m, m' \neq m} \langle a_{\ell m} a_{\ell' m}^* a_{\ell m'} a_{\ell' m'}^* \rangle \\
& = -1 + \frac{1}{(2\ell + 1)^2 C_{\ell}^2} (3C_{\ell}^2(2\ell + 1) \\
& \qquad + 2\ell C_{\ell}^2(2\ell + 1)) \\
& = \frac{2}{2\ell + 1} \qquad (10)
\end{aligned}
$$

where the third equality follows from Wick's theorem.

Cosmic variance is an irreducible source of uncertainty in cosmological measurements of the CMB power spectrum, and one of the major sources of uncertainties especially at the largest scales (low-$\ell$), where we have only a limited number of coefficients $a_{\ell m}$ to average over with respect to the small-scale (high-$\ell$) regime. Equation (10) is valid provided full-sky observations. However, in real data analysis, even if we are able to observe the full sky (e.g., with space missions), we are nevertheless forced to mask a certain fraction of the sky, e.g., to avoid foreground contamination. An approximate estimate of the increase is given by a factor $1/f_{\text{sky}}$, where $f_{\text{sky}}$ but see e.g., [27] for a careful counting of the degrees of freedom available in cut-sky regimes. Current experiments like the Planck satellite are ideally cosmic-variance-limited up to very high multipoles, i.e., $\ell \sim 1500$.

To derive the distribution of the observed $C_\ell$'s, we note that the sum of $\nu = (2\ell + 1)$ standard Gaussian variables follows a $\chi^2$ distribution with $\nu$ degrees of freedom. If we define $\hat{Y}_\ell = \sum_m (|a_{\ell m}|/\sqrt{C_\ell}|^2)$, this new variable has a $\chi^2$ distribution:

$$ p(\hat{Y}_{\ell} | C_{\ell}) = \frac{\hat{Y}_{\ell}^{\nu/2-1}}{\Gamma(\nu/2) 2^{\nu/2}} \exp \left( -\frac{\hat{Y}_{\ell}}{2} \right) \quad (11) $$

The estimator (hereafter observed) $\hat{C}_\ell$ is a multiple of $\hat{Y}_\ell$: $\hat{C}_\ell = C_\ell \hat{Y}_\ell / (2\ell + 1)$, and multiples of $\chi^2$-distributed variables follow a Gamma distribution:

$$ p(\hat{C}_{\ell} | C_{\ell}) \propto C_{\ell}^{-1} \left( \frac{\hat{C}_{\ell}}{C_{\ell}} \right)^{\nu/2-1} \exp \left( -\frac{\nu}{2} \frac{\hat{C}_{\ell}}{C_{\ell}} \right) \quad (12) $$

The previous expression is the probability of the observed power spectrum given the fiducial power, and for fixed data it can be still regarded as a likelihood $\mathcal{L}(C_\ell)$, in which the role of the data is not played by the $a_{\ell m}$'s as in Equation (7), but by the $\hat{C}_\ell$. The mean and variance of the distribution of the $\hat{C}_\ell$'s are $E[\hat{C}_\ell] = C_\ell$ and $\text{Var}[\hat{C}_\ell] = 2C_\ell^2/\nu$. The maximum of the distribution is in $(\nu - 2)/\nu C_\ell$, that does not coincide with the mean of the distribution. As such, the distribution of observed $\hat{C}_\ell$ is skewed. However, in the limit $\nu \to \infty$, the distribution in Equation (12) tends to a Gaussian distribution with same mean and variance, according to the central limit theorem. Note that the variance of the distribution is exactly the cosmic variance introduced in Equation (10). This further stresses the meaning of the cosmic variance as an irreducible source of uncertainty due to the limitation of having access to a single realization of the Universe (i.e., the limitation due to estimating the true power spectrum $C_\ell$ with the observed power spectrum $\hat{C}_\ell$).

## 2.2. Statistics of Joint CMB Temperature and Polarization Fields—Exact Likelihood in Harmonic Space

The above treatment has to be generalized in the case of the joint analysis of temperature and polarization fields T, E, B. In analogy to the temperature case, we can define two sets of spherical harmonics coefficients for E and B:

$$ a_{\ell m}^{E} \equiv -\frac{1}{2} (+2a_{\ell m} + -2 a_{\ell m}) \quad (13a) $$

$$ a_{\ell m}^{B} \equiv \frac{i}{2} (+2a_{\ell m} - -2 a_{\ell m}) \quad (13b) $$

where $\pm 2a_{\ell m}$ are the expansion coefficients of the combinations of Stokes parameters describing the polarization state of the CMB signal—($Q \pm iU$)—in spin-2 spherical harmonics $\pm 2Y_{\ell m}$ (see e.g., [28] for a derivation of the formalism).

The variable $\mathbf{X}_a = (a_{\ell m}^T, a_{\ell m}^E, a_{\ell m}^B)$ is distributed according to a Gaussian multivariate distribution with covariance matrix

$$ \operatorname{cov}[\mathbf{X}_a, \mathbf{X}_a] \equiv \langle \mathbf{X}_a \mathbf{X}_a^\dagger \rangle - (|\mathbf{X}_a|)^2 \quad (14a) $$

$$ \equiv \mathbf{V}_{\ell} = \begin{pmatrix} C_{\ell}^{TT} & C_{\ell}^{TE} & 0 \\ C_{\ell}^{TE} & C_{\ell}^{EE} & 0 \\ 0 & 0 & C_{\ell}^{BB} \end{pmatrix}, \quad (14b) $$

where it is explicitly seen that the temperature and the E-polarization fields are correlated, whereas the parity-even fields (T and E) are uncorrelated with the parity-odd field B (although this is strictly true only in the standard cosmological model when parity violation processes are forbidden in the early Universe).

In analogy to Equation (9), the estimators for the observed power spectra are given by the following matrix:

$$ \mathbf{S}_{\ell} = \frac{1}{2\ell + 1} \sum_{m} \mathbf{X}_{a} \mathbf{X}_{a}^{\dagger} = \begin{pmatrix} \hat{C}_{\ell}^{TT} & \hat{C}_{\ell}^{TE} & \hat{C}_{\ell}^{TB} \\ \hat{C}_{\ell}^{TE} & \hat{C}_{\ell}^{EE} & \hat{C}_{\ell}^{EB} \\ \hat{C}_{\ell}^{TB} & \hat{C}_{\ell}^{EB} & \hat{C}_{\ell}^{BB} \end{pmatrix}, \quad (15) $$

where the observed cross-correlations TB and EB may be non-vanishing as well.

The probability of $\mathbf{X}_a$ at each $\ell$ can therefore be expressed as:

$$ -2 \ln[p(\mathbf{X}_a | \mathbf{V}_\ell)] = \mathbf{X}_a^\dagger \mathbf{V}_\ell^{-1} \mathbf{X}_a + \ln \det[2\pi \mathbf{V}_\ell] \quad (16a) $$

$$ = (2\ell + 1)\text{trace}[\mathbf{S}_\ell \mathbf{V}_\ell^{-1}] + \ln \det[\mathbf{V}_\ell] + \text{const.} \quad (16b) $$

Note that $\mathbf{S}_\ell$ represents sufficient statistics for this likelihood function: in the full-sky regime, Equation (16a) only depends on the data through $\mathbf{S}_\ell$ and therefore information on the CMB sky can be losslessly compressed to a set of power spectrum estimators $\mathbf{S}_\ell = \mathbf{S}_\ell(\hat{C}_\ell^{XY})$, $X, Y = T, E, B$.

The probability of $\mathbf{S}_\ell$ given $\mathbf{V}_\ell = \mathbf{V}_\ell(C_\ell^{XY})$ is obtained by properly normalizing (Equation 16a). In the previous section, we have seen that the single-field $\hat{C}_\ell$ is Gamma-distributed. It is easy to understand that the full set of observed power spectra $\mathbf{S}_\ell$ has
$$ 
\begin{aligned}
& + \sum_{m, m' \neq m} \langle a_{\ell m} a_{\ell' m}^* a_{\ell m'} a_{\ell' m'}^* \rangle \\
& = -1 + \frac{1}{(2\ell + 1)^2 C_{\ell}^2} (3C_{\ell}^2(2\ell + 1) \\
& \qquad + 2\ell C_{\ell}^2(2\ell + 1)) \\
& = \frac{2}{2\ell + 1} \qquad (10)
\end{aligned}
$$

where the third equality follows from Wick's theorem.

Cosmic variance is an irreducible source of uncertainty in cosmological measurements of the CMB power spectrum, and one of the major sources of uncertainties especially at the largest scales (low-$\ell$), where we have only a limited number of coefficients $a_{\ell m}$ to average over with respect to the small-scale (high-$\ell$) regime. Equation (10) is valid provided full-sky observations. However, in real data analysis, even if we are able to observe the full sky (e.g., with space missions), we are nevertheless forced to mask a certain fraction of the sky, e.g., to avoid foreground contamination. An approximate estimate of the increase is given by a factor $1/f_{sky}$, where $f_{sky}$ but see e.g., [27] for a careful counting of the degrees of freedom available in cut-sky regimes. Current experiments like the Planck satellite are ideally cosmic-variance-limited up to very high multipoles, i.e., $\ell \sim 1500$.

To derive the distribution of the observed $C_\ell$'s, we note that the sum of $\nu = (2\ell + 1)$ standard Gaussian variables follows a $\chi^2$ distribution with $\nu$ degrees of freedom. If we define $\hat{Y}_\ell = \sum_m (|a_{\ell m}|/\sqrt{C_\ell}|^2)$, this new variable has a $\chi^2$ distribution:

$$ p(\hat{Y}_{\ell} | C_{\ell}) = \frac{\hat{Y}_{\ell}^{\nu/2-1}}{\Gamma(\nu/2) 2^{\nu/2}} \exp \left( -\frac{\hat{Y}_{\ell}}{2} \right) \quad (11) $$

The estimator (hereafter observed) $\hat{C}_\ell$ is a multiple of $\hat{Y}_\ell$: $\hat{C}_\ell = C_\ell \hat{Y}_\ell / (2\ell + 1)$, and multiples of $\chi^2$-distributed variables follow a Gamma distribution:

$$ p(\hat{C}_{\ell} | C_{\ell}) \propto C_{\ell}^{-1} \left( \frac{\hat{C}_{\ell}}{C_{\ell}} \right)^{\nu/2-1} \exp \left( -\frac{\nu}{2} \frac{\hat{C}_{\ell}}{C_{\ell}} \right) \quad (12) $$

The previous expression is the probability of the observed power spectrum given the fiducial power, and for fixed data it can be still regarded as a likelihood $\mathcal{L}(C_\ell)$, in which the role of the data is not played by the $a_{\ell m}$'s as in Equation (7), but by the $\hat{C}_\ell$. The mean and variance of the distribution of the $\hat{C}_\ell$'s are $E[\hat{C}_\ell] = C_\ell$ and $\text{Var}[\hat{C}_\ell] = 2C_\ell^2/\nu$. The maximum of the distribution is in $(\nu - 2)/\nu C_\ell$, that does not coincide with the mean of the distribution. As such, the distribution of observed $\hat{C}_\ell$ is skewed. However, in the limit $\nu \to \infty$, the distribution in Equation (12) tends to a Gaussian distribution with same mean and variance, according to the central limit theorem. Note that the variance of the distribution is exactly the cosmic variance introduced in Equation (10). This further stresses the meaning of the cosmic variance as an irreducible source of uncertainty due to the limitation of having access to a single realization of the Universe (i.e., the limitation due to estimating the true power spectrum $C_\ell$ with the observed power spectrum $\hat{C}_\ell$).

## 2.2. Statistics of Joint CMB Temperature and Polarization Fields—Exact Likelihood in Harmonic Space

The above treatment has to be generalized in the case of the joint analysis of temperature and polarization fields T, E, B. In analogy to the temperature case, we can define two sets of spherical harmonics coefficients for E and B:

$$ a_{\ell m}^{E} \equiv -\frac{1}{2} (+2a_{\ell m} + -2 a_{\ell m}) \quad (13a) $$

$$ a_{\ell m}^{B} \equiv \frac{i}{2} (+2a_{\ell m} - 2 a_{\ell m}) \quad (13b) $$

where $\pm 2a_{\ell m}$ are the expansion coefficients of the combinations of Stokes parameters describing the polarization state of the CMB signal—($Q \pm iU$)—in spin-2 spherical harmonics $\pm 2Y_{\ell m}$ (see e.g., [28] for a derivation of the formalism).

The variable $\mathbf{X}_a = (a_{\ell m}^T, a_{\ell m}^E, a_{\ell m}^B)$ is distributed according to a Gaussian multivariate distribution with covariance matrix

$$ \text{cov}[\mathbf{X}_a, \mathbf{X}_a] \equiv \langle \mathbf{X}_a \mathbf{X}_a^\dagger \rangle - (|\langle \mathbf{X}_a \rangle|^2) \quad (14a) $$

$$ \equiv \mathbf{V}_\ell = \begin{pmatrix} C_\ell^{TT} & C_\ell^{TE} & 0 \\ C_\ell^{TE} & C_\ell^{EE} & 0 \\ 0 & 0 & C_\ell^{BB} \end{pmatrix}, \quad (14b) $$

where it is explicitly seen that the temperature and the E-polarization fields are correlated, whereas the parity-even fields (T and E) are uncorrelated with the parity-odd field B (although this is strictly true only in the standard cosmological model when parity violation processes are forbidden in the early Universe).

In analogy to Equation (9), the estimators for the observed power spectra are given by the following matrix:

$$ \mathbf{S}_{\ell} = \frac{1}{2\ell + 1} \sum_{m} \mathbf{X}_{a} \mathbf{X}_{a}^{\dagger} = \begin{pmatrix} \hat{C}_{\ell}^{TT} & \hat{C}_{\ell}^{TE} & \hat{C}_{\ell}^{TB} \\ \hat{C}_{\ell}^{TE} & \hat{C}_{\ell}^{EE} & \hat{C}_{\ell}^{EB} \\ \hat{C}_{\ell}^{TB} & \hat{C}_{\ell}^{EB} & \hat{C}_{\ell}^{BB} \end{pmatrix}, \quad (15) $$

where the observed cross-correlations TB and EB may be non-vanishing as well.

The probability of $\mathbf{X}_a$ at each $\ell$ can therefore be expressed as:

$$ -2 \ln[p(\mathbf{X}_a | \mathbf{V}_\ell)] = \mathbf{X}_a^\dagger \mathbf{V}_\ell^{-1} \mathbf{X}_a + \ln \det[2\pi \mathbf{V}_\ell] \quad (16a) $$

$$ = (2\ell + 1)\text{trace}[\mathbf{S}_\ell \mathbf{V}_\ell^{-1}] + \ln \det[\mathbf{V}_\ell] + \text{const.} \quad (16b) $$

Note that $\mathbf{S}_\ell$ represents sufficient statistics for this likelihood function: in the full-sky regime, Equation (16a) only depends on the data through $\mathbf{S}_\ell$ and therefore information on the CMB sky can be losslessly compressed to a set of power spectrum estimators $\mathbf{S}_\ell = \mathbf{S}_\ell(\hat{C}_\ell^{XY})$, $X, Y = T, E, B$.

The probability of $\mathbf{S}_\ell$ given $\mathbf{V}_\ell = \mathbf{V}_\ell(C_\ell^{XY})$ is obtained by properly normalizing (Equation 16a). In the previous section, we have seen that the single-field $\hat{C}_\ell$ is Gamma-distributed. It is easy to understand that the full set of observed power spectra $\mathbf{S}_\ell$ has
a Wishart distribution, i.e., a multi-dimensional generalization of the Gamma distribution, with $\nu = (2\ell + 1)$ degrees of freedom in $p = 3$ dimensions:

$$p(\mathbf{S}_{\ell}|\mathbf{W}_{\ell}) = \mathcal{L}(\mathbf{W}_{\ell}) = \frac{|\mathbf{S}_{\ell}|^{(\nu-p-1)/2} \exp[-\mathrm{trace}(\mathbf{W}_{\ell}^{-1}\mathbf{S}_{\ell}/2)]}{2^{\nu/2}|\mathbf{W}_{\ell}|^{\nu/2}\Gamma_p(\nu/2)} \quad (17)$$

where $\mathbf{W}_{\ell} = \mathbf{V}_{\ell}/\nu$. For given $\hat{C}_{\ell}^{XY}$'s, Equation (17) represents the exact expression of the likelihood function of the $\hat{C}_{\ell}^{XY}$.

Since $\mathbf{V}_{\ell}$ is separable in the two blocks TE and B, we can simplify the problem and consider two separate Wishart distributions for the block TE and for the block B:

$$\mathcal{L}(\mathbf{W}_{\ell}) = \mathcal{L}(\mathbf{W}_{\ell}^{\mathrm{TE}})\mathcal{L}(\mathbf{W}_{\ell}^{\mathrm{B}}) \quad (18)$$

The latter is further simplified since it reduces to the one-dimensional Gamma distribution, as described in details in the previous section. The distribution for the TE block can be fully expanded as:

$$\mathcal{L}(\mathbf{W}_{\ell}^{\mathrm{TE}}) \propto \frac{(C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - (C_{\ell}^{\mathrm{TE}})^2)^{(\nu-3)/2}}{(C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - (C_{\ell}^{\mathrm{TE}})^2)^{\nu/2}} \times \exp \left\{ -\frac{\nu}{2} \left[ \frac{C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} + C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - 2C_{\ell}^{\mathrm{TE}} C_{\ell}^{\mathrm{TE}}}{C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - (C_{\ell}^{\mathrm{TE}})^2} \right] \right\} \quad (19a)$$

The marginal distribution of each individual diagonal element of $\mathbf{S}_{\ell}^{\mathrm{TE}}$ can be obtained by integrating $p(\mathbf{S}_{\ell}|\mathbf{W}_{\ell})^{\mathrm{TE}}$ over $\hat{C}_{\ell}^{\mathrm{TE}}$ and the other diagonal element, and it is again a Gamma distribution as we expect it to be, in analogy to discussion in the previous section. However, the marginal distribution of the off-diagonal terms $\hat{C}_{\ell}^{\mathrm{TE}}$ is not a Gamma distribution, and it is interesting to note that it depends on $C_{\ell}^{\mathrm{TT}}$ and $C_{\ell}^{\mathrm{EE}}$ in addition to $C_{\ell}^{\mathrm{TE}}$ (see [29] for a detailed calculation).

In the limit $\nu \to \infty$, the Wishart distribution of $\hat{\mathbf{X}}_C = (\hat{C}_{\ell}^{\mathrm{TT}}, \hat{C}_{\ell}^{\mathrm{TE}}, \hat{C}_{\ell}^{\mathrm{EE}})$ tends to a multivariate Gaussian distribution with covariance matrix:

$$\begin{align} \mathrm{cov}[\mathbf{X}_C, \mathbf{X}_C] &\equiv \langle \mathbf{X}_C \mathbf{X}_C^\dagger \rangle - \langle \mathbf{X}_C \rangle^2 \equiv \mathbf{C} \tag{20a} \\ &= \frac{1}{2\ell+1} \begin{pmatrix} 2(C_{\ell}^{\mathrm{TT}})^2 & 2C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{TE}} & 2(C_{\ell}^{\mathrm{TE}})^2 \\ 2C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{TE}} & C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} + (C_{\ell}^{\mathrm{TE}})^2 & 2C_{\ell}^{\mathrm{TE}} C_{\ell}^{\mathrm{EE}} \\ 2(C_{\ell}^{\mathrm{TE}})^2 & 2C_{\ell}^{\mathrm{TE}} C_{\ell}^{\mathrm{EE}} & 2(C_{\ell}^{\mathrm{EE}})^2 \end{pmatrix}, \tag{20b} \end{align}$$

The variance of $\hat{C}_{\ell}^{\mathrm{TT}}$ and $\hat{C}_{\ell}^{\mathrm{EE}}$ is the same of the single-field limit, whereas the variance of the cross-correlation $\hat{C}_{\ell}^{\mathrm{TE}}$ reflects the different marginalized distribution of $\hat{C}_{\ell}^{\mathrm{TE}}$ itself.

## 2.3. Statistics of Joint CMB Temperature and Polarization Fields—Exact Likelihood in Real Space

The discussion in sections 2.1 and 2.2 refers to the CMB statistics in harmonic space, i.e., the space in which the CMB fields are

expanded in spherical harmonics and the physical information are encoded in the expansion coefficients $a_{\ell m}$. In this subsection, we will review the basics of CMB statistics in real space.

The starting point are the observed CMB maps of the three Stokes parameters $T$, $Q$, and $U$. These maps can be discretized into $N$ pixels and arranged in $N$-dimensional vectors $\mathbf{T}$, $\mathbf{Q}$, and $\mathbf{U}$. As discussed in the previous section, the statistical properties of these objects are fully encoded in the auto- and cross- power spectra $C_{\ell}^{XY}$, with $X, Y = \{T, E, B\}$ for temperature, E-mode and B-mode polarization.

The exact likelihood function in real space (also called the pixel-based likelihood) is defined as

$$\mathcal{L}(C_{\ell}) = p(\mathbf{m}|C_{\ell}) = \frac{1}{2\pi|M|^{1/2}} \exp\left(-\frac{1}{2}\mathbf{m}^{\mathrm{T}}M^{-1}\mathbf{m}\right), \quad (21)$$

where $\mathbf{m}$ is the vector with $3N$ elements built from the justaxposition of $\mathbf{T}$, $\mathbf{Q}$, and $\mathbf{U}$, and $M$ is the total covariance matrix. The matrix $M$ depends only on the angle between two directions in the sky $\hat{n}_{i,j}$

$$\mathbf{M}(\hat{n}_i \cdot \hat{n}_j) = \begin{pmatrix} \langle T_i T_j \rangle & \langle T_i Q_j \rangle & \langle T_i U_j \rangle \\ \langle T_i Q_j \rangle & \langle Q_i Q_j \rangle & \langle Q_i U_j \rangle \\ \langle T_i U_j \rangle & \langle Q_i U_j \rangle & \langle U_i U_j \rangle \end{pmatrix}. \quad (22)$$

The $(3 \times 3)$ entries in Equation (22) for any given pair of pixels $ij$ depend on the Legendre polynomial $P_{\ell}$ and the fiducial power spectra. As a straightforward example, the entry $\langle T_i T_j \rangle$ is the expression

$$\langle T_i T_j \rangle = \sum_{\ell} \frac{2\ell + 1}{4\pi} P_{\ell}(\hat{r}_i \cdot \hat{r}_j) C_{\ell}^{\mathrm{TT}}. \quad (23)$$

where $P_{\ell}(\hat{r}_i \cdot \hat{r}_j) = \frac{4\pi}{2\ell+1} \sum_m Y_{\ell m}(\hat{r}_i) Y_{\ell m}^*(\hat{r}_j)$. A detailed description of the full procedure to obtain the covariance matrix, together with the expressions of the $(3 \times 3)$ entries, can be found in appendix A of Tegmark and de Oliveira-Costa [30].

It should be noted that the pixel-based likelihood in Equation (21) is exact even in the case of partial sky coverage; this not the case for the likelihood in harmonic space in Equation (17). Note however that it is still possible to derive an exact form for the harmonic-space likelihood even for partial sky coverage [31]. The pixel-based approach ensures mathematical rigor in the evaluation of the likelihood function. Nevertheless, it is highly expensive from a computational point of view. Indeed, the number of pixels needed to retain the information in the first $\ell_{\mathrm{max}}$ multipoles of the power spectrum scales as $\ell_{\mathrm{max}}^2$ and therefore the Cholesky decomposition required to evaluate the inverse of the covariance matrix in Equation (21) scales roughly as $\ell_{\mathrm{max}}^6$, where $\ell_{\mathrm{max}}$ is the highest multipole retained in the analysis. The computational cost is therefore driven by the evaluation of inverse matrix and determinants and becomes prohibitive for $\ell_{\mathrm{max}}$ larger than few hundreds. For this reason, this exact approach is feasible only to study large angular scales, where the information is contained in a relatively small number of multipoles.
a Wishart distribution, i.e., a multi-dimensional generalization of the Gamma distribution, with $\nu = (2\ell + 1)$ degrees of freedom in $p = 3$ dimensions:

$$p(\mathbf{S}_{\ell}|\mathbf{W}_{\ell}) = \mathcal{L}(\mathbf{W}_{\ell}) = \frac{|\mathbf{S}_{\ell}|^{(\nu-p-1)/2} \exp[-\operatorname{trace}(\mathbf{W}_{\ell}^{-1}\mathbf{S}_{\ell}/2)]}{2^{\nu/2}|\mathbf{W}_{\ell}|^{\nu/2}\Gamma_p(\nu/2)} \quad (17)$$

where $\mathbf{W}_{\ell} = \mathbf{V}_{\ell}/\nu$. For given $\hat{C}_{\ell}^{XY}$'s, Equation (17) represents the exact expression of the likelihood function of the $\hat{C}_{\ell}^{XY}$.

Since $\mathbf{V}_{\ell}$ is separable in the two blocks TE and B, we can simplify the problem and consider two separate Wishart distributions for the block TE and for the block B:

$$\mathcal{L}(\mathbf{W}_{\ell}) = \mathcal{L}(\mathbf{W}_{\ell}^{\mathrm{TE}})\mathcal{L}(\mathbf{W}_{\ell}^{\mathrm{B}}) \quad (18)$$

The latter is further simplified since it reduces to the one-dimensional Gamma distribution, as described in details in the previous section. The distribution for the TE block can be fully expanded as:

$$\mathcal{L}(\mathbf{W}_{\ell}^{\mathrm{TE}}) \propto \frac{(C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - (C_{\ell}^{\mathrm{TE}})^2)^{(\nu-3)/2}}{(C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - (C_{\ell}^{\mathrm{TE}})^2)^{\nu/2}} \times \exp \left\{ -\frac{\nu}{2} \left[ \frac{C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} + C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - 2C_{\ell}^{\mathrm{TE}} C_{\ell}^{\mathrm{TE}}}{C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} - (C_{\ell}^{\mathrm{TE}})^2} \right] \right\} \quad (19a)$$

The marginal distribution of each individual diagonal element of $\mathbf{S}_{\ell}^{\mathrm{TE}}$ can be obtained by integrating $p(\mathbf{S}_{\ell}|\mathbf{W}_{\ell})^{\mathrm{TE}}$ over $\hat{C}_{\ell}^{\mathrm{TE}}$ and the other diagonal element, and it is again a Gamma distribution as we expect it to be, in analogy to discussion in the previous section. However, the marginal distribution of the off-diagonal terms $\hat{C}_{\ell}^{\mathrm{TE}}$ is not a Gamma distribution, and it is interesting to note that it depends on $C_{\ell}^{\mathrm{TT}}$ and $C_{\ell}^{\mathrm{EE}}$ in addition to $C_{\ell}^{\mathrm{TE}}$ (see [29] for a detailed calculation).

In the limit $\nu \to \infty$, the Wishart distribution of $\hat{\mathbf{X}}_C = (\hat{C}_{\ell}^{\mathrm{TT}}, \hat{C}_{\ell}^{\mathrm{TE}}, \hat{C}_{\ell}^{\mathrm{EE}})$ tends to a multivariate Gaussian distribution with covariance matrix:

$$\begin{align} \operatorname{cov}[\mathbf{X}_C, \mathbf{X}_C] &\equiv \langle \mathbf{X}_C \mathbf{X}_C^\dagger \rangle - \langle \mathbf{X}_C \rangle^2 &\equiv \mathbf{C} \tag{20a} \\ &= \frac{1}{2\ell+1} \begin{pmatrix} 2(C_{\ell}^{\mathrm{TT}})^2 & 2C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{TE}} & 2(C_{\ell}^{\mathrm{TE}})^2 \\ 2C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{TE}} & C_{\ell}^{\mathrm{TT}} C_{\ell}^{\mathrm{EE}} + (C_{\ell}^{\mathrm{TE}})^2 & 2C_{\ell}^{\mathrm{TE}} C_{\ell}^{\mathrm{EE}} \\ 2(C_{\ell}^{\mathrm{TE}})^2 & 2C_{\ell}^{\mathrm{TE}} C_{\ell}^{\mathrm{EE}} & 2(C_{\ell}^{\mathrm{EE}})^2 \end{pmatrix}, \tag{20b} \end{align}$$

The variance of $\hat{C}_{\ell}^{\mathrm{TT}}$ and $\hat{C}_{\ell}^{\mathrm{EE}}$ is the same of the single-field limit, whereas the variance of the cross-correlation $\hat{C}_{\ell}^{\mathrm{TE}}$ reflects the different marginalized distribution of $\hat{C}_{\ell}^{\mathrm{TE}}$ itself.

## 2.3. Statistics of Joint CMB Temperature and Polarization Fields—Exact Likelihood in Real Space

The discussion in sections 2.1 and 2.2 refers to the CMB statistics in harmonic space, i.e., the space in which the CMB fields are

expanded in spherical harmonics and the physical information are encoded in the expansion coefficients $a_{\ell m}$. In this subsection, we will review the basics of CMB statistics in real space.

The starting point are the observed CMB maps of the three Stokes parameters $T$, $Q$, and $U$. These maps can be discretized into $N$ pixels and arranged in $N$-dimensional vectors $\mathbf{T}$, $\mathbf{Q}$, and $\mathbf{U}$. As discussed in the previous section, the statistical properties of these objects are fully encoded in the auto- and cross- power spectra $C_{\ell}^{XY}$, with $X, Y = \{T, E, B\}$ for temperature, E-mode and B-mode polarization.

The exact likelihood function in real space (also called the pixel-based likelihood) is defined as

$$\mathcal{L}(C_{\ell}) = p(\mathbf{m}|C_{\ell}) = \frac{1}{2\pi|M|^{1/2}} \exp\left(-\frac{1}{2}\mathbf{m}^{\mathrm{T}}M^{-1}\mathbf{m}\right), \quad (21)$$

where $\mathbf{m}$ is the vector with $3N$ elements built from the justaxposition of $\mathbf{T}$, $\mathbf{Q}$, and $\mathbf{U}$, and $M$ is the total covariance matrix. The matrix $M$ depends only on the angle between two directions in the sky $\hat{n}_{i,j}$

$$\mathbf{M}(\hat{n}_i \cdot \hat{n}_j) = \begin{pmatrix} \langle T_i T_j \rangle & \langle T_i Q_j \rangle & \langle T_i U_j \rangle \\ \langle T_i Q_j \rangle & \langle Q_i Q_j \rangle & \langle Q_i U_j \rangle \\ \langle T_i U_j \rangle & \langle Q_i U_j \rangle & \langle U_i U_j \rangle \end{pmatrix}. \quad (22)$$

The $(3 \times 3)$ entries in Equation (22) for any given pair of pixels $ij$ depend on the Legendre polynomial $P_{\ell}$ and the fiducial power spectra. As a straightforward example, the entry $\langle T_i T_j \rangle$ is the expression

$$\langle T_i T_j \rangle = \sum_{\ell} \frac{2\ell + 1}{4\pi} P_{\ell}(\hat{r}_i \cdot \hat{r}_j) C_{\ell}^{\mathrm{TT}}. \quad (23)$$

where $P_{\ell}(\hat{r}_i \cdot \hat{r}_j) = \frac{4\pi}{2\ell+1} \sum_m Y_{\ell m}(\hat{r}_i) Y_{\ell m}^*(\hat{r}_j)$. A detailed description of the full procedure to obtain the covariance matrix, together with the expressions of the $(3 \times 3)$ entries, can be found in appendix A of Tegmark and de Oliveira-Costa [30].

It should be noted that the pixel-based likelihood in Equation (21) is exact even in the case of partial sky coverage; this not the case for the likelihood in harmonic space in Equation (17). Note however that it is still possible to derive an exact form for the harmonic-space likelihood even for partial sky coverage [31]. The pixel-based approach ensures mathematical rigor in the evaluation of the likelihood function. Nevertheless, it is highly expensive from a computational point of view. Indeed, the number of pixels needed to retain the information in the first $\ell_{\max}$ multipoles of the power spectrum scales as $\ell_{\max}^2$ and therefore the Cholesky decomposition required to evaluate the inverse of the covariance matrix in Equation (21) scales roughly as $\ell_{\max}^6$, where $\ell_{\max}$ is the highest multipole retained in the analysis. The computational cost is therefore driven by the evaluation of inverse matrix and determinants and becomes prohibitive for $\ell_{\max}$ larger than few hundreds. For this reason, this exact approach is feasible only to study large angular scales, where the information is contained in a relatively small number of multipoles.
### 3. LIKELIHOOD APPROACHES—SMALL-SCALE REGIME

The exact likelihood of the observed CMB $\hat{C}_{\ell}^{XY}$ as a function of the underlying fiducial CMB $C_{\ell}^{XY}$ is given by Equation (17) in case of full-sky observations:

$$p(\mathbf{S}_{\ell}|\mathbf{W}_{\ell}) = \mathcal{L}(\mathbf{W}_{\ell}) = \frac{|\mathbf{S}_{\ell}|^{(\nu-p-1)/2} \exp[-\text{trace}(\mathbf{W}_{\ell}^{-1}\mathbf{S}_{\ell}/2)]}{2^{p\nu/2}|\mathbf{W}_{\ell}|^{\nu/2}\Gamma_{p}(\nu/2)} \quad (24)$$

However, complications arise in real analysis that make it necessary to replace Equation (17) with a suitable approximation. Complications usually include time-consuming evaluations of Equation (17) due to the inversion of large covariance matrices for each theoretical model.

A standard approach is to develop an approximation of Equation (17) in the full-sky regime that is quadratic in some function of $C_{\ell}^{XY}$, and that can be easily generalized to the cut-sky regime with a proper estimate of the covariance matrix:

$$-2 \ln p(\hat{\mathbf{X}}_C | \mathbf{X}_C) \equiv -2 \ln \mathcal{L}(\mathbf{X}_C) \propto \left[ (\mathbf{Z}_C - \hat{\mathbf{Z}}_C)^T \mathbf{Y}^{-1} (\mathbf{Z}_C - \hat{\mathbf{Z}}_C) + \ln |\mathbf{Y}| \right] \quad (25)$$

where $\mathbf{Z}_C$ ($\hat{\mathbf{Z}}_C$) is the vector containing functions of $C_{\ell}S$ ($\hat{C}_{\ell}S$) and $\mathbf{Y}$ is a suitable choice of the covariance matrix.

In what follows, we introduce a list of the most common approximate forms among those proposed in the literature (see e.g., [29, 32–36]). We further quantify the goodness of the approximation in the full-sky regime following the approach in Percival and Brown [29]: we expand the exact likelihood and the approximate forms along the standard axes (TT, TE, EE) around the maximum $\mathbf{X}_C = \hat{\mathbf{X}}_C$, and compare the expansion coefficients up to a certain order.

As already commented in the previous section, the analytic comparison of the various approximations is carried in absence of noise contaminations and in the limit of infinite angular resolution. We also implicitly assume that the CMB spectra are *unlensed*, i.e., they are the spectra as they would be observed in absence of gravitational lensing effects on the CMB photons. The inclusion of experimental noise, experimental angular resolution, and gravitational lensing effects will be discussed in sections 3.3 and 4.

Before moving to the list of the most common approximate likelihood functions, we would like to mention that it is not trivial to construct an unbiased estimator of the true $C_{\ell}$ in the cut-sky regime. We don't have access to the full-sky set of $a_{\ell m}$ and therefore we cannot directly construct $\hat{C}_{\ell}$. In the case of cut-sky maps, appropriate algorithms have been developed to derive the unbiased estimator $\hat{C}_{\ell}$ to be used in the likelihood analysis. For example, pseudo-$C_{\ell}$ power spectra $\tilde{C}_{\ell}$ can be defined from cut-sky harmonic coefficients $\tilde{a}_{\ell m}$, see section B for further details and references. The pseudo-$C_{\ell}$ are related to the true $C_{\ell}$ in ensemble average as

$$\langle \tilde{C}_{\ell'} \rangle = \sum_{\ell} M_{\ell \ell'} \langle C_{\ell} \rangle \quad (26)$$

where $M_{\ell_1 \ell_2}$ is a coupling matrix that encodes the geometrical effects of cut-sky observations. From Equation (26), it is possible to operatively define an estimator for the $C_{\ell}$ in the cut-sky regime as:

$$\hat{C}_{\ell} = \sum_{\ell'} M_{\ell \ell'}^{-1} \tilde{C}_{\ell'} . \quad (27)$$

The interested reader can find a detailed discussion in **Appendix B**, where we also report alternative methods adopted to construct estimators for the BB spectrum B.2, and for power spectra at large scales via the quadratic maximum likelihood (QML) approach B3. A final remark on the cut-sky case: the compression of information from CMB maps ($\sim (N_{\text{pix}} \times N_{\text{pix}})$ pixels) to CMB spectra ($\sim (\ell_{\text{max}} - \ell_{\text{min}})$ bandpowers) is lossless in the full-sky regime, i.e., the power spectra represent sufficient statistics. In the cut-sky regime, the compression is partly lossy, as the masked regions induce correlations between multipoles which have to be taken into account (see e.g., discussion in **Appendix B**).

#### 3.1. Approximate Forms

The most common approximations are given by quadratic/Gaussian expressions in $C_{\ell}^{XY}$ with different choices for the covariance matrix. Alternatively, quadratic expressions involving more complicated functions $f = f(C_{\ell}^{XY}) \equiv \mathbf{Z}_C$, as well as *ad-hoc* combinations of various approximations have been developed to match the exact likelihood up to a certain order in the perturbative regime (see section 3.2).

* **Symmetric Gaussian.** This approximation is quadratic in $\mathbf{Z}_C = C_{\ell}^{XY}$, with covariance matrix given by the curvature of the Wishart, see Equation (29):

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C) (\mathbf{X}_C - \mathbf{X}_C) + \text{const.} \quad (28)$$

The inverse of the covariance matrix $\mathbf{Y}_C$ is the curvature of the Wishart distribution in Equation (19a), i.e., $\mathbf{Y}_{ij}^{-1} = d^2(-2 \ln p)_{ij}|_{\hat{C}_{\ell} = \hat{C}_{\ell'}}$, computed in $\mathbf{X}_C = \hat{\mathbf{X}}_C$:

$$\mathbf{Y}_{\ell}^{-1} = \frac{\nu}{2 \left[ \hat{C}_{\ell}^{TT} \hat{C}_{\ell}^{EE} - (\hat{C}_{\ell}^{TE})^2 \right]^2} \times \begin{pmatrix} (\hat{C}_{\ell}^{EE})^2 & -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{EE} & (\hat{C}_{\ell}^{TE})^2 \\ -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{EE} & 2[\hat{C}_{\ell}^{TT} \hat{C}_{\ell}^{EE} + (\hat{C}_{\ell}^{TE})^2] & -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{TT} \\ (\hat{C}_{\ell}^{TE})^2 & -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{TT} & (\hat{C}_{\ell}^{TT})^2 \end{pmatrix}, \quad (29)$$

* **Improper Gaussian.** This approximation is similar to the Symmetric Gaussian in Equation (28), with the covariance matrix that appears in the first term replaced by $\mathbf{Y} = \mathbf{Y}(C_{\ell}^{XY})$; i.e., the covariance matrix is given in terms of the model $\hat{C}_{\ell}^{XY}$. This approximation is an improper Gaussian in a sense that there is no determinant term:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C) (\mathbf{X}_C - \mathbf{X}_C) \quad (30)$$

* **Determinant Gaussian.** The expression in Equation (30) can be slightly modified to provide a better fit to the
## 3. LIKELIHOOD APPROACHES—SMALL-SCALE REGIME

The exact likelihood of the observed CMB $\hat{C}_{\ell}^{XY}$ as a function of the underlying fiducial CMB $C_{\ell}^{XY}$ is given by Equation (17) in case of full-sky observations:

$$p(\mathbf{S}_{\ell}|\mathbf{W}_{\ell}) = \mathcal{L}(\mathbf{W}_{\ell}) = \frac{|\mathbf{S}_{\ell}|^{(\nu-p-1)/2} \exp[-\text{trace}(\mathbf{W}_{\ell}^{-1}\mathbf{S}_{\ell}/2)]}{2^{p\nu/2}|\mathbf{W}_{\ell}|^{\nu/2}\Gamma_{p}(\nu/2)} \quad (24)$$

However, complications arise in real analysis that make it necessary to replace Equation (17) with a suitable approximation. Complications usually include time-consuming evaluations of Equation (17) due to the inversion of large covariance matrices for each theoretical model.

A standard approach is to develop an approximation of Equation (17) in the full-sky regime that is quadratic in some function of $C_{\ell}^{XY}$, and that can be easily generalized to the cut-sky regime with a proper estimate of the covariance matrix:

$$-2 \ln p(\hat{\mathbf{X}}_C | \mathbf{X}_C) \equiv -2 \ln \mathcal{L}(\mathbf{X}_C) \propto \left[ (\mathbf{Z}_C - \hat{\mathbf{Z}}_C)^T \mathbf{Y}^{-1} (\mathbf{Z}_C - \hat{\mathbf{Z}}_C) + \ln |\mathbf{Y}| \right] \quad (25)$$

where $\mathbf{Z}_C$ ($\hat{\mathbf{Z}}_C$) is the vector containing functions of $C_{\ell}S$ ($\hat{C}_{\ell}S$) and $\mathbf{Y}$ is a suitable choice of the covariance matrix.

In what follows, we introduce a list of the most common approximate forms among those proposed in the literature (see e.g., [29, 32–36]). We further quantify the goodness of the approximation in the full-sky regime following the approach in Percival and Brown [29]: we expand the exact likelihood and the approximate forms along the standard axes (TT, TE, EE) around the maximum $\mathbf{X}_C = \hat{\mathbf{X}}_C$, and compare the expansion coefficients up to a certain order.

As already commented in the previous section, the analytic comparison of the various approximations is carried in absence of noise contaminations and in the limit of infinite angular resolution. We also implicitly assume that the CMB spectra are *unlensed*, i.e., they are the spectra as they would be observed in absence of gravitational lensing effects on the CMB photons. The inclusion of experimental noise, experimental angular resolution, and gravitational lensing effects will be discussed in sections 3.3 and 4.

Before moving to the list of the most common approximate likelihood functions, we would like to mention that it is not trivial to construct an unbiased estimator of the true $C_{\ell}$ in the cut-sky regime. We don't have access to the full-sky set of $a_{\ell m}$ and therefore we cannot directly construct $\hat{C}_{\ell}$. In the case of cut-sky maps, appropriate algorithms have been developed to derive the unbiased estimator $\hat{C}_{\ell}$ to be used in the likelihood analysis. For example, pseudo-$C_{\ell}$ power spectra $\tilde{C}_{\ell}$ can be defined from cut-sky harmonic coefficients $\tilde{a}_{\ell m}$, see section B for further details and references. The pseudo-$C_{\ell}$ are related to the true $C_{\ell}$ in ensemble average as

$$\langle \tilde{C}_{\ell'} \rangle = \sum_{\ell} M_{\ell' \ell} \langle C_{\ell} \rangle \quad (26)$$

where $M_{\ell_1 \ell_2}$ is a coupling matrix that encodes the geometrical effects of cut-sky observations. From Equation (26), it is possible to operatively define an estimator for the $C_{\ell}$ in the cut-sky regime as:

$$\hat{C}_{\ell} = \sum_{\ell'} M_{\ell \ell'}^{-1} \tilde{C}_{\ell'} . \quad (27)$$

The interested reader can find a detailed discussion in **Appendix B**, where we also report alternative methods adopted to construct estimators for the BB spectrum B.2, and for power spectra at large scales via the quadratic maximum likelihood (QML) approach B3. A final remark on the cut-sky case: the compression of information from CMB maps ($\sim (N_{\text{pix}} \times N_{\text{pix}})$ pixels) to CMB spectra ($\sim (\ell_{\text{max}} - \ell_{\text{min}})$ bandpowers) is lossless in the full-sky regime, i.e., the power spectra represent sufficient statistics. In the cut-sky regime, the compression is partly lossy, as the masked regions induce correlations between multipoles which have to be taken into account (see e.g., discussion in **Appendix B**).

### 3.1. Approximate Forms

The most common approximations are given by quadratic/Gaussian expressions in $C_{\ell}^{XY}$ with different choices for the covariance matrix. Alternatively, quadratic expressions involving more complicated functions $f = f(C_{\ell}^{XY}) \equiv \mathbf{Z}_C$, as well as *ad-hoc* combinations of various approximations have been developed to match the exact likelihood up to a certain order in the perturbative regime (see section 3.2).

* **Symmetric Gaussian.** This approximation is quadratic in $\mathbf{Z}_C = C_{\ell}^{XY}$, with covariance matrix given by the curvature of the Wishart, see Equation (29):

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C) (\mathbf{X}_C - \mathbf{X}_C) + \text{const.} \quad (28)$$

The inverse of the covariance matrix $\mathbf{Y}_C$ is the curvature of the Wishart distribution in Equation (19a), i.e., $\mathbf{Y}_{ij}^{-1} = d^2(-2 \ln p)_{ij}|_{C_{\ell}=\hat{C}_{\ell}}$, computed in $\mathbf{X}_C = \hat{\mathbf{X}}_C$:

$$\mathbf{Y}_{\ell}^{-1} = \frac{\nu}{2 \left[ \hat{C}_{\ell}^{TT} \hat{C}_{\ell}^{EE} - (\hat{C}_{\ell}^{TE})^2 \right]^2} \times \begin{pmatrix} (\hat{C}_{\ell}^{EE})^2 & -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{EE} & (\hat{C}_{\ell}^{TE})^2 \\ -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{EE} & 2[\hat{C}_{\ell}^{TT} \hat{C}_{\ell}^{EE} + (\hat{C}_{\ell}^{TE})^2] & -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{TT} \\ (\hat{C}_{\ell}^{TE})^2 & -2 \hat{C}_{\ell}^{TE} \hat{C}_{\ell}^{TT} & (\hat{C}_{\ell}^{TT})^2 \end{pmatrix}, \quad (29)$$

* **Improper Gaussian.** This approximation is similar to the Symmetric Gaussian in Equation (28), with the covariance matrix that appears in the first term replaced by $\mathbf{Y} = \mathbf{Y}(C_{\ell}^{XY})$; i.e., the covariance matrix is given in terms of the model $\hat{C}_{\ell}^{XY}$. This approximation is an improper Gaussian in a sense that there is no determinant term:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C) (\mathbf{X}_C - \mathbf{X}_C) \quad (30)$$

* **Determinant Gaussian.** The expression in Equation (30) can be slightly modified to provide a better fit to the
exact likelihood approach (see section 3.2). The modification consists in the addition of a $C_{\ell}^{XY}$-dependent determinant term:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\mathbf{X}_C)(\hat{\mathbf{X}}_C - \mathbf{X}_C) + \ln |\mathbf{Y}_C(\mathbf{X}_C)| \quad (31)$$

* **Fiducial Gaussian.** This approximation is similar to Equations (28) and (30), with a constant determinant term (as in Equation 28) and the covariance matrix computed for a given fiducial model (as in Equation 30). The fiducial model for the covariance matrix is however kept fixed, and assumed to be smooth and a close approximation to the underlying model under scrutiny:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\mathbf{X}_{C,\text{fid}})(\hat{\mathbf{X}}_C - \mathbf{X}_C) + \text{const.} \quad (32)$$

The fiducial Gaussian approximation is used in the official analyses of the Planck [37, 38], ACT [9], and SPT [10] collaborations.

* **Log-normal.** This approximation is quadratic in a peculiar function of theoretical and observed spectra, i.e., $Z_C = \hat{C}_{\ell}^{XY} \ln(C_{\ell}^{XY})$, with fixed covariance matrix $\mathbf{Y}_C = \mathbf{Y}_C(\hat{C}_{\ell}^{XY})$:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{Z}}_C - \mathbf{Z}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C)(\hat{\mathbf{Z}}_C - \mathbf{Z}_C) + \text{const.} \quad (33)$$

* **Offset log-normal.** This approximation is a generalization of Equation (33). The data vector is generalized to $Z_C = (1 + a_{XY})\hat{C}_{\ell}^{XY} \ln(C_{\ell}^{XY} + a_{XY}\hat{C}_{\ell}^{XY})$, with $a_{XY}$ a suitable real offset coefficient that may or may be not be the same for every XY pair. The covariance matrix is again as in Equation (33) ($\mathbf{Y}_C = \mathbf{Y}_C(\hat{C}_{\ell}^{XY})$):

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{Z}}_C - \mathbf{Z}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C)(\hat{\mathbf{Z}}_C - \mathbf{Z}_C) + \text{const.} \quad (34)$$

* **One-third-two-thirds.** We briefly mention this approximation as an example of combined likelihood appositely built to match the exact likelihood up to the third order. It is a weighted combination of the improper Gaussian in Equation (30) (with weight 1/3) and of the log-normal approximation in Equation (33) (with weight 2/3). Note that the approximation was explicitly built for the single-field TT-only WMAP analysis [39]:

$$\ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{1}{3} \ln [\mathcal{L}_{\text{improper}}(C_{\ell}^{TT})] + \frac{2}{3} \ln [\mathcal{L}_{\text{lognorm}}(C_{\ell}^{TT})] \quad (35)$$

* **Hamimeche-Lewis.** In Hamimeche and Lewis [32], Hamimeche & Lewis (HL) have developed a form of the likelihood for correlated Gaussian fields (CMB temperature and polarization) that coincides with the exact likelihood in full sky. The authors show with simulations that it provides a very good approximation to the exact likelihood in the cut-sky regime at small scales<fnref n="4" /> ($\ell \ge 30$). The form of the likelihood is quadratic in some peculiar function of the observed, fiducial, and theoretical $C_{\ell}$, as we shall see in section 3.2. The covariance matrix is precomputed for a fixed fiducial model.

The dependence on the fiducial model is negligible. Moreover, should the fiducial fail in matching the true sky, the likelihood is still exact in full sky. The HL likelihood was used in the analysis of the BICEP2/KECK data [11].

In the HL formalism, the likelihood in cut-sky can be approximated as:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \approx \mathbf{X}_{g\ell}^T \mathbf{M}_{f\ell}^{-1} \mathbf{X}_{g\ell} \quad (36)$$

where $\mathbf{X}_g$ is a vector of a specific function of the observed, fiducial, and theoretical $C_{\ell}$, and $\mathbf{M}_f$ is the fiducial model covariance block matrix

$$[\mathbf{M}_f]_{\ell\ell'} = \langle (\hat{\mathbf{X}}_\ell - \mathbf{X}_\ell)(\hat{\mathbf{X}}_{\ell'} - \mathbf{X}_{\ell'})^T \rangle_f \quad (37)$$

with $n(n+1)/2 \times n(n+1)/2$ blocks ($n$ is the number of fields), labeled by $\ell$ and $\ell'$, i.e., we explicitly take into account the possibility that either the cut-sky or anisotropic noise can induce correlations between different multipoles (non-diagonal covariance).

The derivation of Equation (36) is provided in section 3.2, where we will show that it is equivalent to the exact likelihood in the full-sky regime. For more details and a thorough definition of the notation, see Hamimeche and Lewis [32]

## 3.2. Comparison With the Exact Likelihood in the Full-Sky Regime

In this section, we comment on the goodness of the approximations listed in the previous section. The goodness is defined in terms of the ability to match the exact likelihood in the full-sky regime up to a certain order, when both the exact likelihood and the approximate form are expanded around the maximum. This approach is described in Percival and Brown [29].

Let's start by expanding the Wishart distribution in Equation (19a) along the TT direction. In particular, we write Equation (19a) with the following substitutions: $C_{\ell}^{TE} \rightarrow \hat{C}_{\ell}^{TE}$, $C_{\ell}^{EE} \rightarrow \hat{C}_{\ell}^{EE}$, $C_{\ell}^{TT} \rightarrow (1+\epsilon)\hat{C}_{\ell}^{TT}$. We further expand in $\epsilon$. We obtain:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) = \frac{\nu}{(\hat{r}^2 - 1)^2} \left[ \frac{\epsilon^2}{2} + \frac{2}{3(\hat{r}^2 - 1)} \epsilon^3 + \mathcal{O}(\epsilon^4) \right] \quad (38)$$

where $\hat{r} = \hat{C}_{\ell}^{TE}/\sqrt{\hat{C}_{\ell}^{TT}\hat{C}_{\ell}^{EE}}$ and $\nu = 2\ell+1$. It is straightforward to show that the expansion along the EE axis provides the same form of Equation (38) for $-2\ln[\mathcal{L}(C_{\ell}^{EE})]$.

Let's now expand along the TE axis with the following substitutions: $C_{\ell}^{TT} \rightarrow \hat{C}_{\ell}^{TT}$, $C_{\ell}^{EE} \rightarrow \hat{C}_{\ell}^{EE}$, $C_{\ell}^{TE} \rightarrow (1+\epsilon)\hat{C}_{\ell}^{TE}$. The expansion is:

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) = \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left[ \frac{\epsilon^2}{2} - \frac{2\hat{r}^2}{3(\hat{r}^2 - 1)} \epsilon^3 + \mathcal{O}(\epsilon^4) \right] \quad (39)$$

Equation (39) reflects again the difference in the marginal distribution of the cross-correlation spectrum TE with respect to

<footnote n="4">⁴One of the assumptions is that the matrix of the estimators $\hat{C}_{\ell}$ is positive definite. This assumption may break up at large scales.</footnote>
exact likelihood approach (see section 3.2). The modification consists in the addition of a $C_{\ell}^{XY}$-dependent determinant term:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\mathbf{X}_C)(\hat{\mathbf{X}}_C - \mathbf{X}_C) + \ln |\mathbf{Y}_C(\mathbf{X}_C)| \quad (31)$$

* **Fiducial Gaussian.** This approximation is similar to Equations (28) and (30), with a constant determinant term (as in Equation 28) and the covariance matrix computed for a given fiducial model (as in Equation 30). The fiducial model for the covariance matrix is however kept fixed, and assumed to be smooth and a close approximation to the underlying model under scrutiny:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{X}}_C - \mathbf{X}_C)^T \mathbf{Y}_C^{-1} (\mathbf{X}_{C,\text{fid}})(\hat{\mathbf{X}}_C - \mathbf{X}_C) + \text{const.} \quad (32)$$

The fiducial Gaussian approximation is used in the official analyses of the Planck [37, 38], ACT [9], and SPT [10] collaborations.

* **Log-normal.** This approximation is quadratic in a peculiar function of theoretical and observed spectra, i.e., $Z_C = \hat{C}_{\ell}^{XY} \ln(C_{\ell}^{XY})$, with fixed covariance matrix $\mathbf{Y}_C = \mathbf{Y}_C(\hat{C}_{\ell}^{XY})$:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{Z}}_C - \mathbf{Z}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C)(\hat{\mathbf{Z}}_C - \mathbf{Z}_C) + \text{const.} \quad (33)$$

* **Offset log-normal.** This approximation is a generalization of Equation (33). The data vector is generalized to $Z_C = (1 + a_{XY})\hat{C}_{\ell}^{XY} \ln(C_{\ell}^{XY} + a_{XY}\hat{C}_{\ell}^{XY})$, with $a_{XY}$ a suitable real offset coefficient that may or may be not be the same for every XY pair. The covariance matrix is again as in Equation (33) ($\mathbf{Y}_C = \mathbf{Y}_C(\hat{C}_{\ell}^{XY})$):

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \propto (\hat{\mathbf{Z}}_C - \mathbf{Z}_C)^T \mathbf{Y}_C^{-1} (\hat{\mathbf{X}}_C)(\hat{\mathbf{Z}}_C - \mathbf{Z}_C) + \text{const.} \quad (34)$$

* **One-third-two-thirds.** We briefly mention this approximation as an example of combined likelihood appositely built to match the exact likelihood up to the third order. It is a weighted combination of the improper Gaussian in Equation (30) (with weight 1/3) and of the log-normal approximation in Equation (33) (with weight 2/3). Note that the approximation was explicitly built for the single-field TT-only WMAP analysis [39]:

$$\ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{1}{3} \ln [\mathcal{L}_{\text{improper}}(C_{\ell}^{TT})] + \frac{2}{3} \ln [\mathcal{L}_{\text{lognorm}}(C_{\ell}^{TT})] \quad (35)$$

* **Hamimeche-Lewis.** In Hamimeche and Lewis [32], Hamimeche & Lewis (HL) have developed a form of the likelihood for correlated Gaussian fields (CMB temperature and polarization) that coincides with the exact likelihood in full sky. The authors show with simulations that it provides a very good approximation to the exact likelihood in the cut-sky regime at small scales<fnref n="4" /> ($\ell \ge 30$). The form of the likelihood is quadratic in some peculiar function of the observed, fiducial, and theoretical $C_{\ell}$, as we shall see in section 3.2. The covariance matrix is precomputed for a fixed fiducial model.

The dependence on the fiducial model is negligible. Moreover, should the fiducial fail in matching the true sky, the likelihood is still exact in full sky. The HL likelihood was used in the analysis of the BICEP2/KECK data [11].

In the HL formalism, the likelihood in cut-sky can be approximated as:

$$-2 \ln \mathcal{L}(\mathbf{X}_C) \approx \mathbf{X}_{g\ell}^T \mathbf{M}_{f\ell}^{-1} \mathbf{X}_{g\ell} \quad (36)$$

where $\mathbf{X}_g$ is a vector of a specific function of the observed, fiducial, and theoretical $C_{\ell}$, and $\mathbf{M}_f$ is the fiducial model covariance block matrix

$$[\mathbf{M}_f]_{\ell\ell'} = \langle (\hat{\mathbf{X}}_\ell - \mathbf{X}_\ell)(\hat{\mathbf{X}}_{\ell'} - \mathbf{X}_{\ell'})^T \rangle_f \quad (37)$$

with $n(n+1)/2 \times n(n+1)/2$ blocks ($n$ is the number of fields), labeled by $\ell$ and $\ell'$, i.e., we explicitly take into account the possibility that either the cut-sky or anisotropic noise can induce correlations between different multipoles (non-diagonal covariance).

The derivation of Equation (36) is provided in section 3.2, where we will show that it is equivalent to the exact likelihood in the full-sky regime. For more details and a thorough definition of the notation, see Hamimeche and Lewis [32]

## 3.2. Comparison With the Exact Likelihood in the Full-Sky Regime

In this section, we comment on the goodness of the approximations listed in the previous section. The goodness is defined in terms of the ability to match the exact likelihood in the full-sky regime up to a certain order, when both the exact likelihood and the approximate form are expanded around the maximum. This approach is described in Percival and Brown [29].

Let's start by expanding the Wishart distribution in Equation (19a) along the TT direction. In particular, we write Equation (19a) with the following substitutions: $C_{\ell}^{TE} \rightarrow \hat{C}_{\ell}^{TE}$, $C_{\ell}^{EE} \rightarrow \hat{C}_{\ell}^{EE}$, $C_{\ell}^{TT} \rightarrow (1+\epsilon)\hat{C}_{\ell}^{TT}$. We further expand in $\epsilon$. We obtain:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) = \frac{\nu}{(\hat{r}^2 - 1)^2} \left[ \frac{\epsilon^2}{2} + \frac{2}{3(\hat{r}^2 - 1)} \epsilon^3 + \mathcal{O}(\epsilon^4) \right] \quad (38)$$

where $\hat{r} = \hat{C}_{\ell}^{TE}/\sqrt{\hat{C}_{\ell}^{TT}\hat{C}_{\ell}^{EE}}$ and $\nu = 2\ell+1$. It is straightforward to show that the expansion along the EE axis provides the same form of Equation (38) for $-2\ln[\mathcal{L}(C_{\ell}^{EE})]$.

Let's now expand along the TE axis with the following substitutions: $C_{\ell}^{TT} \rightarrow \hat{C}_{\ell}^{TT}$, $C_{\ell}^{EE} \rightarrow \hat{C}_{\ell}^{EE}$, $C_{\ell}^{TE} \rightarrow (1+\epsilon)\hat{C}_{\ell}^{TE}$. The expansion is:

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) = \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left[ \frac{\epsilon^2}{2} - \frac{2\hat{r}^2}{3(\hat{r}^2 - 1)} \epsilon^3 + \mathcal{O}(\epsilon^4) \right] \quad (39)$$

Equation (39) reflects again the difference in the marginal distribution of the cross-correlation spectrum TE with respect to

<footnote n="4">⁴One of the assumptions is that the matrix of the estimators $\hat{C}_{\ell}$ is positive definite. This assumption may break up at large scales.</footnote>
the distributions of the auto-spectra TT and EE (see discussion in section 2.2).

Let's now move to expand each of the approximations reported in section 3.1.

* **Symmetric Gaussian.** The expansion of Equation (28) is of the same form along each of the standard axes TT, TE, EE, with a different normalization factor in the case of the expansion along TE:

$$ -2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2}\right) \quad (40a) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2}\right) \quad (40b) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2}\right) \quad (40c) $$

Note that the expansion is truncated at the second order on $\epsilon$. This result is an exact expansion and it is expected, given the initial form (symmetric Gaussian) of the approximate likelihood. If we compare Equation (40a) with the expansion of the Wishart in Equations (38)–(39), we observe what follows. Firstly, the Symmetric Gaussian matches the Wishart only up to the second order on $\epsilon$. Secondly, the approximate form is, by definition, symmetric in $\epsilon$ and therefore fails in capturing the skewness of the exact likelihood. Finally, it is biased low along the TT, EE axes in a sense that $-2\ln(\mathcal{L}_{\text{gaussian}}(\epsilon)) < -2\ln(\mathcal{L}_{\text{exact}}(\epsilon))$ for $\epsilon > 0$ (i.e., for $C_{\ell} > \hat{C}_{\ell}$). For the opposite reason, it is biased high along the TE axis.

* **Improper Gaussian.** The expansion of Equation (30) along the standard axes are as follows:

$$ -2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right) \quad (41a) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right) \quad (41b) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\hat{r}^2(\hat{r}^2 + 3)\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right) \quad (41c) $$

With respect to the symmetric Gaussian, the improper Gaussian approximation is skewed in the same direction of the Wishart. However, it is still a correct match only up to the second order. With respect to the exact likelihood, Equation (41a) show that the improper Gaussian is biased high along the TT, EE directions, where $-2\ln(\mathcal{L}_{\text{improper}}(\epsilon)) > -2\ln(\mathcal{L}_{\text{exact}}(\epsilon))$ for $\epsilon > 0$ (i.e., for $C_{\ell} > \hat{C}_{\ell}$). For the opposite reason, it is biased low along the TE direction.

* **Determinant Gaussian.** In this case, it is clear from Equation (31) that the likelihood is biased at each multipole $\ell$ as the $C_{\ell}$-dependent determinant term implies that the minimum value for this approximate form is not in $\epsilon = 0$. Indeed, the expansions of Equation (31) along the standard axes include a term of order $\epsilon$:

$$ -2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{1}{(\hat{r}^2 - 1)^2} \left[ -3\epsilon(\hat{r}^2 - 1) + \frac{(\nu - 3)\epsilon^2}{2} \right] $$

$$ + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \quad (42a) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{1}{(\hat{r}^2 - 1)^2} \left[ -3\epsilon(\hat{r}^2 - 1) + \frac{(\nu - 3)\epsilon^2}{2} + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right] \quad (42b) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2}{(\hat{r}^2 - 1)^2} \left[ 3\epsilon(\hat{r}^2 - 1) + \frac{(\nu - 3)(\hat{r}^2 + 1)\epsilon^2}{2} - \frac{(\nu - 1)\hat{r}^2(\hat{r}^2 + 1)\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right] \quad (42c) $$

However, it can be shown that this approximation, although biased at each individual $\ell$, is unbiased “on average”, i.e., reproduces the correct result with reasonable accuracy when summed over a wide-enough range of multiples (see HL).

* **Fiducial Gaussian.** The expansion of Equation (32) is equivalent to the expansion in Equation (40a), only with a different normalization factor. Indeed, the covariance matrix in Equation (32) is fixed to that of a given fiducial model, and therefore the approximation is quadratic in $\epsilon$. Note however that, although the form of the expansion is similar at any $\ell$ between the symmetric Gaussian and the fiducial Gaussian, the latter provides a better approximation of the exact likelihood when summed over a range of multipoles (see discussion in Hamimeche and Lewis [32]).

* **Log-normal.** In this case, we have:

$$ Z_C = \begin{pmatrix} \hat{C}_{\ell}^{TT} \ln[C_{\ell}^{TT}] \\ \hat{C}_{\ell}^{TE} \ln[C_{\ell}^{TE}] \\ \hat{C}_{\ell}^{EE} \ln[C_{\ell}^{EE}] \end{pmatrix}, \quad (43) $$

in Equation (25), with $Y_C = Y_C(\hat{C}_{\ell}^{XY})$ being the curvature matrix. The expansions along the auto- and cross-spectra directions can be easily obtained up to normalization factors:

$$ -2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2} + \mathcal{O}(\epsilon^4) \right) \quad (44a) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2} + \mathcal{O}(\epsilon^4) \right) \quad (44b) $$

$$ -2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2} + \mathcal{O}(\epsilon^4) \right) \quad (44c) $$

Regardless of the normalization factors, a comparison between Equations (44a)–(44c) and Equations (38) and (39) shows that the log-normal distribution provides a good approximation to the exact likelihood up to the second order in the expansion around the maximum. The two distributions have a different shape starting from the third-order term in the series expansion. It is also interesting to note that, normalization
the distributions of the auto-spectra TT and EE (see discussion in section 2.2).

Let's now move to expand each of the approximations reported in section 3.1.

* **Symmetric Gaussian.** The expansion of Equation (28) is of the same form along each of the standard axes TT, TE, EE, with a different normalization factor in the case of the expansion along TE:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2}\right) \quad (40a)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2}\right) \quad (40b)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2}\right) \quad (40c)$$

Note that the expansion is truncated at the second order on $\epsilon$. This result is an exact expansion and it is expected, given the initial form (symmetric Gaussian) of the approximate likelihood. If we compare Equation (40a) with the expansion of the Wishart in Equations (38)–(39), we observe what follows. Firstly, the Symmetric Gaussian matches the Wishart only up to the second order on $\epsilon$. Secondly, the approximate form is, by definition, symmetric in $\epsilon$ and therefore fails in capturing the skewness of the exact likelihood. Finally, it is biased low along the TT, EE axes in a sense that $-2\ln(\mathcal{L}_{\text{gaussian}}(\epsilon)) < -2\ln(\mathcal{L}_{\text{exact}}(\epsilon))$ for $\epsilon > 0$ (i.e., for $C_{\ell} > \hat{C}_{\ell}$). For the opposite reason, it is biased high along the TE axis.

* **Improper Gaussian.** The expansion of Equation (30) along the standard axes are as follows:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2} + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4)\right) \quad (41a)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2} + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4)\right) \quad (41b)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left(\frac{\epsilon^2}{2} - \frac{\hat{r}^2(\hat{r}^2 + 3)\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4)\right) \quad (41c)$$

With respect to the symmetric Gaussian, the improper Gaussian approximation is skewed in the same direction of the Wishart. However, it is still a correct match only up to the second order. With respect to the exact likelihood, Equation (41a) show that the improper Gaussian is biased high along the TT, EE directions, where $-2\ln(\mathcal{L}_{\text{improper}}(\epsilon)) > -2\ln(\mathcal{L}_{\text{exact}}(\epsilon))$ for $\epsilon > 0$ (i.e., for $C_{\ell} > \hat{C}_{\ell}$). For the opposite reason, it is biased low along the TE direction.

* **Determinant Gaussian.** In this case, it is clear from Equation (31) that the likelihood is biased at each multipole $\ell$ as the $C_{\ell}$-dependent determinant term implies that the minimum value for this approximate form is not in $\epsilon = 0$. Indeed, the expansions of Equation (31) along the standard axes include a term of order $\epsilon$:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{1}{(\hat{r}^2 - 1)^2} \left[ -3\epsilon(\hat{r}^2 - 1) + \frac{(\nu - 3)\epsilon^2}{2} \right]$$

$$+ \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \quad (42a)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{1}{(\hat{r}^2 - 1)^2} \left[ -3\epsilon(\hat{r}^2 - 1) + \frac{(\nu - 3)\epsilon^2}{2} + \frac{\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right] \quad (42b)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2}{(\hat{r}^2 - 1)^2} \left[ 3\epsilon(\hat{r}^2 - 1) + \frac{(\nu - 3)(\hat{r}^2 + 1)\epsilon^2}{2} - \frac{(\nu - 1)\hat{r}^2(\hat{r}^2 + 1)\epsilon^3}{(\hat{r}^2 - 1)} + \mathcal{O}(\epsilon^4) \right] \quad (42c)$$

However, it can be shown that this approximation, although biased at each individual $\ell$, is unbiased “on average”, i.e., reproduces the correct result with reasonable accuracy when summed over a wide-enough range of multiples (see HL).

* **Fiducial Gaussian.** The expansion of Equation (32) is equivalent to the expansion in Equation (40a), only with a different normalization factor. Indeed, the covariance matrix in Equation (32) is fixed to that of a given fiducial model, and therefore the approximation is quadratic in $\epsilon$. Note however that, although the form of the expansion is similar at any $\ell$ between the symmetric Gaussian and the fiducial Gaussian, the latter provides a better approximation of the exact likelihood when summed over a range of multipoles (see discussion in Hamimeche and Lewis [32]).

* **Log-normal.** In this case, we have:

$$Z_C = \begin{pmatrix} \hat{C}_{\ell}^{TT} \ln[C_{\ell}^{TT}] \\ \hat{C}_{\ell}^{TE} \ln[C_{\ell}^{TE}] \\ \hat{C}_{\ell}^{EE} \ln[C_{\ell}^{EE}] \end{pmatrix}, \quad (43)$$

in Equation (25), with $Y_C = Y_C(\hat{C}_{\ell}^{XY})$ being the curvature matrix. The expansions along the auto- and cross-spectra directions can be easily obtained up to normalization factors:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2} + \mathcal{O}(\epsilon^4) \right) \quad (44a)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2} + \mathcal{O}(\epsilon^4) \right) \quad (44b)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2} + \mathcal{O}(\epsilon^4) \right) \quad (44c)$$

Regardless of the normalization factors, a comparison between Equations (44a)–(44c) and Equations (38) and (39) shows that the log-normal distribution provides a good approximation to the exact likelihood up to the second order in the expansion around the maximum. The two distributions have a different shape starting from the third-order term in the series expansion. It is also interesting to note that, normalization
factors aside, the expansions along the standard axes are identical. This is a further difference with respect to the case of the Wishart distribution.

* **Offset log-normal.** The log-normal distribution can be slightly modified in a way that it could approximate the exact Wishart distribution up to third order. The modified log-normal, or offset log-normal, is quadratic in:

$$Z_C = \begin{pmatrix} \hat{C}_{\ell}^{TT} (1 + a_{TT}) \ln[C_{\ell}^{TT} + a_{TT}\hat{C}_{\ell}^{TT}] \\ \hat{C}_{\ell}^{TE} (1 + a_{TE}) \ln[C_{\ell}^{TE} + a_{TE}\hat{C}_{\ell}^{TE}] \\ \hat{C}_{\ell}^{EE} (1 + a_{EE}) \ln[C_{\ell}^{EE} + a_{EE}\hat{C}_{\ell}^{EE}] \end{pmatrix}, \quad (45)$$

where the offset factors $a_{TT}$, $a_{TE}$, $a_{EE}$ can be adjusted to match the Wishart distribution up to the third order. The covariance matrix is again assumed to be $Y_C = Y_C(\hat{C}_{\ell}^{XY})$. Expanding the offset log-normal in the usual way, one gets:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2(1 + a_{TT})} + \mathcal{O}(\epsilon^4) \right) \quad (46a)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2(1 + a_{EE})} + \mathcal{O}(\epsilon^4) \right) \quad (46b)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2(1 + a_{TE})} + \mathcal{O}(\epsilon^4) \right) \quad (46c)$$

A comparison between Equations (44a)-(44c) and Equations (38, 39) makes it clear that the offset log-normal distributions is a good approximation to the Wishart distribution up to the third order in the expansion, provided that

$$a_{TT} = a_{EE} = -\frac{1}{4}(1 + 3\hat{r}^2), \quad a_{TE} = -\frac{1}{2}\left(2 + \frac{3(\hat{r}^4 - 1)}{2\hat{r}^2(\hat{r}^2 + 3)}\right) \quad (47)$$

* **One-third two-thirds.** Comparing the TT expansion in Equation (41a) and in Equation (44a) with the TT expansion of the exact likelihood in Equation (38), it is clear that the weighted sum of the improper Gaussian and the lognormal distribution with weights 1/3 and 2/3, respectively matches the Wishart distribution up to the third order in $\epsilon$:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) = \left[ \frac{1}{3}(-2 \ln \mathcal{L}_{\text{improper}}) + \frac{2}{3}(-2 \ln \mathcal{L}_{\text{lognorm}}) \right] \\ \propto \nu \left[ \frac{\epsilon^2}{2} - \frac{2\epsilon^3}{3} + \mathcal{O}(\epsilon^4) \right] \quad (48)$$

* **Hamimeche-Lewis.** By construction, this likelihood approximation matches exactly the Wishart distribution in the full sky regime. Indeed, the true power of this approximation stands in the fact that the exact quadratic form derived from the full-sky exact likelihood result is assumed to be valid also in the cut-sky regime and at high multipoles, where it is faster to evaluate than the exact calculation.

We show here the equivalence between the exact likelihood in Equation (17) and the Hamimeche-Lewis formalism in full sky. In what follows, we make use of the matrix notation adopted in Hamimeche and Lewis [32]. This notation is somehow different from the formalism used in the previous examples. However, it is a more suitable choice to better appreciate the H&L approximation. We assume the matrix of the estimators $\hat{C}_{\ell}$ to be positive definite. In the full-sky limit, given $n$ gaussian fields, the likelihood function is defined as in Equation (16a)

$$-2 \ln \mathcal{L}(C_{\ell}) = (2\ell + 1) \times [\operatorname{trace}[\hat{C}_{\ell} C_{\ell}^{-1}] - \ln |C_{\ell}^{-1} \hat{C}_{\ell}| - n] \quad (49a)$$

$$= (2\ell + 1) \times [\operatorname{trace}[C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}] - \ln |C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}| - n] \quad (49b)$$

$$= (2\ell + 1) \sum_{i}^{n} [D_{\ell,ii} - \ln(D_{\ell,ii}) - 1]. \quad (49c)$$

where<fnref n="5" />, with respect to Equation (16a), $S_{\ell} \rightarrow \hat{C}_{\ell}$ and $V_{\ell} \rightarrow C_{\ell}$. In passing from Equations (49a) to (49c), we consider that the symmetric form is defined using the Hermitian square root and $C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2} = U_{\ell} D_{\ell} U_{\ell}^{T}$, for orthogonal $U_{\ell}$ and diagonal $D_{\ell}$. In other words, we diagonalize the matrix $C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}$.

In order to generalize Equation (49c) to the cut sky, we want to reshuffle it in such a way that it resembles a quadratic form:

$$-2 \ln \mathcal{L}(C_{\ell}) = \frac{2\ell + 1}{2} \sum_{i} [g(D_{\ell,ii})]^2 \quad (50a)$$

$$= \frac{2\ell + 1}{2} \operatorname{Tr}[g(D_{\ell})^2]. \quad (50b)$$

where the function $g(x)$ is defined as

$$g(x) \equiv \operatorname{sign}(x - 1)\sqrt{2(x - \ln(x) - 1)}, \quad (51)$$

and $[g(D_{\ell})]_{ij} = g(D_{\ell,ii})\delta_{ij}$.

In order to transform Equation (50b) in a form that is quadratic also in the matrix elements, we exploit the following matrix identity<fnref n="6" />:

$$\frac{2\ell + 1}{2} \operatorname{Trace}[(C_{f\ell}^{-1/2} C_{g\ell} C_{f\ell}^{-1/2})^2] = X_{g\ell}^T M_{f\ell}^{-1} X_{g\ell}, \quad (53)$$

$$\operatorname{Trace}[g(D_{\ell})] = \operatorname{Trace}[U g(D_{\ell}) U^{T}] \quad (52a)$$

$$= \operatorname{Trace}[C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}] \quad (52b)$$

$$= \operatorname{Trace}[C_{f\ell}^{-1/2} C_{f\ell}^{1/2} (C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}) C_{f\ell}^{1/2} C_{f\ell}^{-1/2}] \quad (52c)$$

$$= \operatorname{Trace}[C_{f\ell}^{-1/2} C_{g\ell} C_{f\ell}^{-1/2}] \quad (52d)$$

<footnote n="5">⁵In passing from Equations (49a) to (49b), we have made use of the properties of the trace[] operator.</footnote>

<footnote n="6">⁶Using the invariance of Tr[A] under diagonalization of A, one has that</footnote>
factors aside, the expansions along the standard axes are identical. This is a further difference with respect to the case of the Wishart distribution.

* **Offset log-normal.** The log-normal distribution can be slightly modified in a way that it could approximate the exact Wishart distribution up to third order. The modified log-normal, or offset log-normal, is quadratic in:

$$Z_C = \begin{pmatrix} \hat{C}_{\ell}^{TT} (1 + a_{TT}) \ln[C_{\ell}^{TT} + a_{TT}\hat{C}_{\ell}^{TT}] \\ \hat{C}_{\ell}^{TE} (1 + a_{TE}) \ln[C_{\ell}^{TE} + a_{TE}\hat{C}_{\ell}^{TE}] \\ \hat{C}_{\ell}^{EE} (1 + a_{EE}) \ln[C_{\ell}^{EE} + a_{EE}\hat{C}_{\ell}^{EE}] \end{pmatrix}, \quad (45)$$

where the offset factors $a_{TT}$, $a_{TE}$, $a_{EE}$ can be adjusted to match the Wishart distribution up to the third order. The covariance matrix is again assumed to be $Y_C = Y_C(\hat{C}_{\ell}^{XY})$. Expanding the offset log-normal in the usual way, one gets:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2(1 + a_{TT})} + \mathcal{O}(\epsilon^4) \right) \quad (46a)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{EE}) \propto \frac{\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2(1 + a_{EE})} + \mathcal{O}(\epsilon^4) \right) \quad (46b)$$

$$-2 \ln \mathcal{L}(C_{\ell}^{TE}) \propto \frac{2\hat{r}^2(\hat{r}^2 + 1)\nu}{(\hat{r}^2 - 1)^2} \left( \frac{\epsilon^2}{2} - \frac{\epsilon^3}{2(1 + a_{TE})} + \mathcal{O}(\epsilon^4) \right) \quad (46c)$$

A comparison between Equations (44a)-(44c) and Equations (38, 39) makes it clear that the offset log-normal distributions is a good approximation to the Wishart distribution up to the third order in the expansion, provided that

$$a_{TT} = a_{EE} = -\frac{1}{4}(1 + 3\hat{r}^2), \quad a_{TE} = -\frac{1}{2}\left(2 + \frac{3(\hat{r}^4 - 1)}{2\hat{r}^2(\hat{r}^2 + 3)}\right) \quad (47)$$

* **One-third two-thirds.** Comparing the TT expansion in Equation (41a) and in Equation (44a) with the TT expansion of the exact likelihood in Equation (38), it is clear that the weighted sum of the improper Gaussian and the lognormal distribution with weights 1/3 and 2/3, respectively matches the Wishart distribution up to the third order in $\epsilon$:

$$-2 \ln \mathcal{L}(C_{\ell}^{TT}) = \left[ \frac{1}{3}(-2 \ln \mathcal{L}_{\text{improper}}) + \frac{2}{3}(-2 \ln \mathcal{L}_{\text{lognorm}}) \right] \\ \propto \nu \left[ \frac{\epsilon^2}{2} - \frac{2\epsilon^3}{3} + \mathcal{O}(\epsilon^4) \right] \quad (48)$$

* **Hamimeche-Lewis.** By construction, this likelihood approximation matches exactly the Wishart distribution in the full sky regime. Indeed, the true power of this approximation stands in the fact that the exact quadratic form derived from the full-sky exact likelihood result is assumed to be valid also in the cut-sky regime and at high multipoles, where it is faster to evaluate than the exact calculation.

We show here the equivalence between the exact likelihood in Equation (17) and the Hamimeche-Lewis formalism in full sky. In what follows, we make use of the matrix notation adopted in Hamimeche and Lewis [32]. This notation is somehow different from the formalism used in the previous examples. However, it is a more suitable choice to better appreciate the H&L approximation. We assume the matrix of the estimators $\hat{C}_{\ell}$ to be positive definite. In the full-sky limit, given $n$ gaussian fields, the likelihood function is defined as in Equation (16a)

$$-2 \ln \mathcal{L}(C_{\ell}) = (2\ell + 1) \times [\operatorname{trace}[\hat{C}_{\ell} C_{\ell}^{-1}] - \ln |C_{\ell}^{-1} \hat{C}_{\ell}| - n] \quad (49a)$$

$$= (2\ell + 1) \times [\operatorname{trace}[C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}] - \ln |C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}| - n] \quad (49b)$$

$$= (2\ell + 1) \sum_{i}^{n} [D_{\ell,ii} - \ln(D_{\ell,ii}) - 1]. \quad (49c)$$

where<fnref n="5" />, with respect to Equation (16a), $S_{\ell} \rightarrow \hat{C}_{\ell}$ and $V_{\ell} \rightarrow C_{\ell}$. In passing from Equations (49a) to (49c), we consider that the symmetric form is defined using the Hermitian square root and $C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2} = U_{\ell} D_{\ell} U_{\ell}^{T}$, for orthogonal $U_{\ell}$ and diagonal $D_{\ell}$. In other words, we diagonalize the matrix $C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}$.

In order to generalize Equation (49c) to the cut sky, we want to reshuffle it in such a way that it resembles a quadratic form:

$$-2 \ln \mathcal{L}(C_{\ell}) = \frac{2\ell + 1}{2} \sum_{i} [g(D_{\ell,ii})]^2 \quad (50a)$$

$$= \frac{2\ell + 1}{2} \operatorname{Tr}[g(D_{\ell})^2]. \quad (50b)$$

where the function $g(x)$ is defined as

$$g(x) \equiv \operatorname{sign}(x - 1) \sqrt{2(x - \ln(x) - 1)}, \quad (51)$$

and $[g(D_{\ell})]_{ij} = g(D_{\ell,ii})\delta_{ij}$.

In order to transform Equation (50b) in a form that is quadratic also in the matrix elements, we exploit the following matrix identity<fnref n="6" />:

$$\frac{2\ell + 1}{2} \operatorname{Trace}[(C_{f\ell}^{-1/2} C_{g\ell} C_{f\ell}^{-1/2})^2] = X_{g\ell}^T M_{f\ell}^{-1} X_{g\ell}, \quad (53)$$

$$\operatorname{Trace}[g(D_{\ell})] = \operatorname{Trace}[U g(D_{\ell}) U^{T}] \quad (52a)$$

$$= \operatorname{Trace}[C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}] \quad (52b)$$

$$= \operatorname{Trace}[C_{f\ell}^{-1/2} C_{f\ell}^{1/2} (C_{\ell}^{-1/2} \hat{C}_{\ell} C_{\ell}^{-1/2}) C_{f\ell}^{1/2} C_{f\ell}^{-1/2}] \quad (52c)$$

$$= \operatorname{Trace}[C_{f\ell}^{-1/2} C_{g\ell} C_{f\ell}^{-1/2}] \quad (52d)$$

<footnote n="5">⁵In passing from Equations (49a) to (49b), we have made use of the properties of the trace[] operator.</footnote>

<footnote n="6">⁶Using the invariance of $\operatorname{Tr}[A]$ under diagonalization of $A$, one has that</footnote>
where $X_{g\ell}$ is the vector of $C_{g\ell} \equiv C_{f\ell}^{1/2} U_{\ell} g(D_{\ell}) U_{\ell}^{T} C_{f\ell}^{1/2}$ elements for a given fiducial model $C_{f\ell}$, with dimension $n(n+1)/2$, and $M_{f\ell}$ is the covariance of $\hat{X}$ evaluated at $C_{f\ell} = C_{f\ell}$. Therefore, Equation (50b) can be rewritten as

$$ -2 \ln \mathcal{L}(C_{\ell}) = \frac{2\ell + 1}{2} \mathrm{Trace}[(C_{f\ell}^{-1/2} C_{g\ell} C_{f\ell}^{-1/2})^2] \quad (54a) \\ = X_{g\ell}^T M_{f\ell}^{-1} X_{g\ell}. \quad (54b) $$

We stress that this formulation is exact in the full sky regime, as it has been obtained by means of matrix identities and no approximations have been adopted so far.

### 3.3. Including the Effects of Noise and Beam Smearing

The signal observed with a real CMB experiment is affected by various sources of experimental contaminations. Here, we focus on two main classes of experimental effects: the noise bias and the beam smearing.

The noise bias is due to the instrumental noise from detectors that adds up to the cosmological signal. It has to be characterized and subtracted from the (overall) signal (an alternative approach posits in cross-correlating different detectors and assuming their individual noise to be uncorrelated, see section 3.4). In the simple case of isotropic noise in real space, the noise level is independent from the direction. This translates in a diagonal noise in harmonic space, i.e., the noise power spectrum $N_{\ell}$ is an additive bias for the CMB power spectrum. A very simple example is the case of isotropic and homogenous noise in real space, i.e., the noise level is the same in each pixel. This translates in a “white noise” in harmonic space, i.e., $N_{\ell}$ is constant in $\ell$. A usual assumption is also to consider the noise in temperature and polarization to be uncorrelated. If the noise is anisotropic (i.e., it changes from pixel to pixel) for example because of a particular scanning strategy that induces anisotropic sky coverage, it may induce correlations between $a_{\ell m}$, and different considerations apply.

The beam smearing is due to the fact that the instrument has a finite angular response. The signal observed along a certain direction takes contributions from all angular directions. These contributions are weighted with the angular response of the instrument. In real space, this effect is described as a convolution of the observed signal with the angular response (hereafter beam) of the instrument $\Theta(\theta, \phi) \propto \int d\Omega' B(\theta - \theta', \phi - \phi')\Theta(\theta', \phi')$. In harmonic space, the convolution becomes a product between the harmonic expansions of the CMB fields and the beam $a'_{\ell m} = b_{\ell m} a_{\ell m}$. In the simple case of gaussian beam of width $\sigma_{FWHM} = 2\sqrt{2\ln 2\sigma}$, the harmonic expression of the beam is independent from $m$ and takes the simple form of $b_{\ell m} \rightarrow B_{\ell} \equiv \exp[-\ell(\ell + 1)\sigma^2] = \exp[-\ell(\ell + 1)\sigma_{FWHM}^2/(8\ln 2)]$. The beam smearing is a multiplicative bias for the CMB power spectrum.

In presence of noise and beam smearing, the observed signal is $d_{\ell m} = B_{\ell} a_{\ell m} + n_{\ell}$, and the estimator in Equation (9) becomes

$$ \hat{C}_{\ell} \rightarrow \hat{D}_{\ell} \equiv \frac{1}{2\ell + 1} \sum_{m} d_{\ell m} d_{\ell m}^{*} = B_{\ell}^{2} \hat{C}_{\ell} + \hat{N}_{\ell} \quad (55) $$

From Equation (55) it is clear that $\hat{D}_{\ell}$ is a biased estimator of the true power spectrum $C_{\ell}$, $\langle \hat{D}_{\ell} \rangle = B_{\ell}^2 C_{\ell} + N_{\ell}$. In addition, the variance of the estimator takes an additional contribution. In presence of noise and beam effects, the variance becomes

$$ \mathrm{var}[\hat{C}_{\ell}] \equiv \langle (C_{\ell} - \hat{C}_{\ell})(C_{\ell} - \hat{C}_{\ell}) \rangle = \frac{2}{2\ell + 1} \left( C_{\ell} + \frac{N_{\ell}}{B_{\ell}^2} \right)^2 \quad (56) $$

The variable $\hat{D}_{\ell} \equiv \hat{C}_{\ell} + \hat{N}_{\ell}/B_{\ell}^2$ still has a $\Gamma$ distribution, and all the considerations made for the $\hat{C}_{\ell}$ estimator still apply to the (slightly) more general case of noise and beam biases, provided that $\hat{C}_{\ell}$ is replaced with $\hat{D}_{\ell}$. More in detail, when both temperature and polarization are considered, the matrix of estimators $S_{\ell}$ still has a Wishart distribution [see Equation (17)] with a revised $W_{\ell} = V_{\ell}/(2\ell + 1)$ matrix$^7$

$$ W_a = \frac{1}{2\ell + 1} \begin{pmatrix} C_{\ell}^{TT} + N_{\ell}^T / (B_{\ell}^T)^2 & C_{\ell}^{TE} & 0 \\ C_{\ell}^{TE} & C_{\ell}^{EE} + N_{\ell}^P / (B_{\ell}^P)^2 & 0 \\ 0 & 0 & C_{\ell}^{BB} + N_{\ell}^P / (B_{\ell}^P)^2 \end{pmatrix}, \quad (57) $$

where we have allowed for different (and uncorrelated) noise levels in temperature and polarization, and for different beam sizes in temperature and polarization.

The covariance matrix of the estimators $\hat{D}_{\ell}$ is equivalent to that of $\hat{C}_{\ell}$ in Equation (20a), provided that $C_{\ell}^{TT}$, $C_{\ell}^{EE}$ are replaced with the power spectra corrected for noise and beam.

### 3.4. Multi-frequency Analysis

The results presented so far have been discussed assuming the special viewpoint of a single-frequency experiment. In reality, CMB experiments often rely on multi-frequency observations to better characterize the cosmological signal and extract it from the multi-component sky-signal observed (see discussions in e.g., [37, 40]). Moreover, multiple detectors sharing the same central frequency are always available, so that the final signal at a certain frequency can be effectively thought as a weighted average of the signals observed with multiple detectors.

In general, if $n$ maps are available, there are $n-1$ combinations that are independent from the signal (if two maps share the same signal and have different noise properties, their difference is independent from the common signal). There is one independent combination defined as the weighted average of the $n$ available maps

$$ a_{\ell m}^{NW} = \sum_{i}^{n} w_i a_{\ell m}^{i} \quad (58) $$

where $w_i$ are the noise weights associated to each of the $n$ maps. In the simple scenario of isotropic noise $N_{\ell}$ for each map, the weights can be defined as $w_i = (N_{\ell}^i)^{-1} / \sum_i (N_{\ell}^i)^{-1}$. The noise-weighted map is a sufficient statistics for the CMB field, and therefore all the considerations above about the

<footnote>$^7$We are dropping the BB part of the distribution, as we have seen that the full Wishart in full-sky is separable into a $T-E$ and $B$ component. For the $B$ component, all considerations in the single-field regime apply.</footnote>
where $\mathbf{X}_{g\ell}$ is the vector of $C_{g\ell} \equiv \mathbf{C}_{f\ell}^{1/2} \mathbf{U}_{\ell} \mathbf{g}(\mathbf{D}_{\ell}) \mathbf{U}_{\ell}^{T} \mathbf{C}_{f\ell}^{1/2}$ elements for a given fiducial model $\mathbf{C}_{f\ell}$, with dimension $n(n+1)/2$, and $\mathbf{M}_{f\ell}$ is the covariance of $\hat{\mathbf{X}}$ evaluated at $\mathbf{C}_{\ell} = \mathbf{C}_{f\ell}$. Therefore, Equation (50b) can be rewritten as

$$ -2 \ln \mathcal{L}(\mathbf{C}_{\ell}) = \frac{2\ell + 1}{2} \mathrm{Trace}[(\mathbf{C}_{f\ell}^{-1/2} \mathbf{C}_{g\ell} \mathbf{C}_{f\ell}^{-1/2})^2] \quad (54a) \\ = \mathbf{X}_{g\ell}^{T} \mathbf{M}_{f\ell}^{-1} \mathbf{X}_{g\ell}. \quad (54b) $$

We stress that this formulation is exact in the full sky regime, as it has been obtained by means of matrix identities and no approximations have been adopted so far.

### 3.3. Including the Effects of Noise and Beam Smearing

The signal observed with a real CMB experiment is affected by various sources of experimental contaminations. Here, we focus on two main classes of experimental effects: the noise bias and the beam smearing.

The noise bias is due to the instrumental noise from detectors that adds up to the cosmological signal. It has to be characterized and subtracted from the (overall) signal (an alternative approach posits in cross-correlating different detectors and assuming their individual noise to be uncorrelated, see section 3.4). In the simple case of isotropic noise in real space, the noise level is independent from the direction. This translates in a diagonal noise in harmonic space, i.e., the noise power spectrum $N_{\ell}$ is an additive bias for the CMB power spectrum. A very simple example is the case of isotropic and homogenous noise in real space, i.e., the noise level is the same in each pixel. This translates in a “white noise” in harmonic space, i.e., $N_{\ell}$ is constant in $\ell$. A usual assumption is also to consider the noise in temperature and polarization to be uncorrelated. If the noise is anisotropic (i.e., it changes from pixel to pixel) for example because of a particular scanning strategy that induces anisotropic sky coverage, it may induce correlations between $a_{\ell m}$, and different considerations apply.

The beam smearing is due to the fact that the instrument has a finite angular response. The signal observed along a certain direction takes contributions from all angular directions. These contributions are weighted with the angular response of the instrument. In real space, this effect is described as a convolution of the observed signal with the angular response (hereafter beam) of the instrument $\Theta(\theta, \phi) \propto \int d\Omega' B(\theta - \theta', \phi - \phi')\Theta(\theta', \phi')$. In harmonic space, the convolution becomes a product between the harmonic expansions of the CMB fields and the beam $a'_{\ell m} = b_{\ell m} a_{\ell m}$. In the simple case of gaussian beam of width $\sigma_{\text{FWHM}} = 2\sqrt{2\ln 2}\sigma$, the harmonic expression of the beam is independent from $m$ and takes the simple form of $b_{\ell m} \rightarrow B_{\ell} \equiv \exp[-\ell(\ell+1)\sigma^2] = \exp[-\ell(\ell+1)\sigma_{\text{FWHM}}^2/(8\ln 2)]$. The beam smearing is a multiplicative bias for the CMB power spectrum.

In presence of noise and beam smearing, the observed signal is $d_{\ell m} = B_{\ell} a_{\ell m} + n_{\ell}$, and the estimator in Equation (9) becomes

$$ \hat{C}_{\ell} \rightarrow \hat{D}_{\ell} \equiv \frac{1}{2\ell + 1} \sum_{m} d_{\ell m} d_{\ell m}^{*} = B_{\ell}^{2} \hat{C}_{\ell} + \hat{N}_{\ell} \quad (55) $$

From Equation (55) it is clear that $\hat{D}_{\ell}$ is a biased estimator of the true power spectrum $C_{\ell}$, $\langle \hat{D}_{\ell} \rangle = B_{\ell}^2 C_{\ell} + N_{\ell}$. In addition, the variance of the estimator takes an additional contribution. In presence of noise and beam effects, the variance becomes

$$ \mathrm{var}[\hat{C}_{\ell}] \equiv \langle (C_{\ell} - \hat{C}_{\ell})(C_{\ell} - \hat{C}_{\ell}) \rangle = \frac{2}{2\ell + 1} \left( C_{\ell} + \frac{N_{\ell}}{B_{\ell}^2} \right)^2 \quad (56) $$

The variable $\hat{D}_{\ell} \equiv \hat{C}_{\ell} + \hat{N}_{\ell}/B_{\ell}^2$ still has a $\Gamma$ distribution, and all the considerations made for the $\hat{C}_{\ell}$ estimator still apply to the (slightly) more general case of noise and beam biases, provided that $\hat{C}_{\ell}$ is replaced with $\hat{D}_{\ell}$. More in detail, when both temperature and polarization are considered, the matrix of estimators $S_{\ell}$ still has a Wishart distribution [see Equation (17)] with a revised $W_{\ell} = V_{\ell}/(2\ell + 1)$ matrix$^7$

$$ W_a = \frac{1}{2\ell + 1} \begin{pmatrix} C_{\ell}^{TT} + N_{\ell}^T / (B_{\ell}^T)^2 & C_{\ell}^{TE} & 0 \\ C_{\ell}^{TE} & C_{\ell}^{EE} + N_{\ell}^P / (B_{\ell}^P)^2 & 0 \\ 0 & 0 & C_{\ell}^{BB} + N_{\ell}^P / (B_{\ell}^P)^2 \end{pmatrix}, \quad (57) $$

where we have allowed for different (and uncorrelated) noise levels in temperature and polarization, and for different beam sizes in temperature and polarization.

The covariance matrix of the estimators $\hat{D}_{\ell}$ is equivalent to that of $\hat{C}_{\ell}$ in Equation (20a), provided that $C_{\ell}^{TT}$, $C_{\ell}^{EE}$ are replaced with the power spectra corrected for noise and beam.

### 3.4. Multi-frequency Analysis

The results presented so far have been discussed assuming the special viewpoint of a single-frequency experiment. In reality, CMB experiments often rely on multi-frequency observations to better characterize the cosmological signal and extract it from the multi-component sky-signal observed (see discussions in e.g., [37, 40]). Moreover, multiple detectors sharing the same central frequency are always available, so that the final signal at a certain frequency can be effectively thought as a weighted average of the signals observed with multiple detectors.

In general, if $n$ maps are available, there are $n-1$ combinations that are independent from the signal (if two maps share the same signal and have different noise properties, their difference is independent from the common signal). There is one independent combination defined as the weighted average of the $n$ available maps

$$ a_{\ell m}^{\text{NW}} = \sum_{i}^{n} w_i a_{\ell m}^{i} \quad (58) $$

where $w_i$ are the noise weights associated to each of the $n$ maps. In the simple scenario of isotropic noise $N_{\ell}$ for each map, the weights can be defined as $w_i = (N_{\ell}^i)^{-1} / \sum_i (N_{\ell}^i)^{-1}$. The noise-weighted map is a sufficient statistics for the CMB field, and therefore all the considerations above about the

<footnote>$^7$We are dropping the BB part of the distribution, as we have seen that the full Wishart in full-sky is separable into a $T-E$ and $B$ component. For the $B$ component, all considerations in the single-field regime apply.</footnote>
choice of the likelihood also apply to the combined map. An estimator for the power spectrum can be constructed from the noise-weighted map.

Another possibility is to build estimators $\hat{C}_{\ell}^{ij} = (1/(2\ell + 1)) \sum_m a_{\ell m}^i (a_{\ell m}^{j})^*$ from pairs of maps and then define a weighted estimator $\hat{C}_{\ell}^{NW} = \sum_{ij} w_{ij} \hat{C}_{\ell}^{ij}$, where the weights $w_{ij}$ depend on the noise levels in the individual maps. It can be shown that the latter solution is equivalent to estimating the observed power spectrum from the noise-weighted map, and again all the considerations about the likelihood choice apply to this case as well (see discussion in Hamimeche and Lewis [32], appendix C).

Regarding the latter solution, a more robust choice is to build the noise-weighted estimator $\hat{C}_{\ell}$ from cross-spectra only, i.e., from pairs $(ij)$ with $i \neq j$. If the noise in the individual maps is uncorrelated, to take cross-spectra is safe with respect to the introduction of possible biases in the final estimator due to unaccounted errors in the noise model<fnref n="8" />. However, the statistics of the estimator obtained from cross-spectra $\hat{C}_{\ell}^{CS}$ may differ from that of the generic noise-weighted estimator. In particular, the cross-spectra might not be positive-definite. Therefore, in principle one should use a distribution for $\hat{C}_{\ell}^{CS}$ other than the Wishart. However, it can be demonstrated (see e.g., appendix C in Hamimeche and Lewis [32]) that, in the limit of many maps available, the distribution of $\hat{C}_{\ell}^{CS}$ approaches that of $\hat{C}_{\ell}^{NW}$, and hence one can use the same approximations developed in the case of the generic noise-weighted estimator.

Before concluding this subsection, a note on the covariance matrix. When multiple maps are available and the estimators are build from a combination of those maps, the expression for the covariance matrix can be generalized as follows:

$$ \mathrm{Cov}[ij C_{\ell_1}^{XY} \quad ab C_{\ell_2}^{WZ}] = \langle ij C_{\ell_1}^{XY} \times ab C_{\ell_2}^{WZ} \rangle - \langle ij C_{\ell_1}^{XY} \rangle \langle ab C_{\ell_2}^{WZ} \rangle \quad (59) $$

where $ij, ab$ denote all possible combinations of pairs of maps, XY, WZ are pairs of fields $T, E, B$, and we have explicitly taken into account the possibility of mode-coupling between $\ell_1, \ell_2$ (e.g., in the cut-sky regime). In the simple case of single-map in full-sky, Equation (59) reduces to Equation (20a).

## 4. GRAVITATIONAL LENSING

In the discussions so far, we have implicitly assumed that the CMB fields are unlensed. This is not the case in the reality. CMB photons traveling from the last scattering surface to the observer feel the gravitational effects of the evolving structures in the Universe. This effect is analogous to the weak lensing effect observed in galaxy surveys, where images of source galaxies are distorted and magnified by foreground structures acting as gravitational lenses. In the CMB case, the CMB as emitted at the last scattering surface is the source and the whole distribution of total (cold and baryonic) matter along the line of sight acts as the foreground lens. In practice, this means that the anisotropy

fields observed at a certain direction in the sky are displaced with respect to the original direction of emission:

$$ X(\hat{n}) = X^{\mathrm{unl}}(\hat{n} + \nabla\phi(\hat{n})) \quad (60) $$

where $X = T, E, B$ and $\phi$ is the lensing potential. The gradient of the lensing potential gives the deflection angle $\alpha = \nabla\phi$. The typical deflection that CMB anisotropies undergo is of order 2.5 arcmin [41]. The lensing potential is given by the integrated contribution of the gravitational potential along the line of sight<fnref n="9" />:

$$ \phi(\hat{n}) = -2 \int_{0}^{\chi_{*}} d\chi W(\chi) \Psi(\hat{n}) \quad (61) $$

where $\Psi$ is the (Weyl) gravitational potential, $\chi$ is the conformal distance and $W(\chi)$ is a geometrical kernel.

Gravitational lensing preserves the total variance of the CMB fields, being a bare displacement of the anisotropy distribution. Very roughly speaking, if we extracted the CMB power spectra from small patches of the sky<fnref n="10" />, we would observe the acoustic peaks to shift to either smaller or larger scales with respect to the full sky average (see discussion in e.g., [42]). The net effect is a smoothing of the acoustic peak structure in the CMB power spectra that can be as high as 20% in the case of the sharper structure in the EE power spectrum with respect to the unlensed case. Another important effect is the generation of spurious (i.e., not primordial) B-modes from the lensing of primordial E-modes, with a power spectrum that resembles a white noise contribution with noise level of ~ 5 µK arcmin at $\ell < 100$, representing a serious contaminant for searches of primordial B-modes.

A detailed description of the effects of gravitational lensing on the CMB spectra is beyond the scope of this manuscript, and can be found in the excellent review by Lewis and Challinor [41]. Here, it is relevant to stress that gravitational lensing modifies the statistical properties of the primary CMB fields in two ways. First, let's consider a fixed lensing potential realization and ensemble average over the CMB realizations. If we Taylor- expand Equation (60), take the harmonic expansion coefficients and compute the covariance of two fields $X, X' = T, E, B$, we get [43, 44]

$$ \langle x_{\ell m} x'_{\ell' m'} \rangle |_{\mathrm{lens}} = C_{\ell}^{\ell\ell'} \delta_{\ell\ell'} \delta_{m-m'} (-1)^m + \sum_{LM} (-1)^M \begin{pmatrix} \ell & \ell' & L \\ m & m' & -M \end{pmatrix} f_{\ell L \ell'}^{\alpha} \phi_{LM} \quad (62) $$

where the term in brackets is the Wigner-3j and $f^{\alpha}$ is a weight of different $\ell\ell'$ pairs depending on the unlensed power spectra and on geometrical factors (a full list can be found in Okamoto and Hu [43]). The $C_{\ell}$ in the first term of the RHS are the lensed power spectra. From Equation (62), it is clear that gravitational lensing not only modifies the structure of the primary CMB spectra

<footnote n="8">⁸If $a_{\ell m}^i = s_{\ell m}^i + n_{\ell}^i$, with noise $n^i$ uncorrelated for any map $i = 1, ..., nmaps$ and signal $s^i$, then $\sum_m |a_{\ell m}^i (a_{\ell m}^j)^*| = \sum_m |s_{\ell m}^i (s_{\ell m}^j)^* + n_{\ell m}^i (n_{\ell m}^j)^*| = (2\ell + 1)\hat{C}_{\ell}^{ij}$.</footnote>

<footnote n="9">⁹The unperturbed line of sight in the Born approximation.</footnote>

<footnote n="10">¹⁰Small enough that the convergence and shear components of the deflection field can be assumed constant in the patch.</footnote>
choice of the likelihood also apply to the combined map. An estimator for the power spectrum can be constructed from the noise-weighted map.

Another possibility is to build estimators $\hat{C}_{\ell}^{ij} = (1/(2\ell + 1)) \sum_m a_{\ell m}^i (a_{\ell m}^{j})^*$ from pairs of maps and then define a weighted estimator $\hat{C}_{\ell}^{NW} = \sum_{ij} w_{ij} \hat{C}_{\ell}^{ij}$, where the weights $w_{ij}$ depend on the noise levels in the individual maps. It can be shown that the latter solution is equivalent to estimating the observed power spectrum from the noise-weighted map, and again all the considerations about the likelihood choice apply to this case as well (see discussion in Hamimeche and Lewis [32], appendix C).

Regarding the latter solution, a more robust choice is to build the noise-weighted estimator $\hat{C}_{\ell}$ from cross-spectra only, i.e., from pairs $(ij)$ with $i \neq j$. If the noise in the individual maps is uncorrelated, to take cross-spectra is safe with respect to the introduction of possible biases in the final estimator due to unaccounted errors in the noise model<fnref n="8" />. However, the statistics of the estimator obtained from cross-spectra $\hat{C}_{\ell}^{CS}$ may differ from that of the generic noise-weighted estimator. In particular, the cross-spectra might not be positive-definite. Therefore, in principle one should use a distribution for $\hat{C}_{\ell}^{CS}$ other than the Wishart. However, it can be demonstrated (see e.g., appendix C in Hamimeche and Lewis [32]) that, in the limit of many maps available, the distribution of $\hat{C}_{\ell}^{CS}$ approaches that of $\hat{C}_{\ell}^{NW}$, and hence one can use the same approximations developed in the case of the generic noise-weighted estimator.

Before concluding this subsection, a note on the covariance matrix. When multiple maps are available and the estimators are build from a combination of those maps, the expression for the covariance matrix can be generalized as follows:

$$ \mathrm{Cov}[ij C_{\ell_1}^{XY} \quad ab C_{\ell_2}^{WZ}] = \langle ij C_{\ell_1}^{XY} \times ab C_{\ell_2}^{WZ} \rangle - \langle ij C_{\ell_1}^{XY} \rangle \langle ab C_{\ell_2}^{WZ} \rangle \quad (59) $$

where $ij, ab$ denote all possible combinations of pairs of maps, XY, WZ are pairs of fields $T, E, B$, and we have explicitly taken into account the possibility of mode-coupling between $\ell_1, \ell_2$ (e.g., in the cut-sky regime). In the simple case of single-map in full-sky, Equation (59) reduces to Equation (20a).

## 4. GRAVITATIONAL LENSING

In the discussions so far, we have implicitly assumed that the CMB fields are unlensed. This is not the case in the reality. CMB photons traveling from the last scattering surface to the observer feel the gravitational effects of the evolving structures in the Universe. This effect is analogous to the weak lensing effect observed in galaxy surveys, where images of source galaxies are distorted and magnified by foreground structures acting as gravitational lenses. In the CMB case, the CMB as emitted at the last scattering surface is the source and the whole distribution of total (cold and baryonic) matter along the line of sight acts as the foreground lens. In practice, this means that the anisotropy

fields observed at a certain direction in the sky are displaced with respect to the original direction of emission:

$$ X(\hat{n}) = X^{\mathrm{unl}}(\hat{n} + \nabla\phi(\hat{n})) \quad (60) $$

where $X = T, E, B$ and $\phi$ is the lensing potential. The gradient of the lensing potential gives the deflection angle $\alpha = \nabla\phi$. The typical deflection that CMB anisotropies undergo is of order 2.5 arcmin [41]. The lensing potential is given by the integrated contribution of the gravitational potential along the line of sight<fnref n="9" />:

$$ \phi(\hat{n}) = -2 \int_{0}^{\chi_{*}} d\chi W(\chi) \Psi(\hat{n}) \quad (61) $$

where $\Psi$ is the (Weyl) gravitational potential, $\chi$ is the conformal distance and $W(\chi)$ is a geometrical kernel.

Gravitational lensing preserves the total variance of the CMB fields, being a bare displacement of the anisotropy distribution. Very roughly speaking, if we extracted the CMB power spectra from small patches of the sky<fnref n="10" />, we would observe the acoustic peaks to shift to either smaller or larger scales with respect to the full sky average (see discussion in e.g., [42]). The net effect is a smoothing of the acoustic peak structure in the CMB power spectra that can be as high as 20% in the case of the sharper structure in the EE power spectrum with respect to the unlensed case. Another important effect is the generation of spurious (i.e., not primordial) B-modes from the lensing of primordial E-modes, with a power spectrum that resembles a white noise contribution with noise level of ~ 5 µK arcmin at $\ell < 100$, representing a serious contaminant for searches of primordial B-modes.

A detailed description of the effects of gravitational lensing on the CMB spectra is beyond the scope of this manuscript, and can be found in the excellent review by Lewis and Challinor [41]. Here, it is relevant to stress that gravitational lensing modifies the statistical properties of the primary CMB fields in two ways. First, let's consider a fixed lensing potential realization and ensemble average over the CMB realizations. If we Taylor- expand Equation (60), take the harmonic expansion coefficients and compute the covariance of two fields $X, X' = T, E, B$, we get [43, 44]

$$ \langle x_{\ell m} x'_{\ell' m'} \rangle |_{\mathrm{lens}} = C_{\ell}^{\ell\ell'} \delta_{\ell\ell'} \delta_{m-m'} (-1)^m + \sum_{LM} (-1)^M \begin{pmatrix} \ell & \ell' & L \\ m & m' & -M \end{pmatrix} f_{\ell L \ell'}^{\alpha} \phi_{LM} \quad (62) $$

where the term in brackets is the Wigner-3j and $f^{\alpha}$ is a weight of different $\ell\ell'$ pairs depending on the unlensed power spectra and on geometrical factors (a full list can be found in Okamoto and Hu [43]). The $C_{\ell}$ in the first term of the RHS are the lensed power spectra. From Equation (62), it is clear that gravitational lensing not only modifies the structure of the primary CMB spectra

<footnote n="8">⁸If $a_{\ell m}^i = s_{\ell m}^i + n_{\ell}^i$, with noise $n^i$ uncorrelated for any map $i = 1, ..., n_{\mathrm{maps}}$ and signal $s^i$, then $\sum_m |a_{\ell m}^i (a_{\ell m}^j)^*| = \sum_m |s_{\ell m}^i (s_{\ell m}^j)^* + n_{\ell m}^i (n_{\ell m}^j)^*| = (2\ell + 1) \hat{C}_{\ell}^{ij}$.</footnote>

<footnote n="9">⁹The unperturbed line of sight in the Born approximation.</footnote>

<footnote n="10">¹⁰Small enough that the convergence and shear components of the deflection field can be assumed constant in the patch.</footnote>
by smearing the acoustic peaks, but also induces off-diagonal covariance terms (the second term in the RHS). Therefore, for a fixed lensing realization, the CMB field becomes anisotropic. This property can be exploited to construct an estimator for the lensing potential $\phi$.

The lensing power spectrum $C_L^{\phi\phi} = (2L+1)^{-1} \sum_M \langle\phi_{LM}\phi_{LM}^*\rangle$ (where we are now taking the ensemble average over the lensing realizations as well) is related to the 4-point correlation function of the primary CMB fields, as it is clear by inspecting Equation (62). In other words, the second modification to the CMB statistics induced by lensing, when the ensemble average is taken over both the CMB and the lensing realizations, is a certain amount of non-Gaussianity measurable from the 4-point function. This property is exploited to construct an estimator for the lensing power spectrum.

The presence of gravitational lensing effects represent a pernicious contaminant for searches for primordial GW. In this case, “delensing” procedures (see e.g., [41] and references therein for a description of delensing procedures) aimed to remove the lensing contamination from the measured CMB sky are crucial to allow for the possibility to detect the primordial tensor BB signal. On the other hand, the presence of gravitational lensing also enriches the amount of information that we can extract from the observations of the CMB sky. As an example, since the gravitational lensing is induced by forming structures, it carries information about the late time evolution of the Universe and can be exploited to constrain those cosmological parameters that govern those stages of the Universe expansion, such as massive neutrinos.

There are two ways in which the information encoded in the gravitational lensing signal can be accessed from CMB measurements. One can exploit the anisotropy and non-Gaussianity properties of the lensed CMB fields to construct estimators for the gravitational potential and for the lensing power spectrum, and use those observables directly in a likelihood analysis. This is what is done e.g., by the Planck collaboration [45], ACT [46], SPT [47], POLARBEAR [48], and this is the goal of future CMB experiments that will be able to reconstruct the lensing signal over a large range of angular scales with exquisite sensitivity (CORE [49], Simons Observatory [20], CMB-S4 [19], PICO [23]). Concerning the choice of the likelihood function for the lensing signal itself, it has been shown [50] that a quadratic expression in the observed lensing power spectrum with a non-diagonal fiducial covariance matrix provides satisfying results.

As seen in Equation (62), cosmological information carried by the gravitational lensing signal are also encoded in the CMB power spectra. In fact, the lensing contribution modifies the primary acoustic structure as emerged from last scattering. Current generation of Boltzmann solvers usually employed in cosmological analysis, such as CAMB [51] and CLASS [52], carefully compute the lensing effect when constructing the theoretical spectra to be compared with the measured spectra in the likelihood analysis. In general, when lensed spectra are considered, the likelihood analysis should reflect the different statistical properties of the lensed CMB. First of all, the non-Gaussian distribution of the lensed field may require to

build a different likelihood function for the lensed fields. Even if the Gaussian approximation can be retained given the current experimental sensitivity, the covariance matrix may still need to reflect the presence of off-diagonal correlations induced by the lensing signal [53]. It has been shown that ignoring such correlations can still be safe for the sensitivity level of Planck [53, 54], but it will become relevant for the next generation of experiments [55, 56].

A final remark concerns cosmological analysis employing the combination of CMB power spectra and the lensing power spectrum. These two data sets are usually treated as independent. However, they are extracted from the same map and, as such, there is a certain level of correlation between the two. As pointed out in Schmittfull et al. [50], there are two sources of correlations: cosmic variance in the CMB field, which may affect the lensing reconstruction; cosmic variance in the lensing field, which affects not only the lensing reconstruction, but can propagate to the lensed CMB power spectra. An alternative solution that might reduce the correlation between fields is the joint analysis of the lensing power spectrum and unlensed CMB spectra. The reconstructed lensing signal—either from CMB observations themselves or from tracers of large-scale structure correlated with the lensing signal such as the cosmic infrared background (CIB)—can be used to delens the primary CMB spectra. This procedure is not only key to unveil the primordial BB tensor spectrum, but it can also improve the sensitivity to those cosmological parameters that are mostly constrained via measurements of the high-$\ell$ damping tail in the *TT*, *TE*, *EE* spectra, such as $N_{\text{eff}}$ [19, 57].

All the above concerns will become much more significant for the next generation of experiments, when more sensitive polarization-based reconstructions will be available.

# 5. LIKELIHOOD APPROACHES—LARGE-SCALE REGIME

In this section, we review the basics of the likelihood approaches at large scales (low multipoles). We consider two generic classes of likelihood methods: pixel-based and simulation-based. When the focus of the likelihood analysis is the sky at large angular scales, the resolution of the map to be analyzed is low enough to make a pixel-based approach computationally feasible. Alternative approaches exploit the information encoded in harmonic space, and build the likelihood function from a simulation-based method or component-separation based method (Blackwell-Rao). For large scale studies the Hamimeche-Lewis likelihood can be also considered, [58, 59], for completeness we will explore later the performance of such likelihood compared with other approximations.

## 5.1. Pixel-Based Approach

The great advantage of the pixel-based approach lies in the fact that the likelihood function so defined is always exact (see section 2.3), including in the cut-sky regime. Equation (21) can be easily generalized in the presence of (Gaussian-distributed) noise, by defining the data vector as $\mathbf{m} = \mathbf{s} + \mathbf{n}$, where $\mathbf{s}$ is the signal per
by smearing the acoustic peaks, but also induces off-diagonal covariance terms (the second term in the RHS). Therefore, for a fixed lensing realization, the CMB field becomes anisotropic. This property can be exploited to construct an estimator for the lensing potential $\phi$.

The lensing power spectrum $C_L^{\phi\phi} = (2L+1)^{-1} \sum_M \langle\phi_{LM}\phi_{LM}^*\rangle$ (where we are now taking the ensemble average over the lensing realizations as well) is related to the 4-point correlation function of the primary CMB fields, as it is clear by inspecting Equation (62). In other words, the second modification to the CMB statistics induced by lensing, when the ensemble average is taken over both the CMB and the lensing realizations, is a certain amount of non-Gaussianity measurable from the 4-point function. This property is exploited to construct an estimator for the lensing power spectrum.

The presence of gravitational lensing effects represent a pernicious contaminant for searches for primordial GW. In this case, “delensing” procedures (see e.g., [41] and references therein for a description of delensing procedures) aimed to remove the lensing contamination from the measured CMB sky are crucial to allow for the possibility to detect the primordial tensor BB signal. On the other hand, the presence of gravitational lensing also enriches the amount of information that we can extract from the observations of the CMB sky. As an example, since the gravitational lensing is induced by forming structures, it carries information about the late time evolution of the Universe and can be exploited to constrain those cosmological parameters that govern those stages of the Universe expansion, such as massive neutrinos.

There are two ways in which the information encoded in the gravitational lensing signal can be accessed from CMB measurements. One can exploit the anisotropy and non-Gaussianity properties of the lensed CMB fields to construct estimators for the gravitational potential and for the lensing power spectrum, and use those observables directly in a likelihood analysis. This is what is done e.g., by the Planck collaboration [45], ACT [46], SPT [47], POLARBEAR [48], and this is the goal of future CMB experiments that will be able to reconstruct the lensing signal over a large range of angular scales with exquisite sensitivity (CORE [49], Simons Observatory [20], CMB-S4 [19], PICO [23]). Concerning the choice of the likelihood function for the lensing signal itself, it has been shown [50] that a quadratic expression in the observed lensing power spectrum with a non-diagonal fiducial covariance matrix provides satisfying results.

As seen in Equation (62), cosmological information carried by the gravitational lensing signal are also encoded in the CMB power spectra. In fact, the lensing contribution modifies the primary acoustic structure as emerged from last scattering. Current generation of Boltzmann solvers usually employed in cosmological analysis, such as CAMB [51] and CLASS [52], carefully compute the lensing effect when constructing the theoretical spectra to be compared with the measured spectra in the likelihood analysis. In general, when lensed spectra are considered, the likelihood analysis should reflect the different statistical properties of the lensed CMB. First of all, the non-Gaussian distribution of the lensed field may require to

build a different likelihood function for the lensed fields. Even if the Gaussian approximation can be retained given the current experimental sensitivity, the covariance matrix may still need to reflect the presence of off-diagonal correlations induced by the lensing signal [53]. It has been shown that ignoring such correlations can still be safe for the sensitivity level of Planck [53, 54], but it will become relevant for the next generation of experiments [55, 56].

A final remark concerns cosmological analysis employing the combination of CMB power spectra and the lensing power spectrum. These two data sets are usually treated as independent. However, they are extracted from the same map and, as such, there is a certain level of correlation between the two. As pointed out in Schmittfull et al. [50], there are two sources of correlations: cosmic variance in the CMB field, which may affect the lensing reconstruction; cosmic variance in the lensing field, which affects not only the lensing reconstruction, but can propagate to the lensed CMB power spectra. An alternative solution that might reduce the correlation between fields is the joint analysis of the lensing power spectrum and unlensed CMB spectra. The reconstructed lensing signal—either from CMB observations themselves or from tracers of large-scale structure correlated with the lensing signal such as the cosmic infrared background (CIB)—can be used to delens the primary CMB spectra. This procedure is not only key to unveil the primordial BB tensor spectrum, but it can also improve the sensitivity to those cosmological parameters that are mostly constrained via measurements of the high-$\ell$ damping tail in the *TT*, *TE*, *EE* spectra, such as $N_{\text{eff}}$ [19, 57].

All the above concerns will become much more significant for the next generation of experiments, when more sensitive polarization-based reconstructions will be available.

# 5. LIKELIHOOD APPROACHES—LARGE-SCALE REGIME

In this section, we review the basics of the likelihood approaches at large scales (low multipoles). We consider two generic classes of likelihood methods: pixel-based and simulation-based. When the focus of the likelihood analysis is the sky at large angular scales, the resolution of the map to be analyzed is low enough to make a pixel-based approach computationally feasible. Alternative approaches exploit the information encoded in harmonic space, and build the likelihood function from a simulation-based method or component-separation based method (Blackwell-Rao). For large scale studies the Hamimeche-Lewis likelihood can be also considered, [58, 59], for completeness we will explore later the performance of such likelihood compared with other approximations.

## 5.1. Pixel-Based Approach

The great advantage of the pixel-based approach lies in the fact that the likelihood function so defined is always exact (see section 2.3), including in the cut-sky regime. Equation (21) can be easily generalized in the presence of (Gaussian-distributed) noise, by defining the data vector as $\mathbf{m} = \mathbf{s} + \mathbf{n}$, where $\mathbf{s}$ is the signal per
pixel in temperature and polarization (s = (T, Q, U)) and n is the instrumental noise. We report here the likelihood function in pixel space for convenience:

$$ \mathcal{L}(C_{\ell}) = \mathcal{P}(\mathbf{m}|C_{\ell}) = \frac{1}{2\pi|M|^{1/2}} \exp\left(-\frac{1}{2}\mathbf{m}^{\mathrm{T}}\mathbf{M}^{-1}\mathbf{m}\right), \quad (63) $$

The full covariance matrix M in Equation (63) is consequently generalized to the sum of the signal and noise covariance matrices $M = S + N$. Furthermore, the effect of beam smearing discussed in section 3.3—and also relevant for the large-scale data—is now taken into account when constructing the full covariance matrix in terms of the beam-weighted sum of Legendre polynomial (see Equation (23)). For example, the explicit expression for the weight $P_{\ell}^{TT}$ for temperature becomes

$$ (P_{\ell}^{TT})_{i,j} = \frac{2\ell + 1}{4\pi} B_{\ell}^{2} P_{\ell}(\hat{\mathbf{n}}_{i} \cdot \hat{\mathbf{n}}_{j}) \quad (64) $$

where $\hat{\mathbf{n}}_i$ is the unit vector pointing toward the ith pixel, $B_\ell$ is given by the product of the instrumental beam Legendre transform and the (HEALPix [60]) pixel window $^{11}$, and $P_\ell$ is the Legendre polynomial of order $\ell$.

As already noted above, the actual feasibility of this mathematically rigorous approach only applies to the very large scales. Even so, massive parallel coding and memory requirements could still be a necessary ingredient. As an example, the Planck collaboration employed this approach in the low-$\ell$ likelihood analysis in temperature and polarization for the 2015 data release [37]. The map resolution was fixed at $N_{side} = 16$ (for comparison, the analysis conducted by the WMAP team employed $N_{side} = 8$ maps) to accommodate the $\ell < 47$ multipoles in the analysis, resulting in $N_{pix} = 3 \times 3072 = 9216$ total number of pixels, further reduced by the application of the analysis mask. In evaluating the likelihood function in Equation (63), the data vector and the noise covariance matrix are fixed, while the signal covariance matrix is recomputed for any given cosmological model to be compared against data, following Equation (23). In practice, only a subsection of S is recomputed, in particular that subsection corresponding to $\ell < 30$. The portion of S corresponding to multipoles $30 \le \ell < 47$ is precomputed from a fixed fiducial model. The choice of the fiducial model does not affect the performance of the likelihood. In fact, at the resolution employed in the large-scale analysis, the sensitivity to multipoles above $\ell \sim 30$ is strongly suppressed.

In 2013, a hybrid approach coupling a MonteCarlo-based approach in temperature (Blackwell-Rao estimator, see section 5.2) to a pixel-based approach in polarization for $\ell < 23$ was adopted [40]. In order to speed up the evaluation of the fully

pixel-based likelihood function for the 2015 release at any given cosmological model, the “brute-force” approach described here has been optimized with the implementation of the Sherman-Morrison-Woodbury formula, described below.

### 5.1.1. Sherman-Morrison-Woodbury Formula

The computational cost required by the pixel-based approach can be dramatically reduced if one consider that only a portion of the signal covariance matrix is reconstructed at any evaluation of the likelihood function. The full covariance matrix can then be decomposed into a varying part, which is function of the theoretical $C_\ell$ (i.e., the power spectra of the theoretical models to be compared against data), and a fixed part given by the fiducial S and the noise covariance matrix. The following step is to further decompose the varying part of the covariance as $S = V^TAV$, via a transformation V that effectively reduces the dimension of the actual evaluation cost from a $N_{pix} \times N_{pix}$ inversion to a $n_\lambda \times n_\lambda$ inversion, where $n_\lambda = 2\ell + 1$ is the dimension of the transformed matrix A. The latter is the only matrix that depends on the theoretical $C_\ell$ and, therefore, it is the only matrix to be recomputed and inverted. The fixed portions of the covariance matrix as well as the transformation matrix V can be pre-computed and stored. For the Planck 2015 data release, the application of such mathematical formalisms allowed to speed-up the likelihood evaluation by an order of magnitude. The mathematical details leading to the application of the Sherman-Morrison-Woodbury formula can be found in Aghanim et al. [37] (appendix B.1) and Hinshaw et al. [61] (appendix A.1).

## 5.2. Blackwell-Rao Estimator

An alternative approach to the likelihood evaluation that overcomes the computational cost of an exact likelihood evaluation in pixel space is provided by the combination of the Gibbs sampling method with the Blackwell-Rao estimator. The Gibbs sampling is a MCMC method applied to the estimation of the observed CMB signal from a raw map containing signal, noise, and foreground contamination in a Bayesian framework [62–64]. The crucial output for the subsequent construction of the likelihood estimator is a set of samples of the CMB sky, or more precisely, a set of sample variances of the sky samples obtained with the Gibbs algorithm. The Blackwell-Rao estimator is then built as an average over the set of sample variances.

Let's assume that the observed map m is composed by the CMB signal s and noise n (the following procedure can be generalized to the case in which foreground f also contribute to the total signal, see e.g., [65]). What we really want to evaluate is the joint probability of having a certain s with a certain $C_\ell$ given the observed sky m, i.e., we want to evaluate $\mathcal{P}(s, C_\ell|\mathbf{m})$. A brute force evaluation by computing a grid of s and $C_\ell$ is computationally prohibitive (it is even more prohibitive once one considers the inclusion of foreground into the equation). Another approach is to sample directly from the distribution via specific algorithms. Although it is again computationally unfeasible to sample directly from the joint distribution, as it would require inverting large and

<footnote>$^{11}$The signal in each pixel is the average over the signal at each point within the surface area of the pixel $\Omega_p$, $f_p = \int d\Omega w_p f(\Omega)$, where $w_p = 1/\Omega_p$ is the weight inside the pixel and zero otherwise. The harmonic transform of the weight $w_p$ is the pixel window function. Including the full dependence of the window function on m, $\ell$ can be computationally demanding. However, if the size of the pixel is small compared to the angular resolution of the experiment, the m-dependence can be neglected and the m-averaged window function $w_\ell$ can be defined. The power spectra computed from the pixelised maps $C_\ell^p$ are related to the unpixelised spectra $C_\ell$ via $C_\ell^p = w_\ell^2 C_\ell$.</footnote>
pixel in temperature and polarization (s = (T, Q, U)) and n is the instrumental noise. We report here the likelihood function in pixel space for convenience:

$$ \mathcal{L}(C_{\ell}) = \mathcal{P}(\mathbf{m}|C_{\ell}) = \frac{1}{2\pi|M|^{1/2}} \exp\left(-\frac{1}{2}\mathbf{m}^{\mathrm{T}}\mathbf{M}^{-1}\mathbf{m}\right), \quad (63) $$

The full covariance matrix M in Equation (63) is consequently generalized to the sum of the signal and noise covariance matrices $M = S + N$. Furthermore, the effect of beam smearing discussed in section 3.3—and also relevant for the large-scale data—is now taken into account when constructing the full covariance matrix in terms of the beam-weighted sum of Legendre polynomial (see Equation (23)). For example, the explicit expression for the weight $P_{\ell}^{TT}$ for temperature becomes

$$ (P_{\ell}^{TT})_{i,j} = \frac{2\ell + 1}{4\pi} B_{\ell}^{2} P_{\ell}(\hat{\mathbf{n}}_{i} \cdot \hat{\mathbf{n}}_{j}) \quad (64) $$

where $\hat{\mathbf{n}}_i$ is the unit vector pointing toward the ith pixel, $B_\ell$ is given by the product of the instrumental beam Legendre transform and the (HEALPix [60]) pixel window $^{11}$, and $P_\ell$ is the Legendre polynomial of order $\ell$.

As already noted above, the actual feasibility of this mathematically rigorous approach only applies to the very large scales. Even so, massive parallel coding and memory requirements could still be a necessary ingredient. As an example, the Planck collaboration employed this approach in the low-$\ell$ likelihood analysis in temperature and polarization for the 2015 data release [37]. The map resolution was fixed at $N_{\text{side}} = 16$ (for comparison, the analysis conducted by the WMAP team employed $N_{\text{side}} = 8$ maps) to accommodate the $\ell < 47$ multipoles in the analysis, resulting in $N_{\text{pix}} = 3 \times 3072 = 9216$ total number of pixels, further reduced by the application of the analysis mask. In evaluating the likelihood function in Equation (63), the data vector and the noise covariance matrix are fixed, while the signal covariance matrix is recomputed for any given cosmological model to be compared against data, following Equation (23). In practice, only a subsection of S is recomputed, in particular that subsection corresponding to $\ell < 30$. The portion of S corresponding to multipoles $30 \le \ell < 47$ is precomputed from a fixed fiducial model. The choice of the fiducial model does not affect the performance of the likelihood. In fact, at the resolution employed in the large-scale analysis, the sensitivity to multipoles above $\ell \sim 30$ is strongly suppressed.

In 2013, a hybrid approach coupling a MonteCarlo-based approach in temperature (Blackwell-Rao estimator, see section 5.2) to a pixel-based approach in polarization for $\ell < 23$ was adopted [40]. In order to speed up the evaluation of the fully

pixel-based likelihood function for the 2015 release at any given cosmological model, the “brute-force” approach described here has been optimized with the implementation of the Sherman-Morrison-Woodbury formula, described below.

### 5.1.1. Sherman-Morrison-Woodbury Formula

The computational cost required by the pixel-based approach can be dramatically reduced if one consider that only a portion of the signal covariance matrix is reconstructed at any evaluation of the likelihood function. The full covariance matrix can then be decomposed into a varying part, which is function of the theoretical $C_\ell$ (i.e., the power spectra of the theoretical models to be compared against data), and a fixed part given by the fiducial S and the noise covariance matrix. The following step is to further decompose the varying part of the covariance as $S = V^TAV$, via a transformation V that effectively reduces the dimension of the actual evaluation cost from a $N_{\text{pix}} \times N_{\text{pix}}$ inversion to a $n_\lambda \times n_\lambda$ inversion, where $n_\lambda = 2\ell + 1$ is the dimension of the transformed matrix A. The latter is the only matrix that depends on the theoretical $C_\ell$ and, therefore, it is the only matrix to be recomputed and inverted. The fixed portions of the covariance matrix as well as the transformation matrix V can be pre-computed and stored. For the Planck 2015 data release, the application of such mathematical formalisms allowed to speed-up the likelihood evaluation by an order of magnitude. The mathematical details leading to the application of the Sherman-Morrison-Woodbury formula can be found in Aghanim et al. [37] (appendix B.1) and Hinshaw et al. [61] (appendix A.1).

## 5.2. Blackwell-Rao Estimator

An alternative approach to the likelihood evaluation that overcomes the computational cost of an exact likelihood evaluation in pixel space is provided by the combination of the Gibbs sampling method with the Blackwell-Rao estimator. The Gibbs sampling is a MCMC method applied to the estimation of the observed CMB signal from a raw map containing signal, noise, and foreground contamination in a Bayesian framework [62–64]. The crucial output for the subsequent construction of the likelihood estimator is a set of samples of the CMB sky, or more precisely, a set of sample variances of the sky samples obtained with the Gibbs algorithm. The Blackwell-Rao estimator is then built as an average over the set of sample variances.

Let's assume that the observed map m is composed by the CMB signal s and noise n (the following procedure can be generalized to the case in which foreground f also contribute to the total signal, see e.g., [65]). What we really want to evaluate is the joint probability of having a certain s with a certain $C_\ell$ given the observed sky m, i.e., we want to evaluate $\mathcal{P}(s, C_\ell|\mathbf{m})$. A brute force evaluation by computing a grid of s and $C_\ell$ is computationally prohibitive (it is even more prohibitive once one considers the inclusion of foreground into the equation). Another approach is to sample directly from the distribution via specific algorithms. Although it is again computationally unfeasible to sample directly from the joint distribution, as it would require inverting large and

<footnote>$^{11}$The signal in each pixel is the average over the signal at each point within the surface area of the pixel $\Omega_p$, $f_p = \int d\Omega w_p f(\Omega)$, where $w_p = 1/\Omega_p$ is the weight inside the pixel and zero otherwise. The harmonic transform of the weight $w_p$ is the pixel window function. Including the full dependence of the window function on m, $\ell$ can be computationally demanding. However, if the size of the pixel is small compared to the angular resolution of the experiment, the m-dependence can be neglected and the m-averaged window function $w_\ell$ can be defined. The power spectra computed from the pixelised maps $C_\ell^p$ are related to the unpixelised spectra $C_\ell$ via $C_\ell^p = w_\ell^2 C_\ell$.</footnote>
dense covariance matrices, it has been proved that the joint distribution can be reconstructed by iteratively sampling over the conditional distributions $P(s|C_\ell, m)$ and $P(C_\ell|s, m)$. In fact, the conditional distributions are known. They are a multivariate Gaussian (posterior distribution of a Wiener filter) for $P(s|C_\ell, m)$ and an inverse Gamma distribution for $P(C_\ell|s, m)$. As for the foregrounds, their marginal distribution does not usually have an analytic representation. However, it can be reconstructed numerically [65]. The iterative sampling is obtained with the implementation of a specific MCMC sampling algorithm, the Gibbs sampling, and follows these steps (we omit foregrounds for simplicity):

* i) start from a guess power spectrum $C_\ell^0$;

* ii) draw a sample $s^1$ from $P(s|C_\ell^0, m)$;

* iii) draw a sample $C_\ell^1$ from $P(C_\ell|s^1, m)$;

* iv) cycle over step ii) and iii) until convergence is reached.

This algorithm lays the ground for the subsequent evaluation of the posterior $P(C_\ell(\theta)|m)$, where $\theta$ is a set of cosmological parameters of a theoretical cosmological model. In principle, one could reconstruct a histogram of the sampled $C_\ell$ from the Gibbs sampling and use the histogram to interpolate for a given theoretical model. However, this procedure is not efficient. Luckily, the Gibbs sampling allows to construct an efficient and arbitrarily exact estimator of the likelihood function, the so called Blackwell-Rao estimator. The basic idea [64, 66] is to note that the $C_\ell$ only depend on the CMB signal and not on the total sky signal (i.e., once we know the CMB sky, there is no additional information coming from the knowledge of other components), $P(C_\ell|s, m) \propto P(C_\ell|s)$. In addition, the $C_\ell$ depend on the CMB signal only through its variance, i.e., $P(C_\ell|s) \propto P(C_\ell|\hat{C}_\ell)$, where $\hat{C}_\ell$ is the power spectrum of the sample CMB map $s$. Note that we are already familiar with the probability distribution $P(C_\ell|\hat{C}_\ell)$ (see Equation (17)). At this point, we can write down the following chain of equivalent probability integrals (we omit the dependence of $C_\ell$ from $\theta$ for brevity):

$$P(C_{\ell}|\mathbf{m}) = \int P(C_{\ell}, \mathbf{s}|\mathbf{m})d\mathbf{s} \quad (65a)$$

$$= \int P(C_{\ell}|\mathbf{s}, \mathbf{m})P(\mathbf{s}|\mathbf{m})d\mathbf{s} \quad (65b)$$

$$= \int P(C_{\ell}|\hat{C}_{\ell})P(\hat{C}_{\ell}|\mathbf{m})d\hat{C}_{\ell} \quad (65c)$$

$$\sim \frac{1}{N_{\text{Gibbs}}} \sum_{i}^{N_{\text{Gibbs}}} P(C_{\ell}|\hat{C}_{\ell}^{i}) \quad (65d)$$

where $N_{\text{Gibbs}}$ is the number of Gibbs samples evaluated in the component separation analysis, and $\hat{C}_\ell^i$ is the power spectrum of the i-th sample. In other words, the posterior on the cosmological parameters of interest given the observed data can be obtained by averaging the individual posteriors over the Gibbs samples, that can be stored during the Gibbs implementation. The evaluation can be made more accurate by increasing the number of samples. It has to be noted that one only needs to store the variance of the sampled CMB maps $s^i$ up to a given multipole $\ell_{\max}$, i.e.,

it is not necessary to store the much more memory-demanding samples themselves.

The accuracy of the Blackwell-Rao estimator depends on the number of samples required to reach convergence. The number of samples grows very steeply with the $\ell_{\max}$ considered in the analysis [66]. This limitation makes the Blackwell-Rao estimator more suited for the likelihood analysis of CMB large scales (low multipoles). Modified version of the Blackwell-Rao estimator have been proposed that allow to reduce the number of samples required for convergence, and therefore allow to extend the viability of this approach to higher $\ell_{\max}$ [67].

The Gibbs method has been employed by the Planck collaboration as a component separation method to obtain maps of the CMB temperature signal and foreground contaminants [68–70]. The Blackwell-Rao estimator, based on the Gibbs samples so obtained, has been employed by the Planck collaboration as an alternative likelihood method for the analysis of the temperature data at large scales [37, 40]. It has been employed by the WMAP collaboration for the likelihood analysis of the large-scale ($\ell < 32$) temperature data [4, 71].

### 5.3. Simulation-Based Approach

In some cases accurate estimate of the noise are not available and/or the probability distribution of residual systematic effects, in map space, is not perfectly Gaussian, and therefore the probability distribution, in map or harmonic space, cannot be expressed analytically and it needs to be learned from simulations.

We describe in this section a simulation-based approach to evaluate the likelihood function in the low-multipole regime. The likelihood distribution is evaluated through realistic simulations of CMB, noise and possible residual systematics. We report here the main steps:

* the initial step is the simulation of n theoretical auto- or cross-power spectra, related to n different cosmological models, represented by a set of cosmological parameters $\theta_j$.

* For each theoretical power spectrum $C_\ell^{\text{th},i}$, m CMB maps are produced, including noise and other residual contaminants.

* For each map, the corresponding power spectrum is evaluated. Therefore, $k = n \times m$ new simulated power spectra $C_\ell^{\text{sim},k}$ are obtained.

* By histogramming the k simulated power spectra $\ell$-by-$\ell$ and $\theta_j$ by $\theta_j$, the probability $P(C_\ell^{\text{sim}}|C_\ell^{\text{th}})$ is built empirically. In evaluating this probability, it is necessary to interpolate it with a suitable function, in order to smooth the scatter due to the limited available number of simulations. In particular, one can define a low-order polynomial interpolation function $f_\ell^i(C_\ell^{\text{sim}}, C_\ell^{\text{th}})$ of the logarithmic of the $C_\ell^{\text{sim}}$ histogram for each i-th initial power spectrum, such that

$$f_{\ell}^{i}(C_{\ell}^{\text{sim}}, C_{\ell}^{\text{th}}) \approx \log P(C_{\ell}^{\text{sim}}|C_{\ell}^{\text{th}}). \quad (66)$$

* Starting from this approximation, the likelihood function for the observed power spectra $C_\ell^{\text{data}}$ can be finally evaluated. The n couples $(C_\ell^{\text{th},i}, f_\ell^i(C_\ell^{\text{th}}, C_\ell^{\text{data}}))$ can be considered as a tabulated version of the log of the joint probability $\log P(C_\ell^{\text{data}}, C_\ell^{\text{th}})$.
dense covariance matrices, it has been proved that the joint distribution can be reconstructed by iteratively sampling over the conditional distributions $P(s|C_\ell, m)$ and $P(C_\ell|s, m)$. In fact, the conditional distributions are known. They are a multivariate Gaussian (posterior distribution of a Wiener filter) for $P(s|C_\ell, m)$ and an inverse Gamma distribution for $P(C_\ell|s, m)$. As for the foregrounds, their marginal distribution does not usually have an analytic representation. However, it can be reconstructed numerically [65]. The iterative sampling is obtained with the implementation of a specific MCMC sampling algorithm, the Gibbs sampling, and follows these steps (we omit foregrounds for simplicity):

* i) start from a guess power spectrum $C_\ell^0$;

* ii) draw a sample $s^1$ from $P(s|C_\ell^0, m)$;

* iii) draw a sample $C_\ell^1$ from $P(C_\ell|s^1, m)$;

* iv) cycle over step ii) and iii) until convergence is reached.

This algorithm lays the ground for the subsequent evaluation of the posterior $P(C_\ell(\theta)|m)$, where $\theta$ is a set of cosmological parameters of a theoretical cosmological model. In principle, one could reconstruct a histogram of the sampled $C_\ell$ from the Gibbs sampling and use the histogram to interpolate for a given theoretical model. However, this procedure is not efficient. Luckily, the Gibbs sampling allows to construct an efficient and arbitrarily exact estimator of the likelihood function, the so called Blackwell-Rao estimator. The basic idea [64, 66] is to note that the $C_\ell$ only depend on the CMB signal and not on the total sky signal (i.e., once we know the CMB sky, there is no additional information coming from the knowledge of other components), $P(C_\ell|s, m) \propto P(C_\ell|s)$. In addition, the $C_\ell$ depend on the CMB signal only through its variance, i.e., $P(C_\ell|s) \propto P(C_\ell|\hat{C}_\ell)$, where $\hat{C}_\ell$ is the power spectrum of the sample CMB map $s$. Note that we are already familiar with the probability distribution $P(C_\ell|\hat{C}_\ell)$ (see Equation (17)). At this point, we can write down the following chain of equivalent probability integrals (we omit the dependence of $C_\ell$ from $\theta$ for brevity):

$$P(C_{\ell}|\mathbf{m}) = \int P(C_{\ell}, \mathbf{s}|\mathbf{m})d\mathbf{s} \quad (65a)$$

$$= \int P(C_{\ell}|\mathbf{s}, \mathbf{m})P(\mathbf{s}|\mathbf{m})d\mathbf{s} \quad (65b)$$

$$= \int P(C_{\ell}|\hat{C}_{\ell})P(\hat{C}_{\ell}|\mathbf{m})d\hat{C}_{\ell} \quad (65c)$$

$$\sim \frac{1}{N_{\text{Gibbs}}} \sum_{i}^{N_{\text{Gibbs}}} P(C_{\ell}|\hat{C}_{\ell}^{i}) \quad (65d)$$

where $N_{\text{Gibbs}}$ is the number of Gibbs samples evaluated in the component separation analysis, and $\hat{C}_\ell^i$ is the power spectrum of the i-th sample. In other words, the posterior on the cosmological parameters of interest given the observed data can be obtained by averaging the individual posteriors over the Gibbs samples, that can be stored during the Gibbs implementation. The evaluation can be made more accurate by increasing the number of samples. It has to be noted that one only needs to store the variance of the sampled CMB maps $s^i$ up to a given multipole $\ell_{\max}$, i.e.,

it is not necessary to store the much more memory-demanding samples themselves.

The accuracy of the Blackwell-Rao estimator depends on the number of samples required to reach convergence. The number of samples grows very steeply with the $\ell_{\max}$ considered in the analysis [66]. This limitation makes the Blackwell-Rao estimator more suited for the likelihood analysis of CMB large scales (low multipoles). Modified version of the Blackwell-Rao estimator have been proposed that allow to reduce the number of samples required for convergence, and therefore allow to extend the viability of this approach to higher $\ell_{\max}$ [67].

The Gibbs method has been employed by the Planck collaboration as a component separation method to obtain maps of the CMB temperature signal and foreground contaminants [68–70]. The Blackwell-Rao estimator, based on the Gibbs samples so obtained, has been employed by the Planck collaboration as an alternative likelihood method for the analysis of the temperature data at large scales [37, 40]. It has been employed by the WMAP collaboration for the likelihood analysis of the large-scale ($\ell < 32$) temperature data [4, 71].

## 5.3. Simulation-Based Approach

In some cases accurate estimate of the noise are not available and/or the probability distribution of residual systematic effects, in map space, is not perfectly Gaussian, and therefore the probability distribution, in map or harmonic space, cannot be expressed analytically and it needs to be learned from simulations.

We describe in this section a simulation-based approach to evaluate the likelihood function in the low-multipole regime. The likelihood distribution is evaluated through realistic simulations of CMB, noise and possible residual systematics. We report here the main steps:

* the initial step is the simulation of n theoretical auto- or cross-power spectra, related to n different cosmological models, represented by a set of cosmological parameters $\theta_j$.

* For each theoretical power spectrum $C_\ell^{\text{th},i}$, m CMB maps are produced, including noise and other residual contaminants.

* For each map, the corresponding power spectrum is evaluated. Therefore, $k = n \times m$ new simulated power spectra $C_\ell^{\text{sim},k}$ are obtained.

* By histogramming the k simulated power spectra $\ell$-by-$\ell$ and $\theta_j$ by $\theta_j$, the probability $P(C_\ell^{\text{sim}}|C_\ell^{\text{th}})$ is built empirically. In evaluating this probability, it is necessary to interpolate it with a suitable function, in order to smooth the scatter due to the limited available number of simulations. In particular, one can define a low-order polynomial interpolation function $f_\ell^i(C_\ell^{\text{sim}}, C_\ell^{\text{th}})$ of the logarithmic of the $C_\ell^{\text{sim}}$ histogram for each i-th initial power spectrum, such that

$$f_{\ell}^{i}(C_{\ell}^{\text{sim}}, C_{\ell}^{\text{th}}) \approx \log P(C_{\ell}^{\text{sim}}|C_{\ell}^{\text{th}}). \quad (66)$$

* Starting from this approximation, the likelihood function for the observed power spectra $C_\ell^{\text{data}}$ can be finally evaluated. The n couples $(C_\ell^{\text{th},i}, f_\ell^i(C_\ell^{\text{th}}, C_\ell^{\text{data}}))$ can be considered as a tabulated version of the log of the joint probability $\log P(C_\ell^{\text{data}}, C_\ell^{\text{th}})$.