Maximum likelihood search results for Data Challenge 04, fixed bandpasses

Introduction

This posting summarizes results from analysis of CMB-S4 Data Challenge 04 using a BICEP/Keck-style parametrized foreground model.
There are a few additional changes compared to the DC02 analysis:
• we fixed the lensing spectrum to the input one. In the previous posting it uses a lensing spectrum that differs from the input,
• we zero out the theory spectrum below ell=30, since that's what was done in the sims.

The current posting is an update over this previous DC04 posting, which we keep separate for comparison.
As in the previous posting, we now use power spectra that have been re-generated to fix a bug, as shown in this posting.

In section 1, we show the main results in the form of figures and histograms, as well as triangle plots including foreground parameters. These rely on new cuts and changes explained in the following sections.
In section 2, we introduce a major improvement for model 00 from the rejection of a few flawed simulations.
In section 3, we show how the use of proper conventions for the banpasses seem to give more correct constraints, especially on the foreground parameters.
These 2 changes concern model 0, 5 and 6.
In section 4, we report tables of r constraints for the different sky models, for different lensing residuals, with and without decorrelation in the ML search.

To-do: In the next run, we will probably release the priors on the foreground parameters.

Note about the model:
In this analysis, 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. However, in principle foreground models may violate them potentially leading to biases (e.g. model 03 where the preferred value of \(\beta_d\) is outside the prior range - see below, Figure 2).

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]
• \(\Delta_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\{log(\Delta_d) \times [\log^2(\nu_1 / \nu_2) / \log^2(217 / 353)] \times f(\ell)\}\). For the \(\ell\) dependence we fix the scaling to take a linear form (pivot scale is \(\ell\)=80).

Note about the simulations:

• Model 0 and 5 provide Gaussian realizations of dust and synchrotron with power law spectrum. These are "uniform" over the sky.
• Model 1, 2, 3, and 6 have a more realistic variation over the sky, but we use only a single realization of each.
• The same if true for model 4 - we use only a single realization.

1: Summary Plots

In Figure 1, we summarize the r results using the new L-cut and bandpass correction fix. The changes are quite minimal, except for a decrease of \(\sigma(r)\) for model 00, which ranges from 3% to 20% depending on the lensing residual level and decorrelation models, as reported in this table.
As we will see in section 3 below the improvement is more important for the foreground parameters.

Note that for model 4, \(\beta_d\) is in slight tension with the imposed prior: Gaussian centered at 1.6 with width 0.11. Similarly for \(\beta_s\) in model 6, (compared to the Gaussian centered at -3.1 with width 0.3).

2: Flawed simulations removal

While analysing the distribution of the ML parameters for the 1000 "04.00" simulations, we found a few outliers. It seems like for some realizations have a flawed foreground component. At least some of these bad realizations have missing synchrotron in 85/95GHz as can be seen in the following pager link.
One way to filter out those realization is to apply a cut based on the goodness of fit. Here we use the log-likelihood (H-L likelihood) value and detect outliers using the modified Z-score method (with a threshold of 5, see here).
In Figure 3, we show the 1000 -log(L) values for sky models 00 to 06. Only model 00 has obvious outliers.
We can also see that somehow the first 100 simulations of model 5 are different from the 900 following. It seems to be due to a different decorrelation parameter value (see this figure).
The effect of these cuts can be seen in Figure 4. The main effect is a slight reduction of the error bars.
Caveat: note that this method relies on the fact that a flawed simulation results in a bad fit, which is not necessarily true.

3: Bandpass conventions

We also found that in the previous XX.00 runs, the level of foregrounds was significantly off in the ML search, as we see in Figure 2, for example in this case. The mean of the recovered ML amplitude parameter is shifted compared to the fiducial model used in the simulations by 20 to 30 \(\sigma\).
We realized that there was a missmatch in the bandpass conventions used in the Maximum Likelihood search. As is states in this definition page:
"What a tophat bandpass actually means is a little tricky and discussed in Bandpass Convention - What does flat mean. Models 00/05/06/09 use nu^-2 in spectral radiance (SR) units, while models 01/02/03 use nu^0 in SR units. (What 04 uses is not known.) "
Until now, the \(\nu^0\) convention was used for all models.
We ran some new ML searches where we updated to the proper convention for models 0, 5 and 6.

Comments on the figure:

• Fixing the bandpass convention seems to improve the recovery of \(A_d\) and \(A_s\), which go from this distribution to this distribution . There is still a \(\gt 3 \sigma\) deviation both for the "r=0" case and the "r=0.003" case.
• In the above cases, while the bias on the amplitude recovery is largely lowered, there is still a significant bias on \(\alpha\) parameters.
• We can also note how the new L-cut gets rid of some extreme parameter values compared to the full distribution.

4: DC4 results table

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\).

Turning on dust decorrelation in the model doesn't cause any bias in \(r\) and the recovered \(\Delta_d\) values are centered around 1 (i.e. analysis recovers zero decorrelation). Adding this parameter does increase \(\sigma(r)\) somewhat.

When comparing these results to the old posting, we get an improvement on \(\sigma(r)\) due to the outlier rejection. We report these decrease in % in Table 00bis.

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 = 4.25 \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.

Note here that the value of \(\sigma(r)\) doesn't change much compared to the 04.00 case despite having a higher dust level. This is probably due to the fact that we only have one realization of the foreground sky (see note in the introduction), thus no impact from cosmic variance.

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). Not sure what type of \(\beta_d\) spatial variations are included in this model, but Colin thinks 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.

Note here that the value of \(\sigma(r)\) doesn't change much compared to the 04.00 case despite having a much higher dust level. This is probably due to the fact that we only have one realization of the foreground sky (see note in the introduction), thus no impact from cosmic variance.

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\).

The recovered means seem quite wacky, and \(A_L\) dependent.

Note here that the value of \(\sigma(r)\) doesn't change much compared to the 04.00 case despite having a higher dust level. This is probably due to the fact that we only have one realization of the foreground sky (see note in the introduction), thus no impact from cosmic variance.

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\). 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.

05: Gaussian decorrelated dust model

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 \(\Delta_d\) = 0.85.

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 \(\Delta_d = 1\) only captures part of the statistical penalty.

06: Flauger MHD foregrounds

Out 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.36\). 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. In general, the results for this model look nearly as good as the simple Gaussian foregrounds (sky model 00).

Appendix A: Table Results for CDT report

As can be seen in Figure 4, we seem to have biases even in the case of the Gaussian foreground simulations, mostly for the foreground parameters. In this posting this last has been improved.

Just as for the CDT report, we 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. As we have seen in , the "algorithmic bias" on r is is now in this case basically irrelevant. Note that we also improved the ones on the foreground parameters.