Maximum likelihood search results for Data Challenge 02

Summary

This is a long posting, so here are the important points:

• A significant part of the bias that we see in this maximum-likelihood search analysis scales with \(A_L\) and seems to just be a mismatch between the lensing \(BB\) spectrum used for the sims and the one used in the analysis. But this also reinforces the idea that we should marginalize over some (small) uncertainty in the lensing spectrum.
• Some of the foreground models, especially Hensley/Draine dust and Ghosh dust, yield low values of \(\beta_d\) which is in tension with the Planck-derived prior used in the analysis. It would be a good idea to rerun this analysis without that prior, since CMB-S4 is quite capable of measuring \(\beta_d\).
• We should try some Fisher forecasts that assume some level of foreground decorrelation across frequencies. While the existing forecasts do include the decorrelation parameter as part of the model, they don't account for the effect of foreground decorrelation in the bandpower covariance matrix.

Introduction

This posting summarizes results from analysis of CMB-S4 Data Challenge 02 using a BICEP/Keck-style parametrized foreground model. For each realization, we find the set of model parameters that maximizes the likelihood multiplied by priors on the dust and sync spectral index parameters (\(\beta_d\) and \(\beta_s\)). These priors are based on Planck data, so they are quite weak in comparison with CMB-S4 sensitivity.

The model includes the following parameters:

• \(r\): tensor-to-scalar ratio
• \(A_d\): \(BB\) power spectrum amplitude of dust, in \(\mu K_{CMB}^2\) units at \(\nu\)=353 GHz and \(\ell\)=80
• \(\beta_d\): dust emissivity spectral index; Gaussian prior on \(\beta_d\) is centered at 1.6 and has width 0.11
• \(A_s\): \(BB\) power spectrum amplitude of synchrotron, in \(\mu K_{CMB}^2\) units at \(\nu\)=23 GHz and \(\ell\)=80
• \(\beta_s\): synchrotron emissivity spectral index; Gaussian prior on \(\beta_s\) is centered at -3.1 and has width 0.3
• \(\alpha_d\): power law index for dust \(\mathcal{D}_\ell\) scaling in \(\ell\); limited to range [-2,2]
• \(\alpha_s\): power law index for synchrotron \(\mathcal{D}_\ell\) scaling in \(\ell\); limited to range [-2,2]
• \(\epsilon\): frequency-independent spatial correlation of dust and synchrotron; limited to range [-1,1]
• \(\rho_d\): dust correlation between 217 and 353 GHz; not included when “Decorrelation model” is set to “none”

For the decorrelation model, we assume that the cross-spectrum of dust between frequencies \(\nu_1\) and \(\nu_2\) is reduced by factor \(\exp\{ - (1-\rho) \times [\log^2(\nu_1 / \nu_2) / \log^2(217 / 353)] \times f(\ell)\}\). For the \(\ell\) dependence, we try out uniform (flat), linear, or quadratic scaling (pivot scale is \(\ell\)=80).

The Figure 1 pager shows histograms of the recovered maximum likelihood parameters for all sky models, four delensing options (residual \(A_L\) = 1.0, 0.3, 0.1, or 0.03), two choices of input tensor signal (\(r\) = 0 or 0.003), and four options for the dust decorrelation model used in the analysis. Here I will highlight the main results for each sky model.

00: Gaussian foregrounds

The mean values and standard deviations of \(r\) for simulations with simple Gaussian foregrounds are summarized in Table 00. With a 10% lensing residual, we don't quite achieve \(\sigma(r) = 5 \times 10^{-3}\) for sims with \(r = 0\).

We see a large (\(\sim \sigma / 2\)) negative bias for the \(A_L\) = 1 case, but the bias scales proportionally with \(A_L\), which suggests that the problem is due to a slightly different lensing \(BB\) spectrum being used for the sims vs the analysis. I thought that we already checked that these matched, but perhaps didn't check at high enough precision; should be an easy fix.

Turning on dust decorrelation in the model doesn't cause any bias in \(r\) and the recovered \(\rho_d\) values are centered around 1 (i.e. analysis recovers zero decorrelation). Adding this parameter does increase \(\sigma(r)\) somewhat. The hit is largest for flat scaling in \(\ell\) and smallest for quadratic scaling, because steeper scaling in \(\ell\) means that high \(\ell\) bandpowers can constrain the amount of dust decorrelation at \(\ell \sim 80\). This seems fairly artificial, so we shouldn't put too much stock in one version or another of the decorrelation model.

01: PySM a1d1s1f1

As has been previously noted, dust power is much higher in this model (\(A_d \sim 12.5 \mu K^2\)) than for the Gaussian foreground sims (\(A_d = 3.75 \mu K^2\)). The PySM d1 dust model does feature a spatially varying spectral index, but we don't find any detectable decorrelation in this analysis. The PySM s1 synchrotron model yields \(A_s \sim 0.5 \mu K^2\) and there is \(\sim 6\)% correlation between dust and sync.

We still see a (negative) bias on \(r\) that is proportional to the residual lensing power but there is an additional positive bias on \(r\), presumably due to foregrounds. Judging by the \(A_L = 0.03\) results, this foreground bias corresponds to \(r \sim 5 \times 10^{-4}\). This bias perhaps goes away when we include dust decorrelation with flat \(\ell\) scaling, and we see \(\sim 1 \sigma\) evidence for dust decorrelation in that case (\(\rho_d = 0.9995 \pm 0.0005\)).

02: PySM a2d4f1s3

The d4 version of PySM dust adds a second dust component (with different blackbody temperature and emissivity power law) based on Meisner & Finkbeiner (2014). I'm not sure what type of \(\beta_d\) spatial variations are included in this model, but I think it is more or less the same as for d1. The s3 synchrotron model adds curvature to the synchrotron spectral index: \(\beta_s \rightarrow \beta_s + C \ln (\nu / \nu_C)\). The a2 AME model uses a 2% polarization fraction for AME, which seems very high, but there is no attempt to model AME in this analysis.

Results for this model show that \(A_d\) is even larger (\(\sim 32.5 \mu K^2\)) than for the d1 dust model. The mean value of \(\beta_d\) decreases from 1.59 (for PySM d1 model) to 1.55, which is probably a sign of the two component dust. The mean value of \(\beta_s\) decreases from -3.05 (for PySM s1 model) to -3.13, which is probably due to synchrotron spectral curvature (and perhaps polarized AME?). Dust–sync correlation is higher, at \(\sim 10\)%, which could be from polarized AME.

The pattern of \(r\) positive biases is basically the same as for the PySM a1d1f1s1 model, including the observation that the bias is minimized for a model that includes decorrelation with flat \(\ell\) scaling. However in this case the significance of \(\rho_d \ne 1\) is down to only \(\sim 0.6 \sigma\).

03: PySM a2d7f1s3

The next PySM version uses the Hensley/Draine dust model, which has additional complexity in the dust SED (perhaps described in arXiv:1709.07897?). The level of dust power is similar to sky model 01 (PySM d1 model), but we find that the emissivity power law is even flatter than the last case, with \(\beta_d \sim 1.44\).

If we focus on the \(A_L = 0.03\) column, we find that the decorrelation model with flat \(\ell\) scaling produces a negative bias in this case and we find the smallest bias for decorrelation with linear \(\ell\) scaling. The significance of \(\rho_d \ne 1\) is a little more than \(1 \sigma\) for both of these options.

04: Ghosh dust model

The Ghosh dust model (described here) is based on GASS HI data with a model for the Galactic magnetic field. For these sims, it is combined with the PySM a2, f1, and s3 components (same as the two previous models).

The analysis of this model yields smaller still values of \(\beta_d \sim 1.3-1.4\) and makes me pretty worried about the Planck-derived \(\beta_d\) prior that was used for this analysis—these values are disfavored by the prior at \(\sim 2.5 \sigma\)! While CMB-S4 bandpowers are so sensitive that I think the likelihood should dominate over the \(\beta_d\) prior, it seems very possible that this level of disagreement between sims and prior could push results around (leading to bias on \(r\)). Dust-sync correlation is still present, but smaller (2–3%), which is probably due to the fact that the Ghosh dust sims don't know anything about the PySM synchrotron or AME components. The fact that they are correlated at all probably happens because both models are based on data at larger scales.

Dust decorrelation is small in absolute terms, but detected at high significance. Using a model without dust decorrelation leads to a large positive bias on \(r\) in the range \(4-5 \times 10^{-3}\). Dust decorrelation with linear \(\ell\) scaling produces the smallest biases, but still quite large compared to other sky models. Also note that the recovered \(\beta_d\) parameters shift around depending on the decorrelation \(\ell\) scaling, which doesn't help my confidence in these results.

05: Gaussian decorrelated dust model

I only had 100 realizations (50 each for \(r\) = 0, 0.003) to analyze for this Gaussian decorrelated model, so the statistics aren't great. However, this model has extremely large dust decorrelation (15% between 217 and 353 GHz at \(\ell\) = 80) and it exactly follows the assumed functional form of decorrelation with linear \(\ell\) scaling, so we can still draw some useful conclusions.

When we choose decorrelation with linear \(\ell\) scaling to match the sims, then we find no bias on \(r\) and recover \(\rho_d\) = 0.85. If we use the wrong \(\ell\) scaling for decorrelation, then we find positive biases on \(r\) in the range \(1-3 \times 10^{-3}\). I didn't bother including in the table the results with no decorrelation model because they are crazy.

An important point to note from this model is that, even for the unbiased case where the decorrelation is correctly modeled in both \(\nu\) and \(\ell\), we find \(\sigma(r) \sim 1.4\), much larger than the target sensitivity of CMB-S4. This shows that, for extreme levels of foreground decorrelation, we lose the ability to clean foregrounds from the maps because the foreground modes are significantly independent between the various CMB-S4 frequencies. Regardless of whether you are doing map-based cleaning or fitting the power spectra as we do here, the only way to improve sensitivity would be use more observing bands that are more closely spaced. It also makes the point that our Fisher forecasts should assume some non-zero level of decorrelation. Adding decorrelation as a free parameter to a forecast that assumes \(\rho_d = 1\) only captures part of the statistical penalty.

06: Flauger MHD foregrounds

My understanding is that this model uses MHD simulations to consistently model polarized dust and synchrotron in the Galactic magnetic field. This makes it quite interesting that this analysis finds negative dust-sync correlation with \(\epsilon \sim -0.34\). The dust power is similar to the Gaussian sims, and \(\beta_d\) matches the Planck value of 1.59. This analysis finds a synchrotron SED power law that is much flatter than usual, \(\beta_s \sim -2.6\), which is inconsistent with the prior at about \(1.5 \sigma\).

This model does not show any significant dust decorrelation and we recover \(r\) values near the sim input level for all choices of decorrelation model. In general, the results for this model look nearly as good as the simple Gaussian foregrounds (sky model 00).

Appendix A: Results for CDT report

As can be seen in Table 00, we seem to have a negative bias even in the case of the Gaussian foreground simulations. This is being investigated, by changing the lensing power spectrum, and removing the priors on \(\beta_d\), and \(\beta_s\), as well as extending the allowed range of \(A_d\) for some models. Work is also on-going to see if these biases can be tracked down to specific bandpowers.

For the CDT report, it was decided to remove this "algorithmic bias" to focus on the bias produced by the different dust simulations. We also chose to report results using the linear \(\ell\) dependence for the decorrelation model. See caption of Table 07.

Appendix B: DC3-DC2 biases from the ML search (cross check with CDT reported numbers)

Here again, we only present the results for \(A_L\) = 0.1, but both for the linear decorrelation and no decorrelation case.

The numbers used in the current draft of the CDT report to show the impact on r recovery of introducing systematics into sims (table 9) are based on the method used by Steve (see e.g. here in a post using sky model 0). However, the other results on r from the parametric analyses shown in the report come from ML searches, as in this posting, so we wish to check whether these methods produce consistent numbers.

The results are obtained by subtracting DC2 results from the DC3, for sky model 3. Note that for DC3, the 1000 simulations were divided in subsamples to which specific systematics were added (see experiment definition), except for the first 125 (case "none"), which are thus the same as in DC2.

For instance, for sky model 3, we take all the even simulations within (0125-0249), which have a fiducial \(r\) = 0.003, and run the analysis both for DC3 and DC2. We then take the difference of the means. These results are reported in Table 08. We highlighted what corresponds to CDT's table 9, i.e. \(r\)=0, for linear decorrelation.