Matrix based E/B separation for DC4 maps.

While full sky linear polarization fields can easily (and uniquely) be decomposed into E and B modes, cut-sky observations with non-uniform coverage and/or mode filtering have aliasing that provoke E to B leakages. The resulting bandpowers then contain pure-E, pure-B, and ambiguous modes that come from a mix of the original E and B modes. We can estimate the bias at the bandpower level by generating many simulations with a known power and computing the mean of these simulations' bandpowers compared to the input. But enhanced sample variance due to the leaked E-modes remains. More advanced estimators have been developed which avoid this increase in variance, by isolating the pure-B (or E) modes from the ambiguous modes.
In the data challenges, our maps are on a partial sky, and our noise is not homogeneous, which is modeled by dividing the noise maps by sqrt(hitcount map), where a hitcount maps depends on the pointing of our telescopes and the size of the field of view. In order to analyse the maps, we currently use the pure-s2hat (or p-s2hat) library, an implementation of the Smith pure estimator for cross-spectra, as introduced in this posting. As we saw in this short posting, while pure-s2hat reduces the E-to-B leakage at intermediate \(\ell\) as expected, there is a spurious low-\(\ell\) power. I think this is expected, coming from the 1/ell^2 dependence of the mixing kernels, and it could be computed and taken care of exactly. This excess of low-\(\ell\) power should be taken care of by our bandpower window functions, which are computed by imputing maps with delta function power of a red spectrum in the pipeline, and computing the bandpower response.

In this posting, we describe another power spectrum method, using matrix-based E/B separation.

Introduction

The BICEP/Keck analysis pipeline uses a method described in the BKVII paper.

The main idea is to build two purification matrices that will project the pure-E and pure-B modes out of an observed map. In BKVII, it is shown that these purification matrices can be constructed by solving the following regularized generalized eigenvalue equation: \begin{align} \left( \mat{\tilde{C}_B} + σ^2\mat I \right) \mat x_i &= λ_i \left( \mat{\tilde{C}_E} + σ^2\mat I \right) \mat x_i , \label{eqn:geneig} \end{align} where \((λ_i,x_i)\) are the eigenvalues and corresponding eigenvectors, \(\mat{\tilde{C}_{(E/B)}}\) are the observed covariance matrices, \(\mat I\) is the identity matrix, and \(σ^2\) is the regularization parameter. The observed covariance matrix is \begin{align} \mat{\tilde{C}} = \mat R \mat C \mat R{}^\top , \label{eqn:obscovmat} \end{align} where the \(\mat R\) is the observation matrix, and \(\mat C\) is the covariance matrix. The observation matrix is here much simpler than in the BK case (where it includes all sorts of effects from data reduction). Here it should simply including the beam smoothing\(^*\), simple filterings\(^*\) and the pixel inverse variance weighting. The pixel-pixel covariance matrix is given by: \begin{align} \mat{C_{E}} \equiv \sum_{\ell m} C_\ell^{EE} \mat Y_{\ell m}^E \mat Y_{\ell m}^{E\dagger} , \label{eqn:covmat} \end{align} and similarly for \(\mat{C_{B}}\). The \(C_\ell\) is somewhat arbitrary here, but would ideally have a dependence similar to the signal. We here use a red spectrum: \(C_\ell = 1/\ell^2\). The covariance matrix is computed using the method from Tegmark & Oliviera-Costa (2001) , Appendix A. We focus only on the polarization.
\({}^*\): The beam smoothing and simple filterings are actually modeled at \(\ell\) level here, since it is a much simpler diagonal multiplication. We instead include it in the \(C_\ell\) used to compute the covariance matrix.

Once we have built these matrices, we can solve for equation \(\ref{eqn:geneig}\). The smallest eigenvalues correspond to the pure-E eigenmodes \(\mat e_i\) which are orthogonal to the B modes, whereas the highest eigenvalues correspond to the pure-B eigenmodes \(\mat b_i\) which are orthogonal to the E modes. The eigenvalues close to 1 are ambiguous modes.
We can then build the purification matrices, the description of which is in this posting (for those who have access to the BK pages). I will try to summarize it in an updated version of the current posting.
These are then applied to input observed maps \(\mat{\tilde{m}} \): \begin{align} \mat{\tilde{m}_{\rm Pure-E}}&= \mat{\Pi_{E} \tilde{m}}, \\ \mat{\tilde{m}_{\rm Pure-B}}&= \mat{\Pi_{B} \tilde{m}}, \\ \label{eqn:maps} \end{align}

The covariance matrix has a dimension \(N_{\rm pix} \times N_{\rm pix}\), so to make the computation tractable we here work with maps at \(N_{\rm side}=128\).

In section 1, we show the observation matrix for the different surveys (DC04 a,b,c and d)
In section 2, we show the observed covariance matrix, focusing on DC04c.
In section 3, we have a look at the eigensystem from equation \(\ref{eqn:geneig}\), and some eigenmodes.
In section 4, we build and plot the purification matrices (pure-E and pure-B projectors).
In section 5, we show pure-E and pure-B maps, as well as their power spectrum.
In section 6, we apply a low-pass filtering to the input maps to reduce high-\(\ell\) aliasing.

1. Observation matrix.

The observing matrix in the case of the S4 pipeline includes the effect from non-uniform coverage, i.e. the hitcount weighting in pixel space. Since we are also downsampling the maps to \(N_{\rm side}=128\), we have to include this in the observation matrix.

Figure 1 shows the different hitcount masks.
As explained here, 04 is the reference mask, with peak value at unity. The other masks have been rescaled to have the same total number of hits (i.e. same sum).
Note that these maps, at \(N_{\rm side}=128\), have respectively 10799, 68616, 9728 and 4915 observed pixels.

Figure 2 shows the different components of the observing matrix:

• The hitcount weighting, at the original resolution with \(N_{\rm side}=512\), which is simply a diagonal matrix.
• The downsampling matrix, which is simply an averaging of 16 pixels into the corresponding pixel in healpix Nested pixelization.

2. Covariance matrix.

We compute the covariance matrix using a EE and BB red spectra as described in \(\ref{eqn:covmat}\).
The covariance matrix is a \(\rm 2 N_{pix} \times 2 N_{pix}\) matrix. Column number n, represents the covariance of pixel number n with the rest of the pixels (both for Q and U).
In figure 3, we show the healpix projection of the column in the middle of the matrix.
Computing the full covariance matrix takes a while, though it's a fully parallelizable task. Some statistics of the runs on the oddyssey cluster are shown here.

Colorbars added on February 17th.. We plot the covariance matrix as is, with a free colorbar, as well as the covariance matrix divided by it maximum, with a symmetric colorbar.

Note that as a first test concerning the fact that we have high ell residuals in the current version with downsampling, we tried computing a covariance matrix where the input red spectrum stops at \(\ell_{\rm max} = 3*128-1=383\) instead of the \(\ell_{\rm max} = 700\) in the original run. In some of the following plots, we present results for both cases.

3. Solving the eigensystem.

Once we have the Q/U covariance matrices, both for E only and B only signals, we can solve the eigensystem shown in equation \(\ref{eqn:geneig}\), using \(\sigma^2 = 0.001 \).
Figure 4 shows the eigensystem, i.e. the sorted eigenvalue distribution, as well as the corresponding eigenvectors.
Figure 5 shows the first, fifth and hundredth eigen modes on the sky, as well as the modes corresponding to the 3 cutoff considered later on.

4. Pure-E and pure-B projectors.

We then build the purification matrices, or pure-E/pure-B projectors. We show such projectors in figure 6.

5. Maps and power spectra.

To test our purification matrices, we use the first unlensed \(\Lambda CDM\) map from the 04 data challenge. These maps don't contain any noise, and by definition only E modes. We expect to have a no power in the output pure-B map, and the input pure-E map to be similar to the input map. Note that we apodize these maps, as we would do in the case of the real analysis, using the inverse variance hitcount map: \begin{align} \mat{\tilde{m}_{\rm Pure-E/B}}&= \mat{\Pi_{E/B} (\tilde{m}*H)/H}., \\ \label{eqn:maps_apod} \end{align} where \(\mat H\) is the apodization map used. We have such maps plotted in figure 7. As we can see, when the beams are tight, we have some high ell residuals in the map.

We then take the power spectra of these Q/U maps, using s2hat, (without purification, so healpix anafast should be giving similar results).
Figure 8 show the EE and BB power spectra for the different maps (input \(N_{\rm side}=512\) maps), pure-E and pure-B maps, as well as the pure-s2hat analysis of the input map and the downsampled input map (at \(N_{\rm side}=128\)).
It seems like the matrix purification is not efficient at high ell. But as can be seen for instance in this figure, s2hat behaves even worse when applied on a downsampled map.
Note also in this figure that the large beam map has a much better purification.
These two effects point towards high-\(\ell\) aliasing.
The two positive things to note are that when doing a fair comparison on the same downsampled input map, the matrix analysis seems to yield better results in terms of purification: first the residual B modes are lower, and the low-\(\ell\) residual in EE and BB is absent.

6. Low-pass filtering of the maps in preprocessing.

As we saw in the previous figures, there seems to be high-\(\ell\) aliasing in the power spectra.
As a first step to avoid this, we low-pass filter the maps. Here I use the original \(a_{\ell m}\) used for the map analysed before and multiply by a beam that has a roll-off at high-\(\ell\). I use a tapered cosine roll-off, for which I can control the \(\ell_{\rm min}\) and \(\ell_{\rm max}\) where the roll-off starts and ends.
We use a roll-off between 373 and 393, and one between 256 and 276.
As we see in Figure 9, the roll-off around \(3*N_{\rm side}-1\) is not strong enough, but the one at \(2*N_{\rm side}\) seems to remove the aliasing.

7. Pipelined analysis of 300 lensed-\(\Lambda\)CDM and noise simulations.

Added on February 17th.
In this section, we compare the mean and std of 300 power spectra using 3 methods: pure-s2hat, s2hat, and the purification matrix method. We analyze the maps after downsampling (and without the low-pass filtering from previous section).

Comments about the figure:

• Note that in all methods, there is still an aliasing from the high ell due to downsampling.
• We see a clear E-to-B leakage in the BB signal in the red spectrum, as already mentioned in this posting.
• While pure-s2hat had a reduced E-to-B leakage, it has a low ell contribution, as mentioned in the introduction.
• The purification matrix method reduces both E-to-B leakage and low ell contribution, which is encouraging.
• Note that the matrix purification yellow curve crosses the pure-s2hat blue curve at \(\ell\simeq 80\), which seems suspicious. It might be due to aliasing effect, so we will have to keep an eye on this in the next iteration using maps at 256.

8. Conclusion and perspective.

In this posting, we show the first implementation of the matrix-based purification method on the S4 maps. It is mostly reusing tools developed for the BK analysis.

• We demonstrate the feasibility of this method, though for now still having some issues from high-\(\ell\) aliasing.
• When compared on the pure-s2hat on the same downsampled maps, the matrix method seems promising.
• It also seems like the low-\(\ell\) excess that we see in pure-s2hat is not present in the matrix analysis.

A few natural follow-ups will be studied:

• While testing the effect of the roll-off in the last figure, we didn't implement such a filter in the observation matrix (or more naturally in the covariance matrix input spectrum, as for the beam).
  In the next run, we will have to carefully chose a roll-off and include it from the first of the purification building.
• Since the roll-off seems tp need to start around \(2*N_{\rm side}\), we might have to push the resolution to \(N_{\rm side}=256\). This might become an issue if we want to analyse the Chile mask with this method.
• Following the last point, it might be interesting to study more carefully what is the highest roll-off we can apply. We would also have to check what is the effect of removing some of the higher bandpowers on the estimation of \(r\).
• For now we only looked at a single map's power spectrum, whereas we are of course mostly interested in the uncertainty on this estimation in terms of \(r\) uncertainty. For this we will need to pipeline the power spectrum taking.