The joint probability can then be interpolated with a suitable low-order polynomial (as done above) for each multipole, $g_{\ell}(C_{\ell}^{\text{data}}, C_{\ell}^{\text{th}})$, such that

$$g_{\ell}(C_{\ell}^{\text{data}}, C_{\ell}^{\text{th}}) \approx \log \mathcal{P}(C_{\ell}^{\text{data}}|C_{\ell}^{\text{th}}). \quad (67)$$

The sum in the multipole range provides the approximation for the log-likelihood, up to a constant,

$$\log \mathcal{L}(\mathbf{C}^{\text{th}}|\mathbf{C}^{\text{data}}) \approx \sum_{\ell=\ell_{\text{min}}}^{\ell_{\text{max}}} g_{\ell}(C_{\ell}^{\text{data}}, C_{\ell}^{\text{th}}) + \text{cost} \quad (68)$$

This approach has been used in the latest analysis of the Planck collaboration [3, 38, 59] to evaluate the low-multipoles polarization likelihood.

We note that this likelihood approach requires the data to be provided in harmonic space as power spectra $\mathbf{C}^{\text{data}}$. A good estimator for the power spectrum at large scales is obtained via the quadratic maximum likelihood (QML) method. We provide a description of the method in section B.3.

As a final remark, we note that the simulation-based approach is very similar in spirit to the class of methods known as Approximate Bayesian Computation (ABC). ABC methods have been initially employed in various fields as a way to bypass the evaluation of the likelihood function with the use of simulated data [72–81]. Recent applications to astrophysics and cosmology include [82–86]. Similarly to ABC, the simulation-based approach approximates the likelihood by comparing simulated data with observed data using the QML approximations to the ML points of the spectra as a distance metric. However, while the likelihood approximation in ABC is usually done by discarding mismatching simulations, the simulation-based approach described here does so by fitting the constructed histogram and then slicing.

# 6. FOREGROUND CONTAMINATION: MODELING AND COMPONENT SEPARATION

We have neglected so far the possibility that the microwave sky contains more components than the CMB. In realistic observations, one has to take into account the fact that the observed signal at a given frequency is a combination of the CMB signal plus additional emissions from so called foreground contaminants. For the range of wavelengths of interest to CMB observations (from few tens to few hundreds GHz), the most relevant contaminants are atmospheric emission (for ground-based experiments in particular) and astrophysical foreground emission. The latter includes Galactic dust, synchrotron and free-free emission and extragalactic contaminants such as clustered and Poisson CIB emission, radio point sources, molecular lines; see e.g. [87, 88] for a description of the multi-component microwave sky and way to create synthetic maps, and [89] for a concise review. The temperature sky is much more composite than the polarization sky. Nevertheless, CMB temperature signal is the dominant component over a wide range of frequencies

and angular scales. On the other hand, polarized foreground can easily dominate over the cosmological signal, especially in the case of CMB B-modes.

To reduce the contamination from atmospheric emission, ground-based CMB experiments are usually placed in specific sites, at high altitudes (to reduce the thickness of the atmosphere above the telescope) and dry locations. Residual contamination can be further removed either by detector pair-differencing (see e.g., [11]) or by modulation of the signal (see e.g., [90]). Balloon-borne and especially satellite missions are less or not at all concerned with atmospheric contaminations. Emission from foreground sources is instead a common issue to all CMB observatories. Of course, it is necessary to account for such foreground contaminations in a proper likelihood analysis, i.e., to decompose the observed map into the individual components. This process is called component separation; various component separation methods exist. They differ for the domain of applicability (pixel space, harmonic space, wavelet space), and for the way data are described. For example, CMB and foregrounds can be parametrized as frequency spectra (parametric methods), can be separated imposing certain conditions to the various components modeled as arbitrary templates (non-parametric), can be separated with no a-priori knowledge of the individual components (blind). A rather complete list of component separation methods can be found in the LAMBDA archive $^{12}$. Here, we cite a few examples that have been also applied to the Planck analyses [91]: ILC [92] (pixel space, internal linear combination of frequency components), SEVEM [93] (pixel space, template-based), NILC [94] (needlet space, internal linear combination), COMMANDER [95] (map domain, parametric), SMICA [96] (harmonic space, non parametric).

After component separation, the cleaned CMB maps may still be affected by residual contamination, which have to be further taken into account. At small scales, this is achieved in various steps, see e.g., the Planck analysis [37]. First, regions of the sky that show an excess contamination from foreground emission are masked away from the analysis. Secondly, the residual contamination in the remaining areas can be modeled at the power spectrum level via specific templates. In other words, one can exploit the fact that the harmonic shape of the foreground is different from the CMB power spectra, and the fact that the spectral dependence of the foreground emission is also different from the CMB spectral dependence. The data vector to be fed in the likelihood function would therefore contain both the CMB signal and the residual contamination, as obtained from observations of the real sky. At the same time, the theoretical spectra to be compared with the data would be given by the sum of the theoretical CMB spectra and the foreground templates. It is crucial that multi-frequency channels are available in order to efficiently fit for both the CMB and the foreground contaminants simultaneously.

The above prescription implicitly assumed that all the information about foreground emission in harmonic space is fully captured by their power spectrum. However, there is no

<footnote>$^{12}$Component separation methods and related references can be accessed at the following url: https://lambda.gsfc.nasa.gov/toolbox/tb_comp_separation.cfm</footnote>
The joint probability can then be interpolated with a suitable low-order polynomial (as done above) for each multipole, $g_{\ell}(C_{\ell}^{\text{data}}, C_{\ell}^{\text{th}})$, such that

$$g_{\ell}(C_{\ell}^{\text{data}}, C_{\ell}^{\text{th}}) \approx \log \mathcal{P}(C_{\ell}^{\text{data}}|C_{\ell}^{\text{th}}). \quad (67)$$

The sum in the multipole range provides the approximation for the log-likelihood, up to a constant,

$$\log \mathcal{L}(\mathbf{C}^{\text{th}}|\mathbf{C}^{\text{data}}) \approx \sum_{\ell=\ell_{\text{min}}}^{\ell_{\text{max}}} g_{\ell}(C_{\ell}^{\text{data}}, C_{\ell}^{\text{th}}) + \text{cost} \quad (68)$$

This approach has been used in the latest analysis of the Planck collaboration [3, 38, 59] to evaluate the low-multipoles polarization likelihood.

We note that this likelihood approach requires the data to be provided in harmonic space as power spectra $C^{\text{data}}$. A good estimator for the power spectrum at large scales is obtained via the quadratic maximum likelihood (QML) method. We provide a description of the method in section B.3.

As a final remark, we note that the simulation-based approach is very similar in spirit to the class of methods known as Approximate Bayesian Computation (ABC). ABC methods have been initially employed in various fields as a way to bypass the evaluation of the likelihood function with the use of simulated data [72–81]. Recent applications to astrophysics and cosmology include [82–86]. Similarly to ABC, the simulation-based approach approximates the likelihood by comparing simulated data with observed data using the QML approximations to the ML points of the spectra as a distance metric. However, while the likelihood approximation in ABC is usually done by discarding mismatching simulations, the simulation-based approach described here does so by fitting the constructed histogram and then slicing.

# 6. FOREGROUND CONTAMINATION: MODELING AND COMPONENT SEPARATION

We have neglected so far the possibility that the microwave sky contains more components than the CMB. In realistic observations, one has to take into account the fact that the observed signal at a given frequency is a combination of the CMB signal plus additional emissions from so called foreground contaminants. For the range of wavelengths of interest to CMB observations (from few tens to few hundreds GHz), the most relevant contaminants are atmospheric emission (for ground-based experiments in particular) and astrophysical foreground emission. The latter includes Galactic dust, synchrotron and free-free emission and extragalactic contaminants such as clustered and Poisson CIB emission, radio point sources, molecular lines; see e.g. [87, 88] for a description of the multi-component microwave sky and way to create synthetic maps, and [89] for a concise review. The temperature sky is much more composite than the polarization sky. Nevertheless, CMB temperature signal is the dominant component over a wide range of frequencies

and angular scales. On the other hand, polarized foreground can easily dominate over the cosmological signal, especially in the case of CMB B-modes.

To reduce the contamination from atmospheric emission, ground-based CMB experiments are usually placed in specific sites, at high altitudes (to reduce the thickness of the atmosphere above the telescope) and dry locations. Residual contamination can be further removed either by detector pair-differencing (see e.g., [11]) or by modulation of the signal (see e.g., [90]). Balloon-borne and especially satellite missions are less or not at all concerned with atmospheric contaminations. Emission from foreground sources is instead a common issue to all CMB observatories. Of course, it is necessary to account for such foreground contaminations in a proper likelihood analysis, i.e., to decompose the observed map into the individual components. This process is called component separation; various component separation methods exist. They differ for the domain of applicability (pixel space, harmonic space, wavelet space), and for the way data are described. For example, CMB and foregrounds can be parametrized as frequency spectra (parametric methods), can be separated imposing certain conditions to the various components modeled as arbitrary templates (non-parametric), can be separated with no a-priori knowledge of the individual components (blind). A rather complete list of component separation methods can be found in the LAMBDA archive $^{12}$. Here, we cite a few examples that have been also applied to the Planck analyses [91]: ILC [92] (pixel space, internal linear combination of frequency components), SEVEM [93] (pixel space, template-based), NILC [94] (needlet space, internal linear combination), COMMANDER [95] (map domain, parametric), SMICA [96] (harmonic space, non parametric).

After component separation, the cleaned CMB maps may still be affected by residual contamination, which have to be further taken into account. At small scales, this is achieved in various steps, see e.g., the Planck analysis [37]. First, regions of the sky that show an excess contamination from foreground emission are masked away from the analysis. Secondly, the residual contamination in the remaining areas can be modeled at the power spectrum level via specific templates. In other words, one can exploit the fact that the harmonic shape of the foreground is different from the CMB power spectra, and the fact that the spectral dependence of the foreground emission is also different from the CMB spectral dependence. The data vector to be fed in the likelihood function would therefore contain both the CMB signal and the residual contamination, as obtained from observations of the real sky. At the same time, the theoretical spectra to be compared with the data would be given by the sum of the theoretical CMB spectra and the foreground templates. It is crucial that multi-frequency channels are available in order to efficiently fit for both the CMB and the foreground contaminants simultaneously.

The above prescription implicitly assumed that all the information about foreground emission in harmonic space is fully captured by their power spectrum. However, there is no

<footnote>$^{12}$Component separation methods and related references can be accessed at the following url: https://lambda.gsfc.nasa.gov/toolbox/tb_comp_separation.cfm</footnote>
reason to believe that foreground are Gaussian distributed, and indeed they are not. The assumption that is usually made is that most of the non-Gaussian contribution is removed from the maps once the foreground-saturated regions are masked away. The non-Gaussianity of the remaining contaminations can be (and actually are) neglected at the likelihood level. This assumption is demonstrated to be reasonably accurate via dedicated simulations (see e.g., discussions in Aghanim et al. [37] and Ade et al. [40]) for the current generation of cosmological surveys. At the same time, a huge effort is devoted to the study of the propagation of uncertainties in the foreground modeling and removal to the final constraints on cosmological parameters in the context of the high-sensitivity next generation CMB surveys (see e.g., [19]).

Foreground contaminations are of course present at large angular scales as well, and must be taken into account when realistic data are analyzed. Taking the results from the Planck collaboration as an example [37, 40], the foreground treatment at large scales is somehow different from the prescription described in the case of small scales. At low multipoles in temperature, the cleaned CMB map is obtained via Gibbs sampling (see section 5.2) implemented in COMMANDER and the bright region along the Galactic plane is further masked for likelihood analysis (7% of the sky). In polarization, the foreground-cleaning is implemented via template-fitting. Residual contamination can be taken into account in the construction of the noise covariance matrix to be employed in the likelihood function.

To conclude this section, we would like to emphasize that the topic of foreground modeling and component separation is extremely wide and the list of references reported in this manuscript is far from being exhaustive. The interested reader is encouraged to consider these references as a mere starting point. Further details can be accessed from the specific collaboration papers that describe the corresponding data processing, and from the overview papers of science forecasts by upcoming collaborations that describe their simulated pipelines of data reduction.

# 7. COMPARING LIKELIHOOD PERFORMANCE

In this section, we compare the performance of various likelihood approaches introduced in section 3 for the high-$\ell$ regime and in section 5 for the low-$\ell$ regime. In particular, we are interested in testing the different likelihood performance with respect to the ability to produce unbiased constraints on cosmological parameters. Where relevant, we also compare the different performance in terms of computational costs. In what follows, we first discuss the comparison of likelihood approaches in the small-scale regime (section 7.1), and then move to the presentation of the results in the large-scale regime (section 7.3).

## 7.1. Small-Scale Regime–Full Sky

In this section we test the performance of the different likelihood approximations at high-$\ell$ on simulated data and compare the results with the assessment presented in section 3.2, which was based on analytic arguments. We mainly want to show whether

adopting a particular power spectrum likelihood approximation, when estimating cosmological parameters, may introduce a bias in the parameters recovered values or/and a misestimation of their associated error bars. In the following, we neglect the impact of foregrounds, and we consider CMB plus noise full-sky maps, for which we derive the temperature and polarization angular power spectra following Equation (15).

The simulated dataset consists of CMB plus noise maps. Specifically, we generate a set of 1000 maps of the CMB sky, $\mathbf{m} = (T, Q, U)$, drawn as Gaussian random realizations of fiducial temperature and polarization power spectra, that correspond to a set of given cosmological parameters. The full-sky maps are generated at a resolution of 3.4 arcmins ($N_{side} = 1024$ in the HEALPix scheme [60]) and smoothed with a symmetric Gaussian beam of FWHM 10 arcmins $^{13}$. To each of these maps we add a simulated realization of white isotropic noise corresponding to a noise level in polarization of $\sigma_n = 1\mu\text{K arcmin}$, which is in the ballpark of what is expected from future CMB experiments. In this way temperature anisotropy maps are signal dominated across almost all the multipoles that are relevant for primary anisotropies, up to $\ell \lesssim 2000$, compatibly with present and forthcoming measurements.

### 7.1.1. Temperature-Only Results

Focusing on temperature alone, we adopt the different likelihood approximations of section 3 to estimate the cosmological parameters from the power spectra of the simulations. For simplicity we only fit for two parameters, specifically $A_s$ and $n_s$, that are respectively the amplitude and the spectral index of the power spectrum of primordial density fluctuations. We generate a grid of model $C_\ell$'s keeping all the other $\Lambda$CDM parameters fixed, while letting $A_s$ and $n_s$ vary in a broad range of values around the input fiducial model, $C_{f\ell}$. Given the angular power spectrum of each simulation, $\hat{C}_\ell$, we evaluate the likelihood of the models at each multipole, $\mathcal{L}(C_\ell)$, accounting for the noise contribution as described in section 3.3. Since here we are considering the small scale regime, the total likelihood for the set of cosmological parameters $\theta = (A_s, n_s)$ is obtained by summing the log-likelihoods over the range of multipoles $\ell = [30, 2000]: \ln\mathcal{L}(\theta) \propto \sum_{\ell} \ln\mathcal{L}(C_\ell)$.

Once we have the total likelihood, assuming flat priors on parameters, we can estimate the expectation value of each parameter following Equation (78). By building histograms of the parameter values obtained from the 1000 simulated maps, we can reconstruct the posterior distribution of each parameter in a “frequentist” fashion. In particular, since the simulations input parameters are known, we can study the distribution of the biases in the estimated parameters, defined as the shifts with respect to the input fiducial values, normalized to the $1\sigma$ marginal error of the posterior distribution. For each simulation, the bias reads $\Delta p = (\hat{p} - p^{\text{in}})/\sigma_p$, where $p$ can be either $A_s$ or $n_s$. The distributions we derived are shown in **Figure 1** for the different likelihood approximations. They should be centered in zero, if

<footnote>$^{13}$Roughly speaking, the resolution identifies the level of the discretization (i.e., the number of pixels). The beam smoothing simulates the angular resolution of the experiment, and exponentially suppresses the scales below the beam width. Note that the beam width is larger than the size of the pixel.</footnote>
reason to believe that foreground are Gaussian distributed, and indeed they are not. The assumption that is usually made is that most of the non-Gaussian contribution is removed from the maps once the foreground-saturated regions are masked away. The non-Gaussianity of the remaining contaminations can be (and actually are) neglected at the likelihood level. This assumption is demonstrated to be reasonably accurate via dedicated simulations (see e.g., discussions in Aghanim et al. [37] and Ade et al. [40]) for the current generation of cosmological surveys. At the same time, a huge effort is devoted to the study of the propagation of uncertainties in the foreground modeling and removal to the final constraints on cosmological parameters in the context of the high-sensitivity next generation CMB surveys (see e.g., [19]).

Foreground contaminations are of course present at large angular scales as well, and must be taken into account when realistic data are analyzed. Taking the results from the Planck collaboration as an example [37, 40], the foreground treatment at large scales is somehow different from the prescription described in the case of small scales. At low multipoles in temperature, the cleaned CMB map is obtained via Gibbs sampling (see section 5.2) implemented in COMMANDER and the bright region along the Galactic plane is further masked for likelihood analysis (7% of the sky). In polarization, the foreground-cleaning is implemented via template-fitting. Residual contamination can be taken into account in the construction of the noise covariance matrix to be employed in the likelihood function.

To conclude this section, we would like to emphasize that the topic of foreground modeling and component separation is extremely wide and the list of references reported in this manuscript is far from being exhaustive. The interested reader is encouraged to consider these references as a mere starting point. Further details can be accessed from the specific collaboration papers that describe the corresponding data processing, and from the overview papers of science forecasts by upcoming collaborations that describe their simulated pipelines of data reduction.

# 7. COMPARING LIKELIHOOD PERFORMANCE

In this section, we compare the performance of various likelihood approaches introduced in section 3 for the high-$\ell$ regime and in section 5 for the low-$\ell$ regime. In particular, we are interested in testing the different likelihood performance with respect to the ability to produce unbiased constraints on cosmological parameters. Where relevant, we also compare the different performance in terms of computational costs. In what follows, we first discuss the comparison of likelihood approaches in the small-scale regime (section 7.1), and then move to the presentation of the results in the large-scale regime (section 7.3).

## 7.1. Small-Scale Regime–Full Sky

In this section we test the performance of the different likelihood approximations at high-$\ell$ on simulated data and compare the results with the assessment presented in section 3.2, which was based on analytic arguments. We mainly want to show whether

adopting a particular power spectrum likelihood approximation, when estimating cosmological parameters, may introduce a bias in the parameters recovered values or/and a misestimation of their associated error bars. In the following, we neglect the impact of foregrounds, and we consider CMB plus noise full-sky maps, for which we derive the temperature and polarization angular power spectra following Equation (15).

The simulated dataset consists of CMB plus noise maps. Specifically, we generate a set of 1000 maps of the CMB sky, $\mathbf{m} = (T, Q, U)$, drawn as Gaussian random realizations of fiducial temperature and polarization power spectra, that correspond to a set of given cosmological parameters. The full-sky maps are generated at a resolution of 3.4 arcmins ($N_{side} = 1024$ in the HEALPix scheme [60]) and smoothed with a symmetric Gaussian beam of FWHM 10 arcmins $^{13}$. To each of these maps we add a simulated realization of white isotropic noise corresponding to a noise level in polarization of $\sigma_n = 1\mu\text{K arcmin}$, which is in the ballpark of what is expected from future CMB experiments. In this way temperature anisotropy maps are signal dominated across almost all the multipoles that are relevant for primary anisotropies, up to $\ell \lesssim 2000$, compatibly with present and forthcoming measurements.

### 7.1.1. Temperature-Only Results

Focusing on temperature alone, we adopt the different likelihood approximations of section 3 to estimate the cosmological parameters from the power spectra of the simulations. For simplicity we only fit for two parameters, specifically $A_s$ and $n_s$, that are respectively the amplitude and the spectral index of the power spectrum of primordial density fluctuations. We generate a grid of model $C_\ell$'s keeping all the other $\Lambda$CDM parameters fixed, while letting $A_s$ and $n_s$ vary in a broad range of values around the input fiducial model, $C_{f\ell}$. Given the angular power spectrum of each simulation, $\hat{C}_\ell$, we evaluate the likelihood of the models at each multipole, $\mathcal{L}(C_\ell)$, accounting for the noise contribution as described in section 3.3. Since here we are considering the small scale regime, the total likelihood for the set of cosmological parameters $\theta = (A_s, n_s)$ is obtained by summing the log-likelihoods over the range of multipoles $\ell = [30, 2000]: \ln\mathcal{L}(\theta) \propto \sum_{\ell} \ln\mathcal{L}(C_\ell)$.

Once we have the total likelihood, assuming flat priors on parameters, we can estimate the expectation value of each parameter following Equation (78). By building histograms of the parameter values obtained from the 1000 simulated maps, we can reconstruct the posterior distribution of each parameter in a “frequentist” fashion. In particular, since the simulations input parameters are known, we can study the distribution of the biases in the estimated parameters, defined as the shifts with respect to the input fiducial values, normalized to the $1\sigma$ marginal error of the posterior distribution. For each simulation, the bias reads $\Delta p = (\hat{p} - p^{\text{in}})/\sigma_p$, where p can be either $A_s$ or $n_s$. The distributions we derived are shown in **Figure 1** for the different likelihood approximations. They should be centered in zero, if

<footnote>$^{13}$Roughly speaking, the resolution identifies the level of the discretization (i.e., the number of pixels). The beam smoothing simulates the angular resolution of the experiment, and exponentially suppresses the scales below the beam width. Note that the beam width is larger than the size of the pixel.</footnote>
there is no bias in the mean recovered parameters value, and they should have unit variance, if $\sigma_p$, as derived from the standard bayesian analysis, is consistent with the *true* error.

Consistently with what has been presented in section 3.2, we find that the Log-Normal together with the symmetric and improper Gaussian approximations show significant bias in the recovered parameters. The bias is due to the fact that the approximations fail in capturing the correct skewness of the likelihood function. The symmetric Gaussian approximation appears to be the one with the largest associated bias and the wider distribution. For this approximation the covariance is proportional to the *measured* power spectrum squared, $\hat{C}_{\ell}^2$. This means that the uncertainties turn out to be larger than they should be if $\hat{C}_{\ell}$ fluctuates upward with respect to the input fiducial model, and viceversa if $\hat{C}_{\ell}$ has a downward fluctuation. As a consequence, upward fluctuations are given less weight in the likelihood than downward fluctuations. This translates into an overall downward bias on the amplitude of the power spectrum, and thus on $A_s$, as it can be easily seen in **Figure 1**. The one-third-two-third approximation is a good match to the exact likelihood at third order, as shown in **Figure 2**, and it appears to provide unbiased parameter results. This is the likelihood approximation adopted by the WMAP team for the analysis of the temperature anisotropies on small angular scales [39]. The offset-lognormal also results in unbiased constraints (bottom left panel in **Figure 1**), as expected since the offset $a_{TT} = -1/4$ has been chosen to match the exact likelihood up to third order in the expansion, see section 3.2. Also the fiducial and determinant Gaussian approximations provide unbiased cosmological parameter estimates, despite having a wrong shape at each single $\ell$ with respect to the exact likelihood. This is because at each particular $\ell$ the shape of the likelihood can randomly be wrong upwards or downwards with respect to the *true* $C_{\ell}$ (see **Figure 2**), so when integrating over the entire range of multipoles there is negligible bias left. An obvious advantage of the fiducial Gaussian approximation with respect to the determinant Gaussian is that the covariance matrix can be pre-computed, inverted, and then kept fixed while sampling the parameter space. Aside from speeding up the computations, this makes it easier to account for additional uncertainties, non-Gaussianities and correlations in the data by simply adding extra terms in the pre-computed covariance matrix. The fiducial Gaussian approximation is the one used for the official small-scale likelihood analysis of the Planck data [40], and it has also been adopted for the cosmological analysis of data collected by the ground-based experiments ACT [97] and SPT [98]. Note, however, that this approximation works well if the models we are trying to fit are smooth, and if the fiducial model used to build the covariance matrix is sufficiently close to the *true* model. In the tests presented here the fiducial model entering the covariance matrix actually coincides with the *true* model, i.e., the power spectrum from which the simulated maps are generated. Obviously, in this case the covariance matrix is an optimal description of the simulated data. However, when dealing with real data we are not in the same fortunate position, since the *true* power spectrum underlying the observed sky is unknown and is exactly what we intend to estimate. In a real situation

we can only deduce a fiducial model for the covariance at the best of our knowledge from both theory and observations. We can also try to improve the accuracy of a first guess covariance matrix by iterating the cosmological parameter extraction a few times, upgrading the covariance matrix at each iteration. Note that, we did not include the Hamimeche-Lewis approximation in the tests, the reason being that on the full-sky it coincides with the exact likelihood.

In order to perform an even more stringent validation of the likelihood approximations, we can estimate the bias of the mean of the recovered parameters from all the simulations: $((\hat{p}) - p^{\text{in}})/(\sigma_p/\sqrt{N_{\text{sims}}})$, where the shift with respect to the input is normalized by the standard deviation of the mean, and $N_{\text{sims}} = 1000$. This test confirms, at high-significance level, that Log-Normal, symmetric and improper Gaussian approximations provide biased results, while the other approximations are unbiased. For the latter likelihood approximations we find that the estimated uncertainties on parameters, $\sigma_p$, agree with those derived from the exact likelihood at the sub-precent level. Furthermore, we find that the standard deviation of the values of $(\hat{p} - p^{\text{in}})/\sigma_p$ is consistent to 1, with a precision well within the accuracy allowed by the finite number of simulations, i.e., $1/\sqrt{2N_{\text{sims}}} = 0.022$. This means that the uncertainties estimated from the bayesian analysis, $\sigma_p$, are consistent with the scatter of the parameters of the simulations.

### 7.1.2. Temperature and Polarization Joint Results

For the fiducial Gaussian approximation, we repeat the test of parameters recovery from the full-sky simulated maps also including the polarization power spectra. As we have already commented, this likelihood approximation works well provided that the fiducial power spectra used to build the covariance matrix are close enough to the *true* power spectra. In order to investigate the impact on the final parameters of misestimating the covariance matrix, we change the fiducial model entering the likelihood with power spectra at the edge of the grid of models introduced above. These power spectra correspond to parameters more than $10\sigma$ away from the input parameters. As expected, results show some sensitivity to the choice of the model, however not for the bias, but rather for the parameters marginal error bars. In fact, using the “wrong” covariance matrix leaves parameters estimates unbiased. This can be noticed in **Figure 3**, where we show the biases on $A_s$ and $n_s$ estimated from the simulations, both using the covariance matrix with the *correct* fiducial model and the one with the modified model. The latter covariance matrix, however, does not match exactly the signal in the simulated dataset, and as a consequence the estimated $1\sigma$ marginal error bars on the cosmological parameters do not agree with the standard deviation from the simulations, showing $a \sim 15\%$ mismatch.

## 7.2. Small-Scale Regime—Cut-Sky Tests

We repeat the above analysis introducing some further level of complication, specifically we assume to work with partial-sky and anisotropic noise data. We generate noise simulations which have a different level of noise variance for each pixel, thus mimicking the effect of a realistic scan strategy for which some regions of the
there is no bias in the mean recovered parameters value, and they should have unit variance, if $\sigma_p$, as derived from the standard bayesian analysis, is consistent with the true error.

Consistently with what has been presented in section 3.2, we find that the Log-Normal together with the symmetric and improper Gaussian approximations show significant bias in the recovered parameters. The bias is due to the fact that the approximations fail in capturing the correct skewness of the likelihood function. The symmetric Gaussian approximation appears to be the one with the largest associated bias and the wider distribution. For this approximation the covariance is proportional to the *measured* power spectrum squared, $\hat{C}_{\ell}^2$. This means that the uncertainties turn out to be larger than they should be if $\hat{C}_{\ell}$ fluctuates upward with respect to the input fiducial model, and viceversa if $\hat{C}_{\ell}$ has a downward fluctuation. As a consequence, upward fluctuations are given less weight in the likelihood than downward fluctuations. This translates into an overall downward bias on the amplitude of the power spectrum, and thus on $A_s$, as it can be easily seen in **Figure 1**. The one-third-two-third approximation is a good match to the exact likelihood at third order, as shown in **Figure 2**, and it appears to provide unbiased parameter results. This is the likelihood approximation adopted by the WMAP team for the analysis of the temperature anisotropies on small angular scales [39]. The offset-lognormal also results in unbiased constraints (bottom left panel in **Figure 1**), as expected since the offset $a_{TT} = -1/4$ has been chosen to match the exact likelihood up to third order in the expansion, see section 3.2. Also the fiducial and determinant Gaussian approximations provide unbiased cosmological parameter estimates, despite having a wrong shape at each single $\ell$ with respect to the exact likelihood. This is because at each particular $\ell$ the shape of the likelihood can randomly be wrong upwards or downwards with respect to the true $C_{\ell}$ (see **Figure 2**), so when integrating over the entire range of multipoles there is negligible bias left. An obvious advantage of the fiducial Gaussian approximation with respect to the determinant Gaussian is that the covariance matrix can be pre-computed, inverted, and then kept fixed while sampling the parameter space. Aside from speeding up the computations, this makes it easier to account for additional uncertainties, non-Gaussianities and correlations in the data by simply adding extra terms in the pre-computed covariance matrix. The fiducial Gaussian approximation is the one used for the official small-scale likelihood analysis of the Planck data [40], and it has also been adopted for the cosmological analysis of data collected by the ground-based experiments ACT [97] and SPT [98]. Note, however, that this approximation works well if the models we are trying to fit are smooth, and if the fiducial model used to build the covariance matrix is sufficiently close to the true model. In the tests presented here the fiducial model entering the covariance matrix actually coincides with the true model, i.e., the power spectrum from which the simulated maps are generated. Obviously, in this case the covariance matrix is an optimal description of the simulated data. However, when dealing with real data we are not in the same fortunate position, since the true power spectrum underlying the observed sky is unknown and is exactly what we intend to estimate. In a real situation

we can only deduce a fiducial model for the covariance at the best of our knowledge from both theory and observations. We can also try to improve the accuracy of a first guess covariance matrix by iterating the cosmological parameter extraction a few times, upgrading the covariance matrix at each iteration. Note that, we did not include the Hamimeche-Lewis approximation in the tests, the reason being that on the full-sky it coincides with the exact likelihood.

In order to perform an even more stringent validation of the likelihood approximations, we can estimate the bias of the mean of the recovered parameters from all the simulations: $((\hat{p}) - p^{\text{in}})/(\sigma_p/\sqrt{N_{\text{sims}}})$, where the shift with respect to the input is normalized by the standard deviation of the mean, and $N_{\text{sims}} = 1000$. This test confirms, at high-significance level, that Log-Normal, symmetric and improper Gaussian approximations provide biased results, while the other approximations are unbiased. For the latter likelihood approximations we find that the estimated uncertainties on parameters, $\sigma_p$, agree with those derived from the exact likelihood at the sub-precent level. Furthermore, we find that the standard deviation of the values of $(\hat{p} - p^{\text{in}})/\sigma_p$ is consistent to 1, with a precision well within the accuracy allowed by the finite number of simulations, i.e., $1/\sqrt{2N_{\text{sims}}} = 0.022$. This means that the uncertainties estimated from the bayesian analysis, $\sigma_p$, are consistent with the scatter of the parameters of the simulations.

### 7.1.2. Temperature and Polarization Joint Results

For the fiducial Gaussian approximation, we repeat the test of parameters recovery from the full-sky simulated maps also including the polarization power spectra. As we have already commented, this likelihood approximation works well provided that the fiducial power spectra used to build the covariance matrix are close enough to the true power spectra. In order to investigate the impact on the final parameters of misestimating the covariance matrix, we change the fiducial model entering the likelihood with power spectra at the edge of the grid of models introduced above. These power spectra correspond to parameters more than $10\sigma$ away from the input parameters. As expected, results show some sensitivity to the choice of the model, however not for the bias, but rather for the parameters marginal error bars. In fact, using the “wrong” covariance matrix leaves parameters estimates unbiased. This can be noticed in **Figure 3**, where we show the biases on $A_s$ and $n_s$ estimated from the simulations, both using the covariance matrix with the correct fiducial model and the one with the modified model. The latter covariance matrix, however, does not match exactly the signal in the simulated dataset, and as a consequence the estimated $1\sigma$ marginal error bars on the cosmological parameters do not agree with the standard deviation from the simulations, showing $a \sim 15\%$ mismatch.

## 7.2. Small-Scale Regime—Cut-Sky Tests

We repeat the above analysis introducing some further level of complication, specifically we assume to work with partial-sky and anisotropic noise data. We generate noise simulations which have a different level of noise variance for each pixel, thus mimicking the effect of a realistic scan strategy for which some regions of the
<figure><img src="image_2.png" /></figure>

**FIGURE 1** | 2D histograms of the biases on the parameters expectation values estimated from 1000 simulated full-sky maps of CMB temperature plus noise. As described in more detail in the text, the biases are computed with respect to the input parameter values and normalized to the 1 σ uncertainty associated to the parameter, $\Delta p = (\hat{p} - p^n)/\sigma_p$, where p can be either $A_s$ or $n_s$. The blue circle indicates a null bias, while the red square is the center of the distribution. Note that for the cases plotted in the first row the two almost completely overlap.

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

**FIGURE 2** | Comparison of some of the likelihood approximations for different multipole regimes, as marked in the plots, where the x-axis corresponds to $D_l = l(l+1)C_l/(2\pi)$.
<figure><img src="image_3.png" /></figure>

**FIGURE 1** | 2D histograms of the biases on the parameters expectation values estimated from 1000 simulated full-sky maps of CMB temperature plus noise. As described in more detail in the text, the biases are computed with respect to the input parameter values and normalized to the 1 σ uncertainty associated to the parameter, $\Delta p = (\hat{p} - p^n)/\sigma_p$, where p can be either $A_s$ or $n_s$. The blue circle indicates a null bias, while the red square is the center of the distribution. Note that for the cases plotted in the first row the two almost completely overlap.

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

**FIGURE 2** | Comparison of some of the likelihood approximations for different multipole regimes, as marked in the plots, where the x-axis corresponds to $D_\ell = \ell(\ell + 1)C_\ell/(2\pi)$.
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** | Histograms of the biases for the parameters expectation values estimated from the power spectra of 1000 simulated full-sky maps of CMB temperature and polarization, plus white isotropic noise. Biases are defined with respect to the input parameter values and normalized to the 1 σ uncertainty on the parameter, $\Delta\rho = (\hat{\rho} - p^n)/\sigma_p$, where $p$ can be either $A_s$ or $n_s$. We adopt the fiducial Gaussian likelihood approximation, in green the case in which the fiducial model entering the covariance matrix matches the input model of the simulations, whereas in blue a case in which the two models differ as explained in the text.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 4** | Histograms of the biases for the parameters expectation values estimated from the power spectra of 1000 simulated partial-sky maps of CMB temperature anisotropies plus anisotropic noise. Biases are defined analogously to **Figure 3**. We compare results from the fiducial gaussian and the Hamimeche-Lewis likelihood approximations.</figcaption></figure>

sky are observed more often than others. To the simulated CMB plus noise maps we then apply a galactic sky mask, that removes the region of the galactic plane where the emission of foregrounds is expected to dominate over the CMB signal. For the sake of this test, we focus on temperature alone and we use a galactic mask that leaves for the analysis 73% of sky. The mask has been apodised with a Gaussian taper that smooths sharp edges and, thus, it helps in localizing the mask power in multipole space. The angular power spectra of the cut-sky maps are extracted using an estimator that corrects for the loss of modes due to the masking and which is based on the pseudo-$C_\ell$ formalism. For completeness the estimator is described in **Appendix B**. The covariance matrix associated to these power spectra has been estimated using the analytic approximation given in Efstathiou [27], and assuming as fiducial model the same model from which the simulations have been generated.

For the parameters recovery test on the cut-sky we explore the fiducial gaussian and Hamimeche-Lewis likelihood approximations. The latter in the single field regime reduces to:

$$-2 \ln \mathcal{L} \approx \sum_{\ell\ell'} \left[ g(\hat{C}_\ell/C_\ell)C_{\ell\ell'} \left[M_f^{-1}\right]_{\ell\ell'} \left[g(\hat{C}_{\ell'}/C_{\ell'})C_{\ell\ell'}\right] \right], \quad (69)$$

where the function $g(x)$ has been introduced in Equation (51), and $[M_f]_{\ell\ell'}$ is the covariance matrix of the $\hat{C}_\ell$ evaluated for $C_{\ell\ell}$.

Following the same procedure described in the previous subsection, we fit for the cosmological parameters $A_s$ and $n_s$. We find that, also under more realistic conditions, the recovered parameters for the fiducial gaussian are unbiased well within the precision allowed by the finite number of Monte Carlo simulations. For the Hamimeche-Lewis approximation, instead, we detect a small bias at the level of 20 and 11% of the sigma on the parameter for $A_s$ and $n_s$ respectively, see **Figure 4**. For both likelihood approximations, however, we find that the parameters marginal errors from the bayesian analysis are consistent with the standard deviation of the simulations, and thus they appear to be a good description of the true uncertainties.

Despite the small average bias, the Hamimeche-Lewis approximation is expected to be more robust against the choice of the fiducial model entering the likelihood, see [32]. In order to verify this point, when computing the likelihood and the covariance matrix of the $\hat{C}_\ell$, we use a fiducial model $C_{\ell\ell'}$ that corresponds to parameters about 10 σ away from those used to generate the simulations. As we have already shown in the previous subsection, in this case the fiducial Gaussian
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Histograms of the biases for the parameters expectation values estimated from the power spectra of 1000 simulated full-sky maps of CMB temperature and polarization, plus white isotropic noise. Biases are defined with respect to the input parameter values and normalized to the 1 σ uncertainty on the parameter, $\Delta\rho = (\hat{\rho} - p^n)/\sigma_p$, where $p$ can be either $A_s$ or $n_s$. We adopt the fiducial Gaussian likelihood approximation, in green the case in which the fiducial model entering the covariance matrix matches the input model of the simulations, whereas in blue a case in which the two models differ as explained in the text.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Histograms of the biases for the parameters expectation values estimated from the power spectra of 1000 simulated partial-sky maps of CMB temperature anisotropies plus anisotropic noise. Biases are defined analogously to **Figure 3**. We compare results from the fiducial gaussian and the Hamimeche-Lewis likelihood approximations.</figcaption></figure>

sky are observed more often than others. To the simulated CMB plus noise maps we then apply a galactic sky mask, that removes the region of the galactic plane where the emission of foregrounds is expected to dominate over the CMB signal. For the sake of this test, we focus on temperature alone and we use a galactic mask that leaves for the analysis 73% of sky. The mask has been apodised with a Gaussian taper that smooths sharp edges and, thus, it helps in localizing the mask power in multipole space. The angular power spectra of the cut-sky maps are extracted using an estimator that corrects for the loss of modes due to the masking and which is based on the pseudo-$C_\ell$ formalism. For completeness the estimator is described in **Appendix B**. The covariance matrix associated to these power spectra has been estimated using the analytic approximation given in Efstathiou [27], and assuming as fiducial model the same model from which the simulations have been generated.

For the parameters recovery test on the cut-sky we explore the fiducial gaussian and Hamimeche-Lewis likelihood approximations. The latter in the single field regime reduces to:

$$ -2 \ln \mathcal{L} \approx \sum_{\ell\ell'} \left[ g(\hat{C}_\ell/C_\ell)C_{\ell\ell'} \left[ M_f^{-1} \right]_{\ell\ell'} \left[ g(\hat{C}_{\ell'}/C_{\ell'})C_{\ell\ell'} \right] \right], \quad (69) $$

where the function $g(x)$ has been introduced in Equation (51), and $[M_f]_{\ell\ell'}$ is the covariance matrix of the $\hat{C}_\ell$ evaluated for $C_{\ell\ell}$.

Following the same procedure described in the previous subsection, we fit for the cosmological parameters $A_s$ and $n_s$. We find that, also under more realistic conditions, the recovered parameters for the fiducial gaussian are unbiased well within the precision allowed by the finite number of Monte Carlo simulations. For the Hamimeche-Lewis approximation, instead, we detect a small bias at the level of 20 and 11% of the sigma on the parameter for $A_s$ and $n_s$ respectively, see **Figure 4**. For both likelihood approximations, however, we find that the parameters marginal errors from the bayesian analysis are consistent with the standard deviation of the simulations, and thus they appear to be a good description of the true uncertainties.

Despite the small average bias, the Hamimeche-Lewis approximation is expected to be more robust against the choice of the fiducial model entering the likelihood, see [32]. In order to verify this point, when computing the likelihood and the covariance matrix of the $\hat{C}_\ell$, we use a fiducial model $C_{\ell\ell'}$ that corresponds to parameters about 10 σ away from those used to generate the simulations. As we have already shown in the previous subsection, in this case the fiducial Gaussian
<figure><img src="image_6.png" /><figcaption>**FIGURE 5** | Histograms of the recovered values of τ for the three likelihoods considered: pixel based (blue), simulation based (orange), and Hamimeche-Lewis (green).</figcaption></figure>

approximation provides still unbiased results, but the bayesian error bars on the parameters are not a good description of the true uncertainties. They differ from the standard deviation of the parameters from the simulations by about 14% on $A_s$ and 11% on $n_s$, respectively. On the contrary, estimates with the Hamimeche-Lewis approximation still show the same level of average bias, but the uncertainties on the parameters are better characterized and they are in agreement with the standard deviation from the simulations.

Furthermore, since it is customary to assess the goodness of the parameters fit with the chi-square statistics, we may use the value of the likelihood at the best-fit to define an effective chi-square as $\chi_{eff}^2 = \min[-2\ln\mathcal{L}(\theta|\hat{C}_{\ell})]$. If we now compare the values we get for each likelihood approximation when simply varying the fiducial model, we find $\Delta\chi_{eff}^2$ of order a few hundreds for the fiducial Gaussian approximation, and $\Delta\chi_{eff}^2$ of a few tens for the Hamimeche-Lewis. As a consequence, trying to assess the goodness of the fit with the latter likelihood can surely provide more stable results, regardless of the fiducial model. However, it is worth stressing that in the present test we used a $C_{f\ell}$ that deviates significantly from the model behind the simulations, had we chosen a model closer to the simulations, also the fiducial Gaussian approximation would have provided sensible results. For further discussion of this topic refer to the **appendix B** of Hamimeche and Lewis [32].

### 7.3. Large-Scale Regime

In this subsection, we focus on the comparison of different likelihood approaches devoted to the analysis of large-scale CMB data. We limit the comparison of large-scale likelihoods to the polarization signal only, since this is the main target of future CMB missions [99]. Moreover, the standard approaches used in temperature, i.e., pixel-based and Blackwell-Rao (see section 5.2)

likelihoods, have been extensively characterized and validated by the Planck Collaboration, see e.g., [37, 38].

We consider three approaches: pixel-based (see section 5.1), HL (see Equations 49a-54a) and simulation-based (see section 5.3) likelihood. As we did for the small scale regime, we show whether adopting a particular likelihood approximation when estimating cosmological parameters may introduce a bias, either in the recovered values of the parameters or on their associated error bars.

We start by generating a set of 1000 maps of the CMB sky drawn as Gaussian random realizations of a single fiducial power spectrum corresponding to a set of known cosmological parameters. The full sky maps are generated on a $N_{side} = 16$ HEALPix grid, which roughly corresponds to a resolution of 3.7 degrees, smoothed with a cosine window function [38, 100]. To each of these maps we add a realization of white isotropic noise corresponding to $\sigma_N = 0.01 \mu\text{K}^2$ on a pixel at our resolution. We choose this particular noise level since it is small but not completely negligible with respect to the typical peak-to-peak CMB signal in a model with a reionization optical depth $\tau = 0.055$, that would be roughly ~ 0.5$\mu\text{K}$. We do not include foreground and systematics effect residuals in our simulations. We process the maps through a QML code computing the auto spectra of all our simulations.

To test the performance of the likelihood approximations, for simplicity, we only fit the reionization optical depth $\tau$ parameter keeping fixed the overall amplitude of the perturbations as parametrized through $A_s \exp(-2\tau)$. For each realization we compute the mean value of $\tau$. The histograms of the $\tau$ values, recovered using the three different approximations, are shown in **Figure 5**. The distributions all show a bias smaller than 5% of $\sigma(\tau)$, perfectly compatible with the resolution of 1000 simulations, (i.e., $1/\sqrt{1000} \sim 3\%$). The Hamimeche-Lewis likelihood and the simulation-based likelihood, both built on the QML estimates, perform similarly and are substantially unbiased. Likewise, the pixel-based likelihood results unbiased as it should, as also already extensively demonstrated by both the WMAP [101] and Planck [37] collaborations. The comparison with the pixel-based likelihood allows also to validate the width of the distributions for the other two likelihood approximations, which are compatible with the pixel based at the 2% level, within the resolution of the Montecarlo (i.e., $1/\sqrt{2 \times 1000} \sim 2\%$).

## 8. CONCLUSIONS

CMB science is nowadays a mature yet still flourishing branch of physics. The ultimate results from the Planck satellite and from large-aperture ground-based experiments (ACT and SPT) fully characterize the CMB temperature field with great accuracy up to very small angular scales, where the contaminations from Galactic and extra-Galactic foreground emissions become dominant. In addition, they have just scratched the surface of the rich information contained in the CMB polarization fields, also the target of small-aperture telescopes (BICEP/Keck, POLARBEAR/Simons Array, SPIDER). A large number of upcoming CMB surveys (Simons Observatory,
<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | Histograms of the recovered values of τ for the three likelihoods considered: pixel based (blue), simulation based (orange), and Hamimeche-Lewis (green).</figcaption></figure>

approximation provides still unbiased results, but the bayesian error bars on the parameters are not a good description of the true uncertainties. They differ from the standard deviation of the parameters from the simulations by about 14% on $A_s$ and 11% on $n_s$, respectively. On the contrary, estimates with the Hamimeche-Lewis approximation still show the same level of average bias, but the uncertainties on the parameters are better characterized and they are in agreement with the standard deviation from the simulations.

Furthermore, since it is customary to assess the goodness of the parameters fit with the chi-square statistics, we may use the value of the likelihood at the best-fit to define an effective chi-square as $\chi_{eff}^2 = \min[-2\ln\mathcal{L}(\theta|\hat{C}_{\ell})]$. If we now compare the values we get for each likelihood approximation when simply varying the fiducial model, we find $\Delta\chi_{eff}^2$ of order a few hundreds for the fiducial Gaussian approximation, and $\Delta\chi_{eff}^2$ of a few tens for the Hamimeche-Lewis. As a consequence, trying to assess the goodness of the fit with the latter likelihood can surely provide more stable results, regardless of the fiducial model. However, it is worth stressing that in the present test we used a $C_{f\ell}$ that deviates significantly from the model behind the simulations, had we chosen a model closer to the simulations, also the fiducial Gaussian approximation would have provided sensible results. For further discussion of this topic refer to the **appendix B** of Hamimeche and Lewis [32].

### 7.3. Large-Scale Regime

In this subsection, we focus on the comparison of different likelihood approaches devoted to the analysis of large-scale CMB data. We limit the comparison of large-scale likelihoods to the polarization signal only, since this is the main target of future CMB missions [99]. Moreover, the standard approaches used in temperature, i.e., pixel-based and Blackwell-Rao (see section 5.2)

likelihoods, have been extensively characterized and validated by the Planck Collaboration, see e.g., [37, 38].

We consider three approaches: pixel-based (see section 5.1), HL (see Equations 49a-54a) and simulation-based (see section 5.3) likelihood. As we did for the small scale regime, we show whether adopting a particular likelihood approximation when estimating cosmological parameters may introduce a bias, either in the recovered values of the parameters or on their associated error bars.

We start by generating a set of 1000 maps of the CMB sky drawn as Gaussian random realizations of a single fiducial power spectrum corresponding to a set of known cosmological parameters. The full sky maps are generated on a $N_{side} = 16$ HEALPix grid, which roughly corresponds to a resolution of 3.7 degrees, smoothed with a cosine window function [38, 100]. To each of these maps we add a realization of white isotropic noise corresponding to $\sigma_N = 0.01 \mu\text{K}^2$ on a pixel at our resolution. We choose this particular noise level since it is small but not completely negligible with respect to the typical peak-to-peak CMB signal in a model with a reionization optical depth $\tau = 0.055$, that would be roughly ~ 0.5$\mu\text{K}$. We do not include foreground and systematics effect residuals in our simulations. We process the maps through a QML code computing the auto spectra of all our simulations.

To test the performance of the likelihood approximations, for simplicity, we only fit the reionization optical depth $\tau$ parameter keeping fixed the overall amplitude of the perturbations as parametrized through $A_s \exp(-2\tau)$. For each realization we compute the mean value of $\tau$. The histograms of the $\tau$ values, recovered using the three different approximations, are shown in **Figure 5**. The distributions all show a bias smaller than 5% of $\sigma(\tau)$, perfectly compatible with the resolution of 1000 simulations, (i.e., $1/\sqrt{1000} \sim 3\%$). The Hamimeche-Lewis likelihood and the simulation-based likelihood, both built on the QML estimates, perform similarly and are substantially unbiased. Likewise, the pixel-based likelihood results unbiased as it should, as also already extensively demonstrated by both the WMAP [101] and Planck [37] collaborations. The comparison with the pixel-based likelihood allows also to validate the width of the distributions for the other two likelihood approximations, which are compatible with the pixel based at the 2% level, within the resolution of the Montecarlo (i.e., $1/\sqrt{2 \times 1000} \sim 2\%$).

## 8. CONCLUSIONS

CMB science is nowadays a mature yet still flourishing branch of physics. The ultimate results from the Planck satellite and from large-aperture ground-based experiments (ACT and SPT) fully characterize the CMB temperature field with great accuracy up to very small angular scales, where the contaminations from Galactic and extra-Galactic foreground emissions become dominant. In addition, they have just scratched the surface of the rich information contained in the CMB polarization fields, also the target of small-aperture telescopes (BICEP/Keck, POLARBEAR/Simons Array, SPIDER). A large number of upcoming CMB surveys (Simons Observatory,
CMB-S4, LiteBIRD) will continue the journey toward the understanding of the early phases of our Universe and its subsequent evolution.

In the highly sophisticated process that aims to efficiently compress the information contained in the raw observed data into maps, spectra, and eventually constraints on cosmological parameters, a key ingredient is the likelihood function $L(d|\theta)$, i.e., the probability of observing a certain data collection d given a certain model $\theta$ (see section A). In the context of CMB analysis, d can be either a CMB map or CMB spectra, while $\theta$ can be the set of cosmological parameters describing the cosmological model under scrutiny. In the simple case of full-sky observations, in absence of foreground contaminations and late-time-Universe effects on the CMB distribution, an exact form of the likelihood function in both real (multivariate Gaussian) and harmonic space (Wishart distribution) can be easily derived, see section 2. Complications arising from realistic observations and data analysis, such as complicated noise properties, foreground obscuration, limited sky coverage, computational costs, etc, may spoil the assumptions on which the derivation of the exact likelihood functions mentioned above is based, or may limit the actual evaluation of the exact functions. Therefore, realistic analysis employ appropriate approximations of the exact likelihood functions, see sections 3, 5. The choice of the likelihood depends on multiple factors, such as sky-fraction retained, data resolution, computational costs, signal-to-noise properties. For example, it may (and actually did) happen that a certain likelihood approximation could work for a certain experiment given its sensitivity, and yet it could turn out to provide biased results for another, more sensitive experiment.

In this review, we have summarized the basics of CMB statistics (section 2) that lead to the definition of the exact likelihoods in real and harmonic space. Then, we have moved to the descriptions of the most common likelihood approximations employed by various CMB collaborations. For the sake of simplicity, we have separated the approximations better suited to the analysis of small angular scales (higher harmonic multipoles) in section 3 and those that better represent the data at large angular scales (lower multipoles) in section 5. Although we have mostly assumed an idealized scenario of isotropic instrumental noise, Gaussian beam, unlensed CMB spectra, and absence of foreground for pedagogical reasons, we have commented about the impact of realistic deviations from the aforementioned scenario, in such a way that the reader is aware of their existence and of the extensive effort devoted to their mitigation. In particular, non-trivial modifications to the primary CMB statistics induced by the gravitational lensing of CMB photons by the evolving large-scale structure in the Universe have been discussed in section 4. In section 6, we have briefly discussed the main foreground contaminants in the microwave sky and mentioned methods of component separation aimed to disentangle the cosmological signal from Galactic and extra-Galactic emissions.

In order to provide concrete examples of the different performance of various likelihood methods, we have reported

results of the comparison between different likelihood approximations in section 7. In particular, we have tested the property of the likelihood function to be unbiased, i.e., to produce a posterior distribution of cosmological parameters that matches the true cosmology (the input values, in the case of simulated data) both in terms of the mean and of the variance of the distribution. Again, we have distinguished between approximations at small scales (section 7.1) and at large scales (section 7.3).

In Appendix, the interested reader can find more details about some basics of statistical inference (section A) and about the most common methods adopted to estimate CMB power spectra in the (realistic) cut-sky regime (section B).

The aim of this review is to bring together and summarize a large amount of information related to a crucial aspect of CMB data analysis. The final scope is to provide a handy, yet self-consistent document to those who are approaching the field and those who are interested in learning some basic aspects of the field, and help them to navigate the vast literature produced over the past decades on these topics. If you have gone so far in your reading, we hope to have reached, at least partly, our goal.

## AUTHOR CONTRIBUTIONS

All the authors contributed in a relevant way to the development and finalization of the manuscript. The ordering of the authors (two tiers, alphabetic order in each tier) reflects the relative contribution, with authors within each tier having contributed equally.

## ACKNOWLEDGMENTS

The authors would like to thank D. Poletti for useful discussions. We acknowledge support from the ASI grant 2016-24-H.0 COSMOS Attività di studio per la comunità scientifica di cosmologia. MG acknowledges support by Argonne National Laboratory (ANL). This work was supported by the U.S. Department of Energy, Office of Science, under contract number DE-AC02-06CH11357. ML and MM acknowledge support from INFN through the InDark and Gruppo IV fundings. MM is supported by the program for young researchers Rita Levi Montalcini year 2015. LS acknowledges support from the ERC-StG ClustersXcosmo grant agreement 716762, and from the postdoctoral grant from Centre National d'Études Spatiales (CNES). LC received funding from the European Union's Horizon 2020 research and innovation programme under grant number 776282.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2020.00015/full#supplementary-material
CMB-S4, LiteBIRD) will continue the journey toward the understanding of the early phases of our Universe and its subsequent evolution.

In the highly sophisticated process that aims to efficiently compress the information contained in the raw observed data into maps, spectra, and eventually constraints on cosmological parameters, a key ingredient is the likelihood function $L(d|\theta)$, i.e., the probability of observing a certain data collection d given a certain model $\theta$ (see section A). In the context of CMB analysis, d can be either a CMB map or CMB spectra, while $\theta$ can be the set of cosmological parameters describing the cosmological model under scrutiny. In the simple case of full-sky observations, in absence of foreground contaminations and late-time-Universe effects on the CMB distribution, an exact form of the likelihood function in both real (multivariate Gaussian) and harmonic space (Wishart distribution) can be easily derived, see section 2. Complications arising from realistic observations and data analysis, such as complicated noise properties, foreground obscuration, limited sky coverage, computational costs, etc, may spoil the assumptions on which the derivation of the exact likelihood functions mentioned above is based, or may limit the actual evaluation of the exact functions. Therefore, realistic analysis employ appropriate approximations of the exact likelihood functions, see sections 3, 5. The choice of the likelihood depends on multiple factors, such as sky-fraction retained, data resolution, computational costs, signal-to-noise properties. For example, it may (and actually did) happen that a certain likelihood approximation could work for a certain experiment given its sensitivity, and yet it could turn out to provide biased results for another, more sensitive experiment.

In this review, we have summarized the basics of CMB statistics (section 2) that lead to the definition of the exact likelihoods in real and harmonic space. Then, we have moved to the descriptions of the most common likelihood approximations employed by various CMB collaborations. For the sake of simplicity, we have separated the approximations better suited to the analysis of small angular scales (higher harmonic multipoles) in section 3 and those that better represent the data at large angular scales (lower multipoles) in section 5. Although we have mostly assumed an idealized scenario of isotropic instrumental noise, Gaussian beam, unlensed CMB spectra, and absence of foreground for pedagogical reasons, we have commented about the impact of realistic deviations from the aforementioned scenario, in such a way that the reader is aware of their existence and of the extensive effort devoted to their mitigation. In particular, non-trivial modifications to the primary CMB statistics induced by the gravitational lensing of CMB photons by the evolving large-scale structure in the Universe have been discussed in section 4. In section 6, we have briefly discussed the main foreground contaminants in the microwave sky and mentioned methods of component separation aimed to disentangle the cosmological signal from Galactic and extra-Galactic emissions.

In order to provide concrete examples of the different performance of various likelihood methods, we have reported

results of the comparison between different likelihood approximations in section 7. In particular, we have tested the property of the likelihood function to be unbiased, i.e., to produce a posterior distribution of cosmological parameters that matches the true cosmology (the input values, in the case of simulated data) both in terms of the mean and of the variance of the distribution. Again, we have distinguished between approximations at small scales (section 7.1) and at large scales (section 7.3).

In Appendix, the interested reader can find more details about some basics of statistical inference (section A) and about the most common methods adopted to estimate CMB power spectra in the (realistic) cut-sky regime (section B).

The aim of this review is to bring together and summarize a large amount of information related to a crucial aspect of CMB data analysis. The final scope is to provide a handy, yet self-consistent document to those who are approaching the field and those who are interested in learning some basic aspects of the field, and help them to navigate the vast literature produced over the past decades on these topics. If you have gone so far in your reading, we hope to have reached, at least partly, our goal.

## AUTHOR CONTRIBUTIONS

All the authors contributed in a relevant way to the development and finalization of the manuscript. The ordering of the authors (two tiers, alphabetic order in each tier) reflects the relative contribution, with authors within each tier having contributed equally.

## ACKNOWLEDGMENTS

The authors would like to thank D. Poletti for useful discussions. We acknowledge support from the ASI grant 2016-24-H.0 COSMOS Attività di studio per la comunità scientifica di cosmologia. MG acknowledges support by Argonne National Laboratory (ANL). This work was supported by the U.S. Department of Energy, Office of Science, under contract number DE-AC02-06CH11357. ML and MM acknowledge support from INFN through the InDark and Gruppo IV fundings. MM is supported by the program for young researchers Rita Levi Montalcini year 2015. LS acknowledges support from the ERC-StG ClustersXcosmo grant agreement 716762, and from the postdoctoral grant from Centre National d'Études Spatiales (CNES). LC received funding from the European Union's Horizon 2020 research and innovation programme under grant number 776282.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2020.00015/full#supplementary-material
REFERENCES

1. Hu W, Dodelson S. Cosmic microwave background anisotropies. *Annu Rev Astron Astrophys.* (2002) **40**:171-216. doi: 10.1146/annurev.astro.40.060401.093926

2. Kolb EW, Turner MS. The early Universe. *Front Phys.* (1990) **69**:1-547.

3. Aghanim N, Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. VI. Cosmological parameters. arXiv [Preprint]. arXiv:1807.06209 (2018).

4. Hinshaw G, Larson D, Komatsu E, Spergel DN, Bennett CL, Dunkley J, et al. Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: cosmological parameter results. Astrophys J. (2012) **208**:19. doi: 10.1088/0067-0049/2 08/2/19

5. Bennett CL, Banday A, Gorski KM, Hinshaw G, Jackson P, Keegstra P, et al. Four year COBE DMR cosmic microwave background observations: maps and basic results. Astrophys J. (1996) **464**:L1-4. doi: 10.1086/310075

6. Fraisse AA, Ade PAR, Amiri M, Benton SJ, Bock JJ, Bond JR, et al. SPIDER: probing the early Universe with a suborbital polarimeter. J Cosmol Astropart Phys. (2013) **1304**:047. doi: 10.1088/1475-7516/2013/04/047

7. Aboobaker AM, Ade P, Araujo D, Aubin F, Baccigalupi C, Bao C, et al. The EBEX balloon-borne experiment? Optics, receiver, and polarimetry. Astrophys J Suppl. (2018) **239**:7. doi: 10.3847/1538-4365/aae434

8. de Bernardis P, Ade PAR, Bock JJ, Bond JR, Borrill J, Boscaleri A, et al. A flat universe from high resolution maps of the cosmic microwave background radiation. Nature. (2000) **404**:955-9. doi: 10.1038/35010035

9. Louis T, Grace E, Hasselfield M, Lungu M, Maurin L, Addison GE, et al. The atacama cosmology telescope: two-season ACTPol spectra and parameters. J Cosmol Astropart Phys. (2017) **1706**:031. doi: 10.1088/1475-7516/2017/06/031

10. Henning JW, Sayre JT, Reichardt CL, Ade PAR, Anderson AJ, Austermann JE, et al. Measurements of the temperature and E-mode polarization of the CMB from 500 square degrees of SPTpol data. Astrophys J. (2018) **852**:97. doi: 10.3847/1538-4357/aa9ff4

11. Ade PAR, Ahmed Z, Aikin RW, Alexander KD, Barkats D, Benton SJ, et al. BICEP2 / Keck Array x: constraints on primordial gravitational waves using Planck, WMAP, and New BICEP2/Keck observations through the 2015 season. Phys Rev Lett. (2018) **121**:221301. doi: 10.1103/PhysRevLett.121.221301

12. Ade PAR, Aguilar M, Akiba Y, Arnold K, Baccigalupi C, Barron D, et al. A measurement of the cosmic microwave background B-Mode polarization power spectrum at sub-degree scales from 2 years of POLARBEAR Data. Astrophys J. (2017) **848**:121. doi: 10.3847/1538-4357/ aa8e9f

13. Essinger-Hileman T, Ali A, Amiri M, Appel JW, Araujo D, Bennett CL, et al. CLASS: the cosmology large angular scale surveyor. Proc SPIE Int Soc Opt Eng. (2014) **9153**:915311. doi: 10.1117/12.2056701

14. Akrami Y, Arroja F, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. I. Overview and the cosmological legacy of Planck arXiv [Preprint]. arXiv:1807.06205 (2018).

15. Henderson SW, Allison R, Austermann J, Baildon T, Battaglia N, Beall JA, et al. Advanced ACTPol cryogenic detector arrays and readout. J Low Temp Phys. (2016) **184**:772-9. doi: 10.1007/s10909-016-1575-z

16. Benson BA, Ade PAR, Ahmed Z, Allen SW, Arnold K, Austermann JE, et al. SPT-3G: a next-generation cosmic microwave background polarization experiment on the south pole telescope. Proc SPIE Int Soc Opt Eng. (2014) **9153**:91531P. doi: 10.1117/12.2057305

17. Suzuki A, Ade PAR, Akiba Y, Aleman C, Arnold K, Baccigalupi C, et al. The POLARBEAR-2 and the simons array experiment. J Low Temp Phys. (2016) **184**:805-10. doi: 10.1007/s10909-015-1425-4

18. Bergman AS, Ade PAR, Akers S, Amiri M, Austermann JA, Beall JA, et al. 280 GHz focal plane unit design and characterization for the spider-2 suborbital polarimeter. J Low Temp Phys. (2018) **193**:1075-84. doi: 10.1007/s10909-018-2065-2

19. Abazajian KN, Adshead P, Ahmed Z, Allen SW, Alonso D, Arnold KS, et al. CMB-S4 Science Book, 1st Edn. arXiv [Preprint]. arXiv:1610.02743 (2016).

20. Ade PAR, Aguirre J, Ahmed Z, Aiola S, Ali A, Alonso D, et al. The simons observatory: science goals and forecasts. J Cosmol Astropart Phys. (2019) **1902**:056. doi: 10.1088/1475-7516/2019/02/056

21. Matsumura T, Akiba Y, Borrill J, Chinone Y, Dobbs M, Fuke H, et al. Mission design of LiteBIRD. J Low Temp Phys. (2014) **176**:733. doi: 10.1007/s10909-013-0996-1

22. Delabrouille J, de Bernardis P, Bouchet FR, Achúcarro A, Ade PAR, Allison R, et al. Exploring cosmic origins with CORE: survey requirements and mission design. J Cosmol Astropart Phys. (2018) **1804**:014. doi: 10.1088/1475-7516/2018/04/014

23. Hanany S, Bennett C, Dodelson S, Page L, Bartlett J, Battaglia N, et al. PICO: probe of inflation and cosmic origins. arXiv [Preprint]. arXiv:1902.10541 (2019).

24. Penzias AA, Wilson RW. A Measurement of excess antenna temperature at 4080-Mc/s. Astrophys J. (1965) **142**:419-21. doi: 10.1086/1 48307

25. Akrami Y, Arroja F, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. IX. Constraints on primordial non-Gaussianity arXiv [Preprint]. arXiv:1905.05697 (2019).

26. Zaldarriaga M, Seljak U. An all sky analysis of polarization in the microwave background. Phys Rev D. (1997) **55**:1830-40. doi: 10.1103/PhysRevD.55.1830

27. Efstathiou G. Myths and truths concerning estimation of power spectra Mon Not R Astron Soc. (2004) **349**:603. doi: 10.1111/j.1365-2966.2004.07530.x

28. Zaldarriaga M. The polarization of the cosmic microwave background. arXiv:astro-ph/0305272 (2003).

29. Percival WJ, Brown ML. Likelihood methods for the combined analysis of CMB temperature and polarisation power spectra Mon Not R Astron Soc. (2006) **372**:1104. doi: 10.1111/j.1365-2966.2006.10910.x

30. Tegmark M, de Oliveira-Costa A. How to measure CMB polarization power spectra without losing information. Phys Rev D. (2001) **64**:063001. doi: 10.1103/PhysRevD.64.063001

31. Upham RE, Whittaker L, Brown ML. Exact joint likelihood of pseudo-$C_\ell$ estimates from correlated Gaussian cosmological fields. Mon Not R Astron Soc. (2020) **491**:3165. doi: 10.1093/mnras/stz3225

32. Hamimeche S, Lewis A. "Likelihood analysis of CMB temperature and polarization power spectra. Phys Rev D. (2008) **77**:103013. doi: 10.1103/PhysRevD.77.103013

33. Bond JR, Jaffe AH, Knox LE. Radical compression of cosmic microwave background data. Astrophys J. (2000) **533**:19. doi: 10.1086/308625

34. Smith S, Challinor A, Rocha G. What can be learned from the lensed cosmic microwave background B-mode polarization power spectrum? Phys Rev D. (2006) **73**:023517. doi: 10.1103/PhysRevD.73.023517

35. Lewis A, Challinor A, Turok N. Analysis of CMB polarization on an incomplete sky. Phys Rev D. (2001) **65**:023505. doi: 10.1103/PhysRevD.65.023505

36. Slosar A, Seljak U, Makarov A. Exact likelihood evaluations and foreground marginalization in low resolution WMAP data. Phys Rev D. (2004) **69**:123003. doi: 10.1103/PhysRevD.69.123003

37. Aghanim N, Arnaud M, Ashdown M, Aumont J, Baccigalupi C, Banday AJ, et al. Planck 2015 results. XI. CMB power spectra, likelihoods, and robustness of parameters. Astron Astrophys. (2016) **594**:A11. doi: 10.1051/0004-6361/201526926

38. Aghanim N, Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. V. CMB power spectra and likelihoods. arXiv [Preprint]. arXiv:1907.12875 (2019).

39. Verde L, Peiris HV, Spergel DN, Nolta M, Bennett CL, Halpern M, et al. First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: parameter estimation methodology. Astrophys J Suppl. (2003) **148**:195. doi: 10.1086/377335

40. Ade PAR, Aghanim N, Armitage-Caplan C, Arnaud M, Ashdown M, Atrio-Barandela F, et al. Planck 2013 results. XV. CMB power spectra and likelihood. Astron Astrophys. (2014) **571**:A15. doi: 10.1015/0004-6361/201321573

41. Lewis A, Challinor A. Weak gravitational lensing of the CMB. Phys Rept. (2006) **429**:1-65. doi: 10.1016/j.physrep.2006.03.002

42. Ade PAR, Aghanim N, Armitage-Caplan C, Arnaud M, Ashdown M, Atrio-Barandela F, et al. Planck 2013 results. XVII. Gravitational
REFERENCES

1. Hu W, Dodelson S. Cosmic microwave background anisotropies. *Annu Rev Astron Astrophys.* (2002) **40**:171-216. doi: 10.1146/annurev.astro.40.060401.093926

2. Kolb EW, Turner MS. The early Universe. *Front Phys.* (1990) **69**:1-547.

3. Aghanim N, Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. VI. Cosmological parameters. arXiv [Preprint]. arXiv:1807.06209 (2018).

4. Hinshaw G, Larson D, Komatsu E, Spergel DN, Bennett CL, Dunkley J, et al. Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: cosmological parameter results. Astrophys J. (2012) **208**:19. doi: 10.1088/0067-0049/2 08/2/19

5. Bennett CL, Banday A, Gorski KM, Hinshaw G, Jackson P, Keegstra P, et al. Four year COBE DMR cosmic microwave background observations: maps and basic results. Astrophys J. (1996) **464**:L1-4. doi: 10.1086/310075

6. Fraisse AA, Ade PAR, Amiri M, Benton SJ, Bock JJ, Bond JR, et al. SPIDER: probing the early Universe with a suborbital polarimeter. J Cosmol Astropart Phys. (2013) **1304**:047. doi: 10.1088/1475-7516/2013/04/047

7. Aboobaker AM, Ade P, Araujo D, Aubin F, Baccigalupi C, Bao C, et al. The EBEX balloon-borne experiment? Optics, receiver, and polarimetry. Astrophys J Suppl. (2018) **239**:7. doi: 10.3847/1538-4365/aae434

8. de Bernardis P, Ade PAR, Bock JJ, Bond JR, Borrill J, Boscaleri A, et al. A flat universe from high resolution maps of the cosmic microwave background radiation. Nature. (2000) **404**:955-9. doi: 10.1038/35010035

9. Louis T, Grace E, Hasselfield M, Lungu M, Maurin L, Addison GE, et al. The atacama cosmology telescope: two-season ACTPol spectra and parameters. J Cosmol Astropart Phys. (2017) **1706**:031. doi: 10.1088/1475-7516/2017/06/031

10. Henning JW, Sayre JT, Reichardt CL, Ade PAR, Anderson AJ, Austermann JE, et al. Measurements of the temperature and E-mode polarization of the CMB from 500 square degrees of SPTpol data. Astrophys J. (2018) **852**:97. doi: 10.3847/1538-4357/aa9ff4

11. Ade PAR, Ahmed Z, Aikin RW, Alexander KD, Barkats D, Benton SJ, et al. BICEP2 / Keck Array x: constraints on primordial gravitational waves using Planck, WMAP, and New BICEP2/Keck observations through the 2015 season. Phys Rev Lett. (2018) **121**:221301. doi: 10.1103/PhysRevLett.121.221301

12. Ade PAR, Aguilar M, Akiba Y, Arnold K, Baccigalupi C, Barron D, et al. A measurement of the cosmic microwave background B-Mode polarization power spectrum at sub-degree scales from 2 years of POLARBEAR Data. Astrophys J. (2017) **848**:121. doi: 10.3847/1538-4357/ aa8e9f

13. Essinger-Hileman T, Ali A, Amiri M, Appel JW, Araujo D, Bennett CL, et al. CLASS: the cosmology large angular scale surveyor. Proc SPIE Int Soc Opt Eng. (2014) **9153**:915311. doi: 10.1117/12.2056701

14. Akrami Y, Arroja F, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. I. Overview and the cosmological legacy of Planck arXiv [Preprint]. arXiv:1807.06205 (2018).

15. Henderson SW, Allison R, Austermann J, Baildon T, Battaglia N, Beall JA, et al. Advanced ACTPol cryogenic detector arrays and readout. J Low Temp Phys. (2016) **184**:772-9. doi: 10.1007/s10909-016-1575-z

16. Benson BA, Ade PAR, Ahmed Z, Allen SW, Arnold K, Austermann JE, et al. SPT-3G: a next-generation cosmic microwave background polarization experiment on the south pole telescope. Proc SPIE Int Soc Opt Eng. (2014) **9153**:91531P. doi: 10.1117/12.2057305

17. Suzuki A, Ade PAR, Akiba Y, Aleman C, Arnold K, Baccigalupi C, et al. The POLARBEAR-2 and the simons array experiment. J Low Temp Phys. (2016) **184**:805-10. doi: 10.1007/s10909-015-1425-4

18. Bergman AS, Ade PAR, Akers S, Amiri M, Austermann JA, Beall JA, et al. 280 GHz focal plane unit design and characterization for the spider-2 suborbital polarimeter. J Low Temp Phys. (2018) **193**:1075-84. doi: 10.1007/s10909-018-2065-2

19. Abazajian KN, Adshead P, Ahmed Z, Allen SW, Alonso D, Arnold KS, et al. CMB-S4 Science Book, 1st Edn. arXiv [Preprint]. arXiv:1610.02743 (2016).

20. Ade PAR, Aguirre J, Ahmed Z, Aiola S, Ali A, Alonso D, et al. The simons observatory: science goals and forecasts. J Cosmol Astropart Phys. (2019) **1902**:056. doi: 10.1088/1475-7516/2019/02/056

21. Matsumura T, Akiba Y, Borrill J, Chinone Y, Dobbs M, Fuke H, et al. Mission design of LiteBIRD. J Low Temp Phys. (2014) **176**:733. doi: 10.1007/s10909-013-0996-1

22. Delabrouille J, de Bernardis P, Bouchet FR, Achúcarro A, Ade PAR, Allison R, et al. Exploring cosmic origins with CORE: survey requirements and mission design. J Cosmol Astropart Phys. (2018) **1804**:014. doi: 10.1088/1475-7516/2018/04/014

23. Hanany S, Bennett C, Dodelson S, Page L, Bartlett J, Battaglia N, et al. PICO: probe of inflation and cosmic origins. arXiv [Preprint]. arXiv:1902.10541 (2019).

24. Penzias AA, Wilson RW. A Measurement of excess antenna temperature at 4080-Mc/s. Astrophys J. (1965) **142**:419-21. doi: 10.1086/1 48307

25. Akrami Y, Arroja F, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. IX. Constraints on primordial non-Gaussianity arXiv [Preprint]. arXiv:1905.05697 (2019).

26. Zaldarriaga M, Seljak U. An all sky analysis of polarization in the microwave background. Phys Rev D. (1997) **55**:1830-40. doi: 10.1103/PhysRevD.55.1830

27. Efstathiou G. Myths and truths concerning estimation of power spectra Mon Not R Astron Soc. (2004) **349**:603. doi: 10.1111/j.1365-2966.2004.07530.x

28. Zaldarriaga M. The polarization of the cosmic microwave background. arXiv:astro-ph/0305272 (2003).

29. Percival WJ, Brown ML. Likelihood methods for the combined analysis of CMB temperature and polarisation power spectra Mon Not R Astron Soc. (2006) **372**:1104. doi: 10.1111/j.1365-2966.2006.10910.x

30. Tegmark M, de Oliveira-Costa A. How to measure CMB polarization power spectra without losing information. Phys Rev D. (2001) **64**:063001. doi: 10.1103/PhysRevD.64.063001

31. Upham RE, Whittaker L, Brown ML. Exact joint likelihood of pseudo-$C_\ell$ estimates from correlated Gaussian cosmological fields. Mon Not R Astron Soc. (2020) **491**:3165. doi: 10.1093/mnras/stz3225

32. Hamimeche S, Lewis A. "Likelihood analysis of CMB temperature and polarization power spectra. Phys Rev D. (2008) **77**:103013. doi: 10.1103/PhysRevD.77.103013

33. Bond JR, Jaffe AH, Knox LE. Radical compression of cosmic microwave background data. Astrophys J. (2000) **533**:19. doi: 10.1086/308625

34. Smith S, Challinor A, Rocha G. What can be learned from the lensed cosmic microwave background B-mode polarization power spectrum? Phys Rev D. (2006) **73**:023517. doi: 10.1103/PhysRevD.73.023517

35. Lewis A, Challinor A, Turok N. Analysis of CMB polarization on an incomplete sky. Phys Rev D. (2001) **65**:023505. doi: 10.1103/PhysRevD.65.023505

36. Slosar A, Seljak U, Makarov A. Exact likelihood evaluations and foreground marginalization in low resolution WMAP data. Phys Rev D. (2004) **69**:123003. doi: 10.1103/PhysRevD.69.123003

37. Aghanim N, Arnaud M, Ashdown M, Aumont J, Baccigalupi C, Banday AJ, et al. Planck 2015 results. XI. CMB power spectra, likelihoods, and robustness of parameters. Astron Astrophys. (2016) **594**:A11. doi: 10.1051/0004-6361/201526926

38. Aghanim N, Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. V. CMB power spectra and likelihoods. arXiv [Preprint]. arXiv:1907.12875 (2019).

39. Verde L, Peiris HV, Spergel DN, Nolta M, Bennett CL, Halpern M, et al. First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: parameter estimation methodology. Astrophys J Suppl. (2003) **148**:195. doi: 10.1086/377335

40. Ade PAR, Aghanim N, Armitage-Caplan C, Arnaud M, Ashdown M, Atrio-Barandela F, et al. Planck 2013 results. XV. CMB power spectra and likelihood. Astron Astrophys. (2014) **571**:A15. doi: 10.1015/0004-6361/201321573

41. Lewis A, Challinor A. Weak gravitational lensing of the CMB. Phys Rept. (2006) **429**:1-65. doi: 10.1016/j.physrep.2006.03.002

42. Ade PAR, Aghanim N, Armitage-Caplan C, Arnaud M, Ashdown M, Atrio-Barandela F, et al. Planck 2013 results. XVII. Gravitational
lensing by large-scale structure. *Astron Astrophys.* (2014) **571**:A17.
doi: 10.1051/0004-6361/201321543

43. Okamoto T, Hu W. CMB lensing reconstruction on the full sky. *Phys Rev D.* (2003) **67**:083002. doi: 10.1103/PhysRevD.67.083002

44. Hu W, Okamoto T. Mass reconstruction with cmb polarization. *Astrophys J.* (2002) **574**:566–74. doi: 10.1086/341110

45. Aghanim N, Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. VIII. Gravitational lensing. *arXiv [Preprint]*. arXiv:1807.06210 (2018).

46. Sherwin BD, van Engelen A, Sehgal N, Madhavacheril M, Addison GE, Aiola S, et al. Two-season Atacama Cosmology Telescope polarimeter lensing power spectrum. *Phys Rev D.* (2017) **95**:123529. doi: 10.1103/PhysRevD.95.123529

47. Omori Y, Chown R, Simard G, Story KT, Aylor K, Baxter EJ, et al. A 2500 deg² CMB lensing map from combined south pole telescope and planck data. *Astrophys J.* (2017) **849**:124. doi: 10.3847/1538-4357/aa8d1d

48. Ade PAR, Akiba Y, Anthony AE, Arnold K, Barron D, Boettger D, et al. Evidence for gravitational lensing of the cosmic microwave background polarization from cross-correlation with the cosmic infrared background. *Phys Rev Lett.* (2014) **112**:131302. doi: 10.1103/PhysRevLett.112.131302

49. Challinor A, Allison R, Carron J, Errard J, Feeney S, Kitching T, et al. Exploring cosmic origins with CORE: gravitational lensing of the CMB. *J Cosmol Astropart Phys.* (2018) **1804**:018. doi: 10.1088/1475-7516/2018/04/018

50. Schmittfull MM, Challinor A, Hanson D, Lewis A. Joint analysis of CMB temperature and lensing-reconstruction power spectra. *Phys Rev D.* (2013) **88**:063012. doi: 10.1103/PhysRevD.88.063012

51. Lewis A, Challinor A, Lasenby A. Efficient computation of CMB anisotropies in closed FRW models. *Astrophys J.* (2000) **538**:473–6. doi: 10.1086/309179

52. Lesgourgues J. The Cosmic Linear Anisotropy Solving System (CLASS) I: overview. *arXiv [Preprint]*. arXiv:1104.2932 (2011).

53. Benoit-Lévy A, Smith KM, Hu W. Non-Gaussian structure of the lensed CMB power spectra covariance matrix. *Phys Rev D.* (2012) **86**:123008. doi: 10.1103/PhysRevD.86.123008

54. Lewis A. Lensed CMB simulation and parameter estimation. *Phys Rev D.* (2005) **71**:083008. doi: 10.1103/PhysRevD.71.083008

55. Motloch P, Hu W. Lens covariance effects on likelihood analyses of CMB power spectra. *Phys Rev D.* (2017) **96**:103517. doi: 10.1103/PhysRevD.96.103517

56. Motloch P, Hu W, Benoit-Lévy A. CMB lens sample covariance and consistency relations. *Phys Rev D.* (2017) **95**:043518. doi: 10.1103/PhysRevD.95.043518

57. Green D, Meyers J, van Engelen A. CMB delensing beyond the B modes. *J Cosmol Astropart Phys.* (2017) **1712**:005. doi: 10.1088/1475-7516/2017/12/005

58. Vanneste S, Henrot-Versillé S, Louis T, Tristram M. Quadratic estimator for CMB cross-correlation. *Phys Rev D.* (2018) **98**:103526. doi: 10.1103/PhysRevD.98.103526

59. Aghanim N, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, Banday AJ, et al. Planck intermediate results. XLVI. Reduction of large-scale systematic effects in HFI polarization maps and estimation of the reionization optical depth. *Astron Astrophys.* (2016) **596**:A107. doi: 10.1051/0004-6361/201628890

60. Gorski KM, Hivon E, Banday AJ, Wandelt BD, Hansen FK, Reinecke M, et al. HEALPix - A Framework for high resolution discretization, and fast analysis of data distributed on the sphere. *Astrophys J.* (2005) **622**:759–71. doi: 10.1086/427976

61. Hinshaw G, Nolta MR, Bennett CL, Bean R, Dore O, Greason, MR, et al. Three-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: temperature analysis. *Astrophys J Suppl.* (2007) **170**:288. doi: 10.1086/513698

62. Eriksen HK, O’Dwyer IJ, Jewell JB, Wandelt BD, Larson DL, Gorski KM, et al. Power spectrum estimation from high-resolution maps by Gibbs sampling. *Astrophys J Suppl.* (2004) **155**:227–41. doi: 10.1086/425219

63. Jewell J, Levin S, Anderson CH. Application of Monte Carlo algorithms to the Bayesian analysis of the cosmic microwave background. *Astrophys J.* (2004) **609**:1–14. doi: 10.1086/383515

64. Wandelt BD, Larson DL, Lakshminarayanan A. Global, exact cosmic microwave background data analysis using Gibbs sampling. *Phys Rev D.* (2004) **70**:083511. doi: 10.1103/PhysRevD.70.083511

65. Eriksen HK, Jewell JB, Dickinson C, Banday AJ, Gorski KM, Lawrence CR. Joint Bayesian component separation and CMB power spectrum estimation. *Astrophys J.* (2008) **676**:10–32. doi: 10.1086/525277

66. Chu M, Eriksen HK, Knox L, Gorski KM, Jewell JB, Larson DL, et al. Cosmological parameter constraints as derived from the Wilkinson Microwave Anisotropy Probe data via Gibbs sampling and the Blackwell-Rao estimator. *Phys Rev D.* (2005) **71**:103002. doi: 10.1103/PhysRevD.71.103002

67. Rudjord O, Groeneboom NE, Eriksen HK, Huey G, Gorski KM, Jewell JB. CMB likelihood approximation by a Gaussianized Blackwell-Rao estimator. *Astrophys J.* (2009) **692**:1669–77. doi: 10.1088/0004-637X/692/2/1669

68. Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, Banday AJ, et al. Planck 2018 results. IV. Diffuse component separation. *arXiv [Preprint]*. arXiv:1807.06208 (2018).

69. Ade PAR, Aghanim N, Alves MIR, Arnaud M, Ashdown M, Aumont J, et al. Planck 2015 results. XXV. Diffuse low-frequency Galactic foregrounds. *Astron Astrophys.* (2016) **594**:A25. doi: 10.1051/0004-6361/201526803

70. Ade PAR, Aghanim N, Armitage-Caplan C, Arnaud M, Ashdown M, Atrio-Barandela F, et al. Planck 2013 results. XII. Diffuse component separation. *Astron Astrophys.* (2014) **571**: A12. doi: 10.1051/0004-6361/201321580

71. Larson D, Dunkley J, Hinshaw G, Komatsu E, Nolta MR, Bennett CL, et al. Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) observations: power spectra and WMAP-derived parameters. *Astrophys J Suppl.* (2011) **192**:16. doi: 10.1088/0067-0049/192/2/16

72. Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW. Population growth of human y chromosomes: a study of y chromosome microsatellites. *Mol Biol Evol.* (1999) **16**:1791–8.

73. Beaumont MA, Zhang W, Balding DJ. Approximate bayesian computation in population genetics. *Genetics.* (2002) **162**:2025–35.

74. Marjoram P, Molitor J, Plagnol V, Tavaré S. Markov chain monte carlo without likelihoods. *Proc Natl Acad Sci USA.* (2003) **100**:15324–8. doi: 10.1073/pnas.0306899100

75. Sisson SA, Fan Y, Tanaka MM. Sequential monte carlo without likelihoods. *Proc Natl Acad Sci USA.* (2007) **104**:1760–5. doi: 10.1073/pnas.0607208104

76. Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MP. Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems. *J R Soc Interface.* (2009) **6**:187–202. doi: 10.1098/rsif.2008.0172

77. Beaumont MA, Cornuët JM, Marin JM, Robert CP. Adaptive approximate bayesian computation. *Biometrika.* (2009) **96**:983–90. doi: 10.1093/biomet/asp052

78. McKinley T, Cook AR, Deardon R. Inference in epidemic models without likelihoods. *Int J Biostat.* **5**:24. doi: 10.2202/1557-4679.1171

79. Beaumont MA. Approximate bayesian computation in evolution and ecology. *Annu Rev Ecol Evol Syst.* (2010) **41**:379–406. doi: 10.1146/annurev-ecolsys-102209-144621

80. Liepe J, Kirk P, Filippi S, Toni T, Barnes CP, Stumpf MP. A framework for parameter estimation and model selection from experimental data in systems biology using approximate bayesian computation. *Nat Protoc.* (2014) **9**:439–56. doi: 10.1038/nprot.2014.025

81. Turner BM, Van Zandt T. A tutorial on approximate bayesian computation. *J Math Psychol.* (2012) **56**:69–85. doi: 10.1016/j.jmp.2012.02.005

82. Akeret J, Refregier A, Amara A, Seehars S, Hasner C, Approximate Bayesian computation for forward modeling in cosmology. *J Cosmol Astropart Phys.* (2015) **1508**:043. doi: 10.1088/1475-7516/2015/08/043

83. Alsing J, Wandelt B. Nuisance hardened data compression for fast likelihood-free inference. *Mon Not R Astron Soc.* (2019) **488**:5093. doi: 10.1093/mnras/stz1900

84. Alsing J, Charnock T, Feeney S, Wandelt B. Fast likelihood-free cosmology with neural density estimators and active learning. *Mon Not R Astron Soc.* (2019) **488**:4440. doi: 10.1093/mnras/stz1960

85. Leclercq F. Bayesian optimization for likelihood-free cosmological inference. *Phys Rev D.* (2018) **98**:063511. doi: 10.1103/PhysRevD.98.063511
lensing by large-scale structure. *Astron Astrophys.* (2014) **571**:A17.
doi: 10.1051/0004-6361/201321543

43. Okamoto T, Hu W. CMB lensing reconstruction on the full sky. *Phys Rev D.* (2003) **67**:083002. doi: 10.1103/PhysRevD.67.083002

44. Hu W, Okamoto T. Mass reconstruction with cmb polarization. *Astrophys J.* (2002) **574**:566–74. doi: 10.1086/341110

45. Aghanim N, Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, et al. Planck 2018 results. VIII. Gravitational lensing. arXiv [Preprint]. arXiv:1807.06210 (2018).

46. Sherwin BD, van Engelen A, Sehgal N, Madhavacheril M, Addison GE, Aiola S, et al. Two-season Atacama Cosmology Telescope polarimeter lensing power spectrum. Phys Rev D. (2017) 95:123529. doi: 10.1103/PhysRevD.95.123529

47. Omori Y, Chown R, Simard G, Story KT, Aylor K, Baxter EJ, et al. A 2500 deg² CMB lensing map from combined south pole telescope and planck data. Astrophys J. (2017) 849:124. doi: 10.3847/1538-4357/aa8d1d

48. Ade PAR, Akiba Y, Anthony AE, Arnold K, Barron D, Boettger D, et al. Evidence for gravitational lensing of the cosmic microwave background polarization from cross-correlation with the cosmic infrared background. Phys Rev Lett. (2014) 112:131302. doi: 10.1103/PhysRevLett.112.131302

49. Challinor A, Allison R, Carron J, Errard J, Feeney S, Kitching T, et al. Exploring cosmic origins with CORE: gravitational lensing of the CMB. J Cosmol Astropart Phys. (2018) 1804:018. doi: 10.1088/1475-7516/2018/04/018

50. Schmittfull MM, Challinor A, Hanson D, Lewis A. Joint analysis of CMB temperature and lensing-reconstruction power spectra. Phys Rev D. (2013) **88**:063012. doi: 10.1103/PhysRevD.88.063012

51. Lewis A, Challinor A, Lasenby A. Efficient computation of CMB anisotropies in closed FRW models. Astrophys J. (2000) 538:473–6. doi: 10.1086/309179

52. Lesgourgues J. The Cosmic Linear Anisotropy Solving System (CLASS) I: overview. arXiv [Preprint]. arXiv:1104.2932 (2011).

53. Benoit-Lévy A, Smith KM, Hu W. Non-Gaussian structure of the lensed CMB power spectra covariance matrix. Phys Rev D. (2012) 86:123008. doi: 10.1103/PhysRevD.86.123008

54. Lewis A. Lensed CMB simulation and parameter estimation. Phys Rev D. (2005) **71**:083008. doi: 10.1103/PhysRevD.71.083008

55. Motloch P, Hu W. Lens covariance effects on likelihood analyses of CMB power spectra. Phys Rev D. (2017) 96:103517. doi: 10.1103/PhysRevD.96.103517

56. Motloch P, Hu W, Benoit-Lévy A. CMB lens sample covariance and consistency relations. Phys Rev D. (2017) 95:043518. doi: 10.1103/PhysRevD.95.043518

57. Green D, Meyers J, van Engelen A. CMB delensing beyond the B modes. J Cosmol Astropart Phys. (2017) 1712:005. doi: 10.1088/1475-7516/2017/12/005

58. Vanneste S, Henrot-Versillé S, Louis T, Tristram M. Quadratic estimator for CMB cross-correlation. Phys Rev D. (2018) 98:103526. doi: 10.1103/PhysRevD.98.103526

59. Aghanim N, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, Banday AJ, et al. Planck intermediate results. XLVI. Reduction of large-scale systematic effects in HFI polarization maps and estimation of the reionization optical depth. Astron Astrophys. (2016) 596:A107. doi: 10.1051/0004-6361/201628890

60. Gorski KM, Hivon E, Banday AJ, Wandelt BD, Hansen FK, Reinecke M, et al. HEALPix - A Framework for high resolution discretization, and fast analysis of data distributed on the sphere. Astrophys J. (2005) 622:759–71. doi: 10.1086/427976

61. Hinshaw G, Nolta MR, Bennett CL, Bean R, Dore O, Greason, MR, et al. Three-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: temperature analysis. Astrophys J Suppl. (2007) 170:288. doi: 10.1086/513698

62. Eriksen HK, O'Dwyer IJ, Jewell JB, Wandelt BD, Larson DL, Gorski KM, et al. Power spectrum estimation from high-resolution maps by Gibbs sampling. Astrophys J Suppl. (2004) 155:227–41. doi: 10.1086/425219

63. Jewell J, Levin S, Anderson CH. Application of Monte Carlo algorithms to the Bayesian analysis of the cosmic microwave background. Astrophys J. (2004) 609:1–14. doi: 10.1086/383515

64. Wandelt BD, Larson DL, Lakshminarayanan A. Global, exact cosmic microwave background data analysis using Gibbs sampling. Phys Rev D. (2004) 70:083511. doi: 10.1103/PhysRevD.70.083511

65. Eriksen HK, Jewell JB, Dickinson C, Banday AJ, Gorski KM, Lawrence CR. Joint Bayesian component separation and CMB power spectrum estimation. Astrophys J. (2008) 676:10–32. doi: 10.1086/525277

66. Chu M, Eriksen HK, Knox L, Gorski KM, Jewell JB, Larson DL, et al. Cosmological parameter constraints as derived from the Wilkinson Microwave Anisotropy Probe data via Gibbs sampling and the Blackwell-Rao estimator. Phys Rev D. (2005) 71:103002. doi: 10.1103/PhysRevD.71.103002

67. Rudjord O, Groeneboom NE, Eriksen HK, Huey G, Gorski KM, Jewell JB. CMB likelihood approximation by a Gaussianized Blackwell-Rao estimator. Astrophys J. (2009) 692:1669–77. doi: 10.1088/0004-637X/692/2/1669

68. Akrami Y, Ashdown M, Aumont J, Baccigalupi C, Ballardini M, Banday AJ, et al. Planck 2018 results. IV. Diffuse component separation. arXiv [Preprint]. arXiv:1807.06208 (2018).

69. Ade PAR, Aghanim N, Alves MIR, Arnaud M, Ashdown M, Aumont J, et al. Planck 2015 results. XXV. Diffuse low-frequency Galactic foregrounds. Astron Astrophys. (2016) 594:A25. doi: 10.1051/0004-6361/201526803

70. Ade PAR, Aghanim N, Armitage-Caplan C, Arnaud M, Ashdown M, Atrio-Barandela F, et al. Planck 2013 results. XII. Diffuse component separation. Astron Astrophys. (2014) 571: A12. doi: 10.1051/0004-6361/201321580

71. Larson D, Dunkley J, Hinshaw G, Komatsu E, Nolta MR, Bennett CL, et al. Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) observations: power spectra and WMAP-derived parameters. Astrophys J Suppl. (2011) 192:16. doi: 10.1088/0067-0049/192/2/16

72. Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW. Population growth of human y chromosomes: a study of y chromosome microsatellites. Mol Biol Evol. (1999) 16:1791–8.

73. Beaumont MA, Zhang W, Balding DJ. Approximate bayesian computation in population genetics. Genetics. (2002) 162:2025–35.

74. Marjoram P, Molitor J, Plagnol V, Tavaré S. Markov chain monte carlo without likelihoods. Proc Natl Acad Sci USA. (2003) 100:15324–8. doi: 10.1073/pnas.0306899100

75. Sisson SA, Fan Y, Tanaka MM. Sequential monte carlo without likelihoods. Proc Natl Acad Sci USA. (2007) 104:1760–5. doi: 10.1073/pnas.0607208104

76. Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MP. Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems. J R Soc Interface. (2009) 6:187–202. doi: 10.1098/rsif.2008.0172

77. Beaumont MA, Cornuet JM, Marin JM, Robert CP. Adaptive approximate bayesian computation. Biometrika. (2009) 96:983–90. doi: 10.1093/biomet/asp052

78. McKinley T, Cook AR, Deardon R. Inference in epidemic models without likelihoods. Int J Biostat. 5:24. doi: 10.2202/1557-4679.1171

79. Beaumont MA. Approximate bayesian computation in evolution and ecology. Annu Rev Ecol Evol Syst. (2010) 41:379–406. doi: 10.1146/annurev-ecolsys-102209-144621

80. Liepe J, Kirk P, Filippi S, Toni T, Barnes CP, Stumpf MP. A framework for parameter estimation and model selection from experimental data in systems biology using approximate bayesian computation. Nat Protoc. (2014) 9:439–56. doi: 10.1038/nprot.2014.025

81. Turner BM, Van Zandt T. A tutorial on approximate bayesian computation. J Math Psychol. (2012) 56:69–85. doi: 10.1016/j.jmp.2012.02.005

82. Akeret J, Refregier A, Amara A, Seehars S, Hasner C, Approximate Bayesian computation for forward modeling in cosmology. J Cosmol Astropart Phys. (2015) 1508:043. doi: 10.1088/1475-7516/2015/08/043

83. Alsing J, Wandelt B. Nuisance hardened data compression for fast likelihood-free inference. Mon Not R Astron Soc. (2019) 488:5093. doi: 10.1093/mnras/stz1900

84. Alsing J, Charnock T, Feeney S, Wandelt B. Fast likelihood-free cosmology with neural density estimators and active learning. Mon Not R Astron Soc. (2019) 488:4440. doi: 10.1093/mnras/stz1960

85. Leclercq F. Bayesian optimization for likelihood-free cosmological inference. Phys Rev D. (2018) 98:063511. doi: 10.1103/PhysRevD.98.063511
86. Alsing J, Wandelt B, Feeney S. Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology. *Mon Not R Astron Soc.* (2018) **477**:2874. doi: 10.1093/mnras/sty819

87. Delabrouille J, Betoule M, Melin JB, Miville-Deschenes MA, Gonzalez-Nuevo J, Le Jeune M, et al. The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths. *Astron. Astrophys.* (2013) **553**:A96. doi: 10.1051/0004-6361/201220019

88. Sehgal N, Bode P, Das S, Hernandez-Monteagudo C, Huffenberger K, Lin Y, et al. Simulations of the microwave sky. *Astrophys J.* (2010) **709**:920–36. doi: 10.1088/0004-637X/709/2/920

89. Ichiki K. CMB foreground: a concise review. *Prog Theor Exp Phys.* **6**:2014. doi: 10.1093/ptep/ptu065

90. Takakura S, Aguilar M, Akiba Y, Arnold K, Baccigalupi C, Barron D, et al. Performance of a continuously rotating half-wave plate on the POLARBEAR telescope *J Cosmol Astropart Phys.* (2017) **1705**:008. doi: 10.1088/1475-7516/2017/05/008

91. Leach SM, Cardoso JF, Baccigalupi C, Barreiro RB, Betoule M, Bobin J, et al. Component separation methods for the Planck mission. *Astron Astrophys.* (2008) **491**:597–615. doi: 10.1051/0004-6361:200810116

92. Bennett CL, Smoot GF, Hinshaw G, Wright EL, Kogut A, de Amici G, et al. Preliminary separation of galactic and cosmic microwave emission for the COBE Differential Microwave Radiometer. *Astrophys J.* (1992) **396**:L7–12. doi: 10.1086/186505

93. Fernandez-Cobos R, Vielva P, Barreiro RB, Martinez-Gonzalez E. Multi-resolution internal template cleaning: an application to the Wilkinson Microwave Anisotropy Probe 7-yr polarization data. *Mon Not R Astron Soc.* (2012) **420**:2162. doi: 10.1111/j.1365-2966.2011.20182.x

94. Delabrouille J, Cardoso JF, Jeune ML, Betoule M, Fay G, Guilloux F. A full sky, low foreground, high resolution CMB map from WMAP. *Astron Astrophys.* (2009) **493**:835. doi: 10.1051/0004-6361:200810514

95. Eriksen HK, Dickinson C, Lawrence CR, Baccigalupi C, Banday AJ, Gorski KM, et al. CMB component separation by parameter estimation. *Astrophys J.* (2006) **641**:665–82. doi: 10.1086/500499

96. Cardoso JF, Martin M, Delabrouille J, Betoule M, Patanchon G. Component separation with flexible models. Application to the separation of astrophysical emissions. *arXiv [Preprint]*. arXiv:0803.1814 (2008).

97. Dunkley J, Calabrese E, Sievers J, Addison GE, Battaglia N, Battistelli ES, et al. The Atacama Cosmology Telescope: likelihood for small-scale CMB data" *J Cosmol Astropart Phys.* (2013) **1307**:025. doi: 10.1088/1475-7516/2013/07/025

98. Keisler R, Reichardt CL, Aird KA, Benson BA, Bleem LE, Carlstrom JE, et al. A measurement of the damping tail of the cosmic microwave background power spectrum with the south pole telescope. *Astrophys J.* (2011) **743**:28. doi: 10.1088/0004-637X/743/1/28

99. Suzuki A, Ade PAR, Akiba Y, Alonso D, Arnold K, Aumont J, et al. The LiteBIRD Satellite Mission - Sub-Kelvin Instrument. *J Low Temp Phys.* (2018) **193**:1048–56. doi: 10.1007/s10909-018-1947-7

100. Benabed K, Cardoso JF, Prunet S, Hivon E. TEASING: a fast and accurate approximation for the low multipole likelihood of the Cosmic Microwave Background temperature. *Mon Not R Astron Soc.* (2009) **400**:219. doi: 10.1111/j.1365-2966.2009.15202.x

101. Page L, Hinshaw G, Komatsu E, Nolta MR, Spergel DN, Bennett CL, et al. Three year Wilkinson Microwave Anisotropy Probe (WMAP) observations: polarization analysis. *Astrophys J Suppl.* (2007) **170**:335. doi: 10.1086/513699

102. Hivon E, Gorski KM, Netterfield CB, Crill BP, Prunet S, Hansen F. Master of the cosmic microwave background anisotropy power spectrum: a fast method for statistical analysis of large and complex cosmic microwave background data sets. *Astrophys J.* (2002) **567**:2. doi: 10.1086/338126

103. Wandelt BD, Hivon E, Gorski KM. The pseudo-$C_\ell$ method: cosmic microwave background anisotropy power spectrum statistics for high precision cosmology. *Phys Rev D.* (2001) **64**:083003. doi: 10.1103/PhysRevD.64.083003

104. Brown ML, Castro PG, Taylor AN. Cosmic microwave background temperature and polarization pseudo-$C_\ell$ estimators and covariances. *Mon Not R Astron Soc.* (2005) **360**:4. doi: 10.1111/j.1365-2966.2005.09111.x

105. Bunn EF, Zaldarriaga M, Tegmark M, de Oliveira-Costa A. E/B decomposition of finite pixelized CMB maps. *Phys Rev D.* (2003) **67**:023501. doi: 10.1103/PhysRevD.67.023501

106. Challinor A, Chon G. Error analysis of quadratic power spectrum estimates for CMB polarization: sampling covariance. *Mon Not R Astron Soc.* (2005) **360**:509. doi: 10.1111/j.1365-2966.2005.09076.x

107. Grain J, Tristram M, Stompor R. Polarized CMB power spectrum estimation using the pure pseudo-cross-spectrum approach. *Phys Rev D.* (2009) **79**:123515. doi: 10.1103/PhysRevD.79.123515

108. Smith KM, Zaldarriaga M. A general solution to the E-B mixing problem. *Phys Rev D.* (2007) **76**:043001. doi: 10.1103/PhysRevD.76.043001

109. Smith KM. Pseudo-$C_\ell$ estimators which do not mix E and B modes. *Phys Rev D.* (2006) **74**:083002. doi: 10.1103/PhysRevD.74.083002

110. Zaldarriaga M. Fluctuations in the cosmic microwave background. *arXiv:astro-ph/9806122* (1998).

111. Tegmark M. How to measure CMB power spectra without losing information. *Phys Rev D.* (1997)55:5895. doi: 10.1103/PhysRevD.55.5895

112. Bond JR, Jaffe AH, Knox L. Estimating the power spectrum of the cosmic microwave background. *Phys Rev D.* (1998) **57**:2117. doi: 10.1103/PhysRevD.57.2117

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Gerbino, Lattanzi, Migliaccio, Pagano, Salvati, Colombo, Gruppuso, Natoli and Polenta. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
86. Alsing J, Wandelt B, Feeney S. Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology. *Mon Not R Astron Soc.* (2018) **477**:2874. doi: 10.1093/mnras/sty819

87. Delabrouille J, Betoule M, Melin JB, Miville-Deschenes MA, Gonzalez-Nuevo J, Le Jeune M, et al. The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths. *Astron. Astrophys.* (2013) **553**:A96. doi: 10.1051/0004-6361/201220019

88. Sehgal N, Bode P, Das S, Hernandez-Monteagudo C, Huffenberger K, Lin Y, et al. Simulations of the microwave sky. *Astrophys J.* (2010) **709**:920–36. doi: 10.1088/0004-637X/709/2/920

89. Ichiki K. CMB foreground: a concise review. *Prog Theor Exp Phys.* **6**:2014. doi: 10.1093/ptep/ptu065

90. Takakura S, Aguilar M, Akiba Y, Arnold K, Baccigalupi C, Barron D, et al. Performance of a continuously rotating half-wave plate on the POLARBEAR telescope *J Cosmol Astropart Phys.* (2017) **1705**:008. doi: 10.1088/1475-7516/2017/05/008

91. Leach SM, Cardoso JF, Baccigalupi C, Barreiro RB, Betoule M, Bobin J, et al. Component separation methods for the Planck mission. *Astron Astrophys.* (2008) **491**:597–615. doi: 10.1051/0004-6361:200810116

92. Bennett CL, Smoot GF, Hinshaw G, Wright EL, Kogut A, de Amici G, et al. Preliminary separation of galactic and cosmic microwave emission for the COBE Differential Microwave Radiometer. *Astrophys J.* (1992) **396**:L7–12. doi: 10.1086/186505

93. Fernandez-Cobos R, Vielva P, Barreiro RB, Martinez-Gonzalez E. Multi-resolution internal template cleaning: an application to the Wilkinson Microwave Anisotropy Probe 7-yr polarization data. *Mon Not R Astron Soc.* (2012) **420**:2162. doi: 10.1111/j.1365-2966.2011.20182.x

94. Delabrouille J, Cardoso JF, Jeune ML, Betoule M, Fay G, Guilloux F. A full sky, low foreground, high resolution CMB map from WMAP. *Astron Astrophys.* (2009) **493**:835. doi: 10.1051/0004-6361:200810514

95. Eriksen HK, Dickinson C, Lawrence CR, Baccigalupi C, Banday AJ, Gorski KM, et al. CMB component separation by parameter estimation. *Astrophys J.* (2006) **641**:665–82. doi: 10.1086/500499

96. Cardoso JF, Martin M, Delabrouille J, Betoule M, Patanchon G. Component separation with flexible models. Application to the separation of astrophysical emissions. *arXiv [Preprint]*. arXiv:0803.1814 (2008).

97. Dunkley J, Calabrese E, Sievers J, Addison GE, Battaglia N, Battistelli ES, et al. The Atacama Cosmology Telescope: likelihood for small-scale CMB data" *J Cosmol Astropart Phys.* (2013) **1307**:025. doi: 10.1088/1475-7516/2013/07/025

98. Keisler R, Reichardt CL, Aird KA, Benson BA, Bleem LE, Carlstrom JE, et al. A measurement of the damping tail of the cosmic microwave background power spectrum with the south pole telescope. *Astrophys J.* (2011) **743**:28. doi: 10.1088/0004-637X/743/1/28

99. Suzuki A, Ade PAR, Akiba Y, Alonso D, Arnold K, Aumont J, et al. The LiteBIRD Satellite Mission - Sub-Kelvin Instrument. *J Low Temp Phys.* (2018) **193**:1048–56. doi: 10.1007/s10909-018-1947-7

100. Benabed K, Cardoso JF, Prunet S, Hivon E. TEASING: a fast and accurate approximation for the low multipole likelihood of the Cosmic Microwave Background temperature. *Mon Not R Astron Soc.* (2009) **400**:219. doi: 10.1111/j.1365-2966.2009.15202.x

101. Page L, Hinshaw G, Komatsu E, Nolta MR, Spergel DN, Bennett CL, et al. Three year Wilkinson Microwave Anisotropy Probe (WMAP) observations: polarization analysis. *Astrophys J Suppl.* (2007) **170**:335. doi: 10.1086/513699

102. Hivon E, Gorski KM, Netterfield CB, Crill BP, Prunet S, Hansen F. Master of the cosmic microwave background anisotropy power spectrum: a fast method for statistical analysis of large and complex cosmic microwave background data sets. *Astrophys J.* (2002) **567**:2. doi: 10.1086/338126

103. Wandelt BD, Hivon E, Gorski KM. The pseudo-$C_\ell$ method: cosmic microwave background anisotropy power spectrum statistics for high precision cosmology. *Phys Rev D.* (2001) **64**:083003. doi: 10.1103/PhysRevD.64.083003

104. Brown ML, Castro PG, Taylor AN. Cosmic microwave background temperature and polarization pseudo-$C_\ell$ estimators and covariances. *Mon Not R Astron Soc.* (2005) **360**:4. doi: 10.1111/j.1365-2966.2005.09111.x

105. Bunn EF, Zaldarriaga M, Tegmark M, de Oliveira-Costa A. E/B decomposition of finite pixelized CMB maps. *Phys Rev D.* (2003) **67**:023501. doi: 10.1103/PhysRevD.67.023501

106. Challinor A, Chon G. Error analysis of quadratic power spectrum estimates for CMB polarization: sampling covariance. *Mon Not R Astron Soc.* (2005) **360**:509. doi: 10.1111/j.1365-2966.2005.09076.x

107. Grain J, Tristram M, Stompor R. Polarized CMB power spectrum estimation using the pure pseudo-cross-spectrum approach. *Phys Rev D.* (2009) **79**:123515. doi: 10.1103/PhysRevD.79.123515

108. Smith KM, Zaldarriaga M. A general solution to the E-B mixing problem. *Phys Rev D.* (2007) **76**:043001. doi: 10.1103/PhysRevD.76.043001

109. Smith KM. Pseudo-$C_\ell$ estimators which do not mix E and B modes. *Phys Rev D.* (2006) **74**:083002. doi: 10.1103/PhysRevD.74.083002

110. Zaldarriaga M. Fluctuations in the cosmic microwave background. *arXiv:astro-ph/9806122* (1998).

111. Tegmark M. How to measure CMB power spectra without losing information. *Phys Rev D.* (1997)55:5895. doi: 10.1103/PhysRevD.55.5895

112. Bond JR, Jaffe AH, Knox L. Estimating the power spectrum of the cosmic microwave background. *Phys Rev D.* (1998) **57**:2117. doi: 10.1103/PhysRevD.57.2117

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Gerbino, Lattanzi, Migliaccio, Pagano, Salvati, Colombo, Gruppuso, Natoli and Polenta. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.