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

HAL
open science

# Pure E and B polarization maps via Wiener filtering

Emory F. Bunn, Benjamin Wandelt

► To cite this version:

Emory F. Bunn, Benjamin Wandelt. Pure E and B polarization maps via Wiener filtering. Physical Review D, 2017, 96 (4), pp.043523. 10.1103/PhysRevD.96.043523. hal-01584648

HAL Id: hal-01584648

https://hal.science/hal-01584648v1

Submitted on 4 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Pure *E* and *B* polarization maps via Wiener filtering

Emory F. Bunn*

*Physics Department, University of Richmond, Richmond, Virginia 23173, USA*

Benjamin Wandelt

Sorbonne Université, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris,
98 bis bd Arago, 75014 Paris, France;

Sorbonne Université, Institut Lagrange de Paris (ILP), 98 bis Boulevard Arago, 75014 Paris, France;
and Departments of Physics and Astronomy, University of Illinois at Urbana-Champaign,
Urbana, Illinois 61801, USA

(Received 12 October 2016; published 30 August 2017)

In order to draw scientific conclusions from observations of cosmic microwave background (CMB) polarization, it is necessary to separate the contributions of the *E* and *B* components of the data. For data with incomplete sky coverage, there are ambiguous modes, which can be sourced by either *E* or *B* signals. Techniques exist for producing “pure” *E* and *B* maps, which are guaranteed to be free of cross-contamination, although the standard method, which involves constructing an eigenbasis, has a high computational cost. We show that such pure maps can be thought of as resulting from the application of a Wiener filter to the data. This perspective leads to far more efficient methods of producing pure maps. Moreover, by expressing the idea of purification in the general framework of Wiener filtering (i.e., maximization of a posterior probability), it leads to a variety of generalizations of the notion of pure *E* and *B* maps, e.g., accounting for noise or other contaminants in the data as well as correlations with temperature anisotropy.

DOI: 10.1103/PhysRevD.96.043523

## I. INTRODUCTION

Characterization of the polarization of the cosmic micro-
wave background (CMB) radiation is an extremely high
priority in cosmology. Polarization has been measured by
the Planck [1,2] and WMAP [3] satellites, as well as by a
number of suborbital telescopes (e.g., [4–7]).

Polarization maps can constrain cosmological models in
a variety of ways, perhaps the most exciting of which is the
potential detection of a stochastic background of gravita-
tional waves produced during an inflationary epoch (e.g.,
[8]). Detection of this signal would constitute a direct
confirmation of inflation and a measurement of the infla-
tionary energy scale, with revolutionary effects on the field.

Detection of this signal, as well as other phenomena such
as gravitational lensing, in CMB polarization maps depends
on the fact that a polarization map can be represented as
the sum of a scalar *E* component and a pseudoscalar *B*
component (the spin-2 analogues of curl-free and diver-
gence-free vector fields respectively) [9–12]. To linear
order in perturbation theory, scalar density perturbations
produce only the *E* component, leaving the *B* channel clear
as a probe of other phenomena such as inflationary
gravitational waves. The *E*-*B* decomposition is also impor-
tant in the analysis of surveys of weak gravitational lensing
(e.g., [13] and references therein).

Because the *B* component is predicted to be at least an
order of magnitude smaller than the *E* component over all
angular scales, detection of *B* modes requires clean
separation of the two components. For data with complete
sky coverage, leakage-free separation of the two compo-
nents is straightforward. However, in a data set with partial
sky coverage, there are “ambiguous modes” that cannot be
uniquely assigned to either *E* or *B* [14–18].

Although the *E*-*B* decomposition is not unique for data
with incomplete sky coverage, one can uniquely decom-
pose any polarization map into three components, dubbed
“pure *E*,” “pure *B*,” and “ambiguous” [15,17]. The pure *E*
(*B*) component lies in the orthogonal complement of the
vector space of all *B* (*E*) modes. The ambiguous compo-
nent is orthogonal to both the pure *E* and pure *B* spaces.
This decomposition has the advantage that any signal found
in the pure *B* map is guaranteed to have come from actual
*B* modes.

Original work on the pure-ambiguous decomposition
involved the construction of orthonormal bases for the
various subspaces, which is a slow and cumbersome
procedure. An alternative view emerges when considering
the CMB polarization analysis problem in the context of
Gibbs sampling [19–21]. Since Gibbs sampling infers the
posterior statistics of the all-sky signal given the data,
separating *E* and *B* modes is trivial for every sample.
The resulting set of *T*, *E*, and *B* map samples represents
the information the data contain. Gibbs sampling does not
require mode decompositions but relies purely on optimal

<footnote>*ebunn@richmond.edu</footnote>
filtering of the data augmented to cover the entire sky. It is therefore natural to ask whether the E/B mode separation problem can be approached from the filtering point of view.

In this paper, we show that the decomposition can be thought of as an application of the familiar Wiener filter, which allows for a much faster implementation. In particular, because the filter can be expressed in terms of operations that are diagonal in either pixel space or spherical harmonic space, efficient techniques such as conjugate gradient solution and the messenger method [22,23] can be applied.

In addition to providing an efficient method of computation, the Wiener filter approach places the pure-ambiguous decomposition into a context that allows it to be generalized in a variety of ways. This approach can be used to generate the original pure and ambiguous maps, but it can also be used to generate more useful generalizations of them that treat noise and other contaminants and the mask together in a unified way, with the mask regarded simply as a region of infinite noise. The resulting filtered pure maps thus suppress noise-dominated modes while simultaneously accounting for the mask.

The method also has a natural generalization to include correlations between temperature and polarization. In particular, it provides a natural way of “purifying” the E map from the portion that is correlated with temperature, giving a clear view of the new information contained in such a map.

There are many other methods of dealing with the problem of EB leakage. Some work at the power spectrum level, without producing real-space maps of the E and B modes [24–28]. Others produce maps in real space of estimates of the derivatives of the purified polarization maps [29–32]. While these are potentially quite valuable, they differ from the approach taken herein, which aims to produce filtered maps of the actual polarization measurements, rather than their derivatives. Ref. [33] assesses the merits of these methods for power-spectrum estimation. Wavelet-based methods [34–36] have also been developed. These approaches require a certain amount of tuning (e.g., careful choice of scale-dependent masks), whereas the Wiener approach allows all filtering to be performed in a principled way from the data and a theoretical model.

The remainder of this paper is organized as follows. Section II establishes notation and provides a brief review of some aspects of the E-B decomposition. Section III shows how the pure-ambiguous decomposition can be implemented as a quadratic minimization problem involving extension of the data into the masked region, allowing us to take advantage of efficient all-sky E-B separation. In Sec. IV, we show that this approach generalizes in a natural way to the application of a Wiener filter. Section V generalizes the previous results to include temperature-correlation. Section VI shows some examples of the method, and Sec. VII contains some brief concluding remarks.

## II. BASIC SETUP AND NOTATION

Our data set consists of measurements of the linear polarization Stokes parameters Q, U at a set of $N_{\text{obs}}$ pixels covering part of the sky. We will consider the correlation with temperature (Stokes I) measurements in Sec. V. As usual, each measurement contains both signal and noise:

$$d_j = s_j + n_j. \quad (1)$$

We assume Gaussian noise with covariance matrix N. In addition to instrumental noise, $n_j$ can include the effects of residual foregrounds or other contaminants in the data.

The index j labels both the pixel location and the Stokes parameter, so the vector $\vec{d}$ has dimension $2N_{\text{obs}}$. The signal can be expressed as a spherical harmonic expansion

$$s_j = s_j^E + s_j^B = \sum_{l,m} (a_{E,lm} Y_{E,lm}^Z(\hat{r}_j) + a_{B,lm} Y_{B,lm}^Z(\hat{r}_j)). \quad (2)$$

Here $\hat{r}_j$ is the location of the pixel corresponding to measurement j and $Z \in \{Q, U\}$ is the Stokes parameter of that measurement. The functions Y are related to the spin-2 spherical harmonics:

$$Y_{E,lm}^Q = Y_{B,lm}^U = -\frac{1}{2}(2Y_{lm} + -2Y_{lm}), \quad (3)$$

$$Y_{B,lm}^Q = -Y_{E,lm}^U = -\frac{1}{2}(2Y_{lm} - -2Y_{lm}). \quad (4)$$

If the data cover a small enough region that the flat-sky approximation is appropriate, then the spherical harmonics can be replaced with plane waves and fast Fourier transforms may be used. See [17] and references therein for further details.

The signal vector can thus be written

$$\vec{s} = \mathbf{Y}_E \vec{e} + \mathbf{Y}_B \vec{b}, \quad (5)$$

where the vectors $\vec{e}$ and $\vec{b}$ contain the coefficients $a_{E,lm}$ and $a_{B,lm}$ respectively, and the matrices $\mathbf{Y}_Z$ contain the spherical harmonics evaluated at the pixel locations. [To be specific, for each element $(\mathbf{Y}_Z)_{ja}$, j labels the pixel location and Stokes parameter, and $\alpha$ labels the index pair $(l, m)$.]

If we make the usual assumption that the data are derived from a statistically isotropic, parity-respecting Gaussian random process, then the theory is completely described by the signal covariance matrix,

$$\mathbf{S} \equiv \langle \vec{s} \vec{s}^\dagger \rangle = \langle \vec{s}^E \vec{s}^{E\dagger} \rangle + \langle \vec{s}^B \vec{s}^{B\dagger} \rangle \equiv \mathbf{S}_E + \mathbf{S}_B, \quad (6)$$

where
$$S_Z = Y_Z C_Z Y_Z^\dagger. \qquad (7)$$

Here $Z \in \{E, B\}$, and the diagonal matrices $C_{ZZ}$ contain the power spectra. Specifically, $(C_{ZZ})_{\alpha\beta} = \delta_{\alpha\beta}C_l^{ZZ}$, where $l$ is the multipole index corresponding to $\alpha$ and $C_l^{ZZ} \equiv \langle a_{Z,lm}^* a_{Z,lm} \rangle$ is the usual power spectrum.

If the data vector $\vec{d}$ covers the whole sky, then the matrices $Y_E, Y_B$ span orthogonal subspaces and, hence, the $E$ and $B$ signals can be estimated independently of one another. In maps with incomplete sky coverage, however, this is not the case: there are “ambiguous modes” that lie in both subspaces simultaneously. It is impossible to say whether power in such a mode came from $E$ or $B$ modes. However, we can divide the signal vector space into three orthogonal subspaces, dubbed pure $E$, pure $B$, and ambiguous spaces. The pure $E$ space is the orthogonal complement of the space spanned by $Y_B$: that is, it consists of modes that are orthogonal (over the observed region) to all possible $B$ maps. The pure $E$ space is defined similarly, and the ambiguous space is orthogonal to both of these. If the data set $\vec{d}$ is projected onto the pure $B$ subspace, then all signal from $E$ modes will be mapped to zero; that is, any power seen in this pure $B$ map (beyond the noise) is guaranteed to have come from $B$ modes.

### III. FAST PURIFICATION

Early work focused on finding bases for the pure and ambiguous spaces by solving an eigenvalue problem, and this method was implemented in the analysis of the BICEP2 data [37]. However, this procedure is slow and cumbersome. In this section, we show an efficient way of finding the pure maps corresponding to a given data set.

Let us ignore noise for the moment and consider just the effects of incomplete sky coverage. Our data contains observations at $N_{obs}$ points on the sky. Let us embed the $2N_{obs}$-dimensional data vector into a larger $2N_{pix}$-dimensional vector space whose pixels cover the entire sky. Let $M$ be the “mask operator” which orthogonally projects onto the space of observed pixels:

$$M_{jk} = \begin{cases} \delta_{jk} & \text{if } j \text{ corresponds to an observed pixel} \\ 0 & \text{otherwise.} \end{cases} \qquad (8)$$

We set the data vector to zero for all unobserved pixels, so $M\vec{d} = \vec{d}$.

One way to produce the pure $B$ map corresponding to $\vec{d}$ is to find the extension of $\vec{d}$ into the unobserved region that minimizes the total $B$ power. We will now describe this procedure in detail and prove that it is equivalent to the usual definition of the pure $B$ map. Naturally, an equivalent statement applies to the pure $E$ map.

Let $\vec{\delta}$ be an extension of $\vec{d}$ into the observed region, so that

$$M\vec{\delta} = M\vec{d}. \qquad (9)$$

Because $\vec{\delta}$ is an all-sky map, we can unambiguously decompose it into $E$ and $B$ components by applying projection operators $P_E$ and $P_B$. We choose $\vec{\delta}$ to minimize the quantity

$$\phi = (P_B \vec{\delta})^2, \qquad (10)$$

subject to the constraint (9). The result is

$$P_B \vec{\delta} = M\vec{\lambda}, \qquad (11)$$

where $\vec{\lambda}$ is a vector of Lagrange multipliers. In other words, $P_B \vec{\delta}$ lies entirely in the observed (unmasked) part of the sky. Because $M^2 = M$,

$$M P_B \vec{\delta} = P_B \vec{\delta}. \qquad (12)$$

We know that

$$\vec{\delta} = P_E \vec{\delta} + P_B \vec{\delta}. \qquad (13)$$

We wish to show that the first term in this expression contains the $E$ modes (pure and ambiguous), and the second contains the pure $B$ modes. The two terms lie in the $E$ and $B$ subspaces, respectively, so all we have to show is that the second term is “pure,” i.e., that it is orthogonal to all possible $E$ modes over the observed region. Let $\vec{e}$ be an arbitrary $E$ mode. The dot product of $\vec{e}$ with $P_B \vec{\delta}$, taken only over the observed pixels, is

$$\vec{e}^\dagger M P_B \vec{\delta} = \vec{e}^\dagger P_B \vec{\delta} = (P_B \vec{e})^\dagger \vec{\delta} = 0, \qquad (14)$$

using Eq. (12), then the Hermiticity of $P_B$, then the fact that $\vec{e}$ is an $E$ mode. So $P_B \vec{\delta}$ is indeed a pure $B$ map.

For full-sky data, implemented in HEALPix for example [38], we can transform easily and quickly back and forth between the pixel basis and the spherical harmonic basis. The projections $P_E, P_B$ are trivial in the latter basis: they are diagonal matrices with ones and zeros along the diagonal. Because of this, the problem above can be solved efficiently using conjugate gradient minimization.

We can generalize this procedure by acknowledging the existence of noise in the maps. Instead of requiring strict agreement between $\vec{\delta}$ and $\vec{d}$ in the observed pixels, we can penalize disagreement in the usual way by introducing an inverse noise matrix $N^{-1}$ and finding $\vec{\delta}_N$ that minimizes

$$\chi^2 = (P_B \vec{\delta}_N)^2 + (\vec{d} - \vec{\delta}_N) N^{-1} (\vec{d} - \vec{\delta}_N). \qquad (15)$$

The inverse noise matrix is diagonal, with $(N^{-1})_{jj} = 0$ for masked pixels. For observed pixels, we set $N^{-1}_{jj} = \sigma_j^{-2}$ for some noise levels $\sigma_j$. In the no-noise limit $\sigma_j \to 0$,
discrepancies in the unmasked region are infinitely penalized, and the result corresponds to the pure map described above:

$$ \vec{\delta} = \lim_{\sigma_j \to 0} \vec{\delta}_N. \qquad (16) $$

One way to prove this formally is to write out $\partial\chi^2/\partial\delta_{NJ} = 0$ for both masked and unmasked pixels. For unmasked pixels, the limit $\sigma_j \to 0$ enforces $\delta_{NJ} \to \delta_j$. For masked pixels, it require $(\mathbf{P}_B\vec{\delta})_j = 0$ as required by Eq. (11).

As we will see in the next section, this correspondence leads to a natural way of thinking of the pure E and B modes as applications of a Wiener filter.

## IV. WIENER FILTER

We begin by recalling the general features of the Wiener filter. Assume that we have a data vector $\vec{d} = \vec{s} + \vec{n}$, where the signal $\vec{s}$ and noise $\vec{n}$ are Gaussian random vectors with zero mean and covariance matrices S and N. The Wiener-filtered (WF) map is the signal $\vec{s}$ that maximizes the posterior probability. The probability density of $\vec{s}$ is proportional to $e^{-\chi^2/2}$, with

$$ \chi^2 = \vec{s}^\dagger \mathbf{S}^{-1} \vec{s} + (\vec{d} - \vec{s})^\dagger \mathbf{N}^{-1} (\vec{d} - \vec{s}). \qquad (17) $$

The WF map is the vector $\vec{s}_W$ that minimizes this:

$$ \vec{s}_W = (\mathbf{S}^{-1} + \mathbf{N}^{-1})^{-1} \mathbf{N}^{-1} \vec{d} = \mathbf{S}(\mathbf{S} + \mathbf{N})^{-1} \vec{d}. \qquad (18) $$

If S contains the covariances of both E and B signal, then the WF map contains the maximum-likelihood estimates for the combined E and B signal. One way to isolate one of these components is to treat the other as if it were a source of noise. For instance, we can get the WF B map by replacing S by S_B and N by S_E + N:

$$ \vec{s}_W^B = (\mathbf{S}_B^{-1} + (\mathbf{S}_E + \mathbf{N})^{-1})^{-1} (\mathbf{S}_E + \mathbf{N})^{-1} \vec{d} \qquad (19) $$

$$ = \mathbf{S}_B(\mathbf{S}_B + (\mathbf{S}_E + \mathbf{N}))^{-1}\vec{d}. \qquad (20) $$

A similar expression applies to the E signal. Note that $\vec{s}_W = \vec{s}_W^E + \vec{s}_W^B$.

These maps include the ambiguous modes. In particular, in the limit of low noise, $\vec{s}_W = \vec{s}_W^E + \vec{s}_W^B$ approaches $\vec{d}$. In practice, since our theory typically assumes that there is much more E power than B power, ambiguous modes with high signal-to-noise are assigned mostly but not entirely to the E map.

The idea behind “pure” modes is to be absolutely sure that there is no cross-contamination. That is, a pure B map is one whose power cannot possibly have come from E modes. One way to produce a WF pure B map is to let the E-mode power tend to infinity. Define

$$ \mathbf{S}(\alpha) = \mathbf{S}_B + \alpha \mathbf{S}_E \qquad (21) $$

to be the signal matrix with E power inflated by a factor $\alpha$. Replace $\mathbf{S}_B + \mathbf{S}_E$ with $\mathbf{S}(\alpha)$ in Eq. (20) to get

$$ \vec{s}_W^B(\alpha) = \mathbf{S}_B(\mathbf{S}(\alpha) + \mathbf{N})^{-1}\vec{d} \qquad (22) $$

$$ = \mathbf{S}_B \mathbf{S}(\alpha)^{-1} (\mathbf{S}(\alpha)^{-1} + \mathbf{N}^{-1})^{-1} \mathbf{N}^{-1} \vec{d}, \qquad (23) $$

In the limit $\alpha \to \infty$, the only modes that survive are those that lie in the null space of $\mathbf{S}_E$. These are the “pure B” modes.

The easy way to evaluate expressions like these is to extend the vector $\vec{d}$ to cover the entire sky, so that we can quickly convert back and forth between the pixel and spherical harmonic bases. We assign infinite noise to the unobserved pixels, so $(\mathbf{N}^{-1})_{jj} = 0$ for such pixels.

In the spherical harmonic basis, the signal matrices are diagonal:

$$ \mathbf{S}_E = \mathrm{diag}(C_2^E, \dots, C_{l_{\max}}^E, 0, \dots, 0) \qquad (24) $$

$$ \mathbf{S}_B = \mathrm{diag}(0, \dots, C_2^B, \dots, C_{l_{\max}}^B), \qquad (25) $$

where we have ordered the spherical harmonic basis to have all E modes first, and there are $2l + 1$ copies of each $C_l^Z$. Therefore,

$$ \mathbf{S}(\alpha)^{-1} = \mathrm{diag}((\alpha C_2^E)^{-1}, \dots, (\alpha C_{l_{\max}}^E)^{-1}, \\ (C_2^B)^{-1}, \dots, (C_{l_{\max}}^B)^{-1}). \qquad (26) $$

In the limit $\alpha \to \infty$, the E terms go to zero. The resulting matrix is the pseudoinverse $\mathbf{S}_B^+$, which is the inverse of $\mathbf{S}_B$ within the subspace of B modes and zero in the orthogonal subspace of E modes.

The WF pure B map is therefore

$$ \vec{s}_W^{\mathrm{pB}} = \lim_{\alpha \to \infty} \vec{s}_W^B(\alpha) = \mathbf{S}_B \mathbf{S}_B^+ (\mathbf{S}_B^+ + \mathbf{N}^{-1})^{-1} \mathbf{N}^{-1} \vec{d}. \qquad (27) $$

Note that $\mathbf{S}_B\mathbf{S}_B^+ = \mathbf{P}_B$ is the operator that projects onto the B subspace. We can describe this procedure in the following way: to get the WF pure B map, we apply the filter (18), assuming infinite signal in the E sector, and then apply the projection $\mathbf{P}_B$ to the result.

The WF pure B map deserves the name “pure,” in the sense that it is derived entirely from the B modes of the true all-sky signal and is independent of any E modes. To see this, consider the unprojected filtered map $(\mathbf{S}_B^+ + \mathbf{N}^{-1})\mathbf{N}^{-1}\vec{d}$. This map minimizes $\chi^2 = \vec{s}^\dagger \mathbf{S}_B^+ \vec{s} + (\vec{d} - \vec{s})^\dagger \mathbf{N}^{-1} (\vec{d} - \vec{s})$. Now, suppose that $\vec{d}$ is derived from only E modes—that is, it can be extended into the unobserved region in a way that has only E power. Then this extension will have $\chi^2 = 0$ and, hence, will be the map that minimizes $\chi^2$. When we then
apply the projection operator $\mathbf{P}_B$ onto the $B$ subspace, the result will be zero. By linearity, therefore, any contribution to $\vec{d}$ that can be derived from $E$ modes contributes zero to the WF pure $B$ map.

Assuming uncorrelated noise, the matrix $\mathbf{N}^{-1}$ is diagonal in the pixel basis, while $\mathbf{S}_B$ and $\mathbf{S}_B^+$ are diagonal in the spherical harmonic basis, leading to efficient ways of evaluating this expression. In particular, the application of the matrix $(\mathbf{S}_B^+ + \mathbf{N}^{-1})^{-1}$ can be performed by either the conjugate gradient or the messenger method.

The condition number of the matrix $\mathbf{S}_B^+ + \mathbf{N}^{-1}$ is formally infinite: there are $E$ modes that lie entirely within the masked region, which lie in the null space of this matrix. This can lead to slow convergence when solving with conjugate gradient methods. One solution is to make the condition number finite, by making $\alpha$ large but finite instead of taking the limit $\alpha \to \infty$ and/or by treating the masked regions as having large (as opposed to infinite) noise. Alternatively, one can use the messenger method, which gracefully handles the case of formally singular matrices.

Suppose that we take the no-noise limit of the above procedure—that is, $(\mathbf{N}^{-1})_{jj} \to \infty$ for observed pixels and zero for unobserved pixels. This forces perfect agreement over the observed region, so it corresponds to a constrained minimization problem similar to the previous section. Indeed, if we also adopt a white-noise (flat) power spectrum for $\mathbf{S}_B$, then the WF pure $B$ map reduces precisely to the original pure $B$ map.

If we take the no-noise limit but do not adopt a white-noise power spectrum, we get different generalizations of the original pure $B$ procedure, with different multipole-space weightings of the output. These maps are not “pure” by the original definition of ref. [17]: they are not orthogonal to all $E$ modes with respect to the usual inner product. However, they do share the more important property of the original pure maps: they are sourced only by true $B$ modes, with no contribution from $E$ modes.

One useful application of the no-noise limit is to provide an initial guess for the full WF calculation. The no-noise problem is lower-dimensional than the full WF (because observed pixels are fixed), so it can be solved quickly. For modes with high signal-to-noise, the result will be close to the WF map, so the output of this calculation can provide a good starting point for the conjugate gradient or messenger minimization.

## V. CORRELATION WITH TEMPERATURE

Thus far, we have assumed that the only relevant data are measurements of the linear polarization (Stokes $Q$, $U$). The polarization is correlated with the CMB temperature anisotropy (Stokes $I$). We now generalize the earlier results to account for this.

Suppose that the data consists of both temperature and polarization measurements. We will write the data vector as $\vec{d} = (\vec{d}_T \ \vec{d}_P)$, with all of the temperature measurements in $\vec{d}_T$ and all of the polarization measurements in $\vec{d}_P$.

As before, we will assume that these vectors cover the entire sky, assigning infinite noise to unobserved pixels. Then the signal covariance matrix can be written simply in the spherical harmonic basis. If the spherical harmonic coefficients are in the order $T, E, B$, then $\mathbf{S}$ can be written in block form as

$$ \mathbf{S} = \begin{pmatrix} \hat{\mathbf{S}}_T & \hat{\mathbf{S}}_X & 0 \\ \hat{\mathbf{S}}_X^\dagger & \hat{\mathbf{S}}_E & 0 \\ 0 & 0 & \hat{\mathbf{S}}_B \end{pmatrix}. \quad (28) $$

The matrix $\hat{\mathbf{S}}_T = \text{diag}(C_2^{TT}, \dots, C_{\text{max}}^{TT})$, and $\hat{\mathbf{S}}_E, \hat{\mathbf{S}}_B, \hat{\mathbf{S}}_X$ similarly contain the $EE, BB$, and $TE$ power spectra. (The hats indicate that these are sub-blocks of the full covariance matrix. For instance, $\hat{\mathbf{S}}_B$ is the nonzero block of the larger matrix $\mathbf{S}_B$ seen in the previous section.)

We can find the pure $B$ map just as before, by replacing $\hat{\mathbf{S}}_E$ by $\alpha \hat{\mathbf{S}}_E$ and taking the limit $\alpha \to \infty$. If we do this while leaving the other spectra, including the $TE$ covariance, fixed, then it is straightforward to show that

$$ \lim_{\alpha \to \infty} \mathbf{S}(\alpha)^{-1} = \begin{pmatrix} \hat{\mathbf{S}}_T^{-1} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \hat{\mathbf{S}}_B^{-1} \end{pmatrix}. \quad (29) $$

The temperature and polarization sectors decouple, and the pure $B$ map is identical to the original expression (27).<fnref n="1" />

For the pure $E$ map, the situation is somewhat different. If our goal is to produce an $E$ map that is free of $B$ contamination, then we take the limit as the $B$ power tends to infinity, and the resulting matrix is

$$ \lim_{\alpha \to \infty} \mathbf{S}(\alpha)^{-1} = \begin{pmatrix} (\hat{\mathbf{S}}_T \quad \hat{\mathbf{S}}_X)^{-1} & 0 \\ (\hat{\mathbf{S}}_X^\dagger \quad \hat{\mathbf{S}}_E)^{-1} & 0 \\ 0 & 0 \end{pmatrix} \quad (30) $$

This says that the natural way to find the pure $E$ map in the presence of temperature data is to apply the WF to the $T$

<footnote n="1">¹One might wonder whether it is correct to hold $\hat{\mathbf{S}}_X$ fixed while taking this limit. If we imagine that the “extra” $E$ power associated with the inflated power spectrum came from a source uncorrelated with $T$, then this is the correct procedure. If, on the other hand, we take the limit while holding the $TE$ correlation fixed, by multiplying $\hat{\mathbf{S}}_X$ by $\alpha^{1/2}$, then the upper left block of $\mathbf{S}(\alpha)^{-1}$ changes to $(\hat{\mathbf{S}}_T - \hat{\mathbf{S}}_X\hat{\mathbf{S}}_E^{-1}\hat{\mathbf{S}}_X^\dagger)^{-1}$. The $TE$ and $B$ sectors still decouple, and the pure $B$ filter is the same.</footnote>
and E signals jointly. In the notation of the previous section, we find a combined “pure TE map”

$$ \vec{s}_W^{pTE} = S_{TE}S_{TE}^{+}(S_{TE}^{+} + N^{-1})N^{-1}\vec{d}, \quad (31) $$

where $S_{TE}$ is the signal covariance matrix containing both T and E signals and their covariance, and $S_{TE}^{+}$ is the pseudoinverse found by restricting to the subspace of T and E modes (i.e., all B modes lie in the null space of both $S_{TE}$ and $S_{TE}^{+}$).

Of course, “purifying” the TE map from contamination by B modes is typically of limited interest. A more useful pure map E map would be one that had been purified with respect to T—that is, a map that contains only the part of the E polarization that is not predicted by the temperature data. The procedure here would be to take the T power to infinity, holding the TE correlation fixed. This corresponds to replacing $\hat{S}_T$, $\hat{S}_X$ by $\alpha\hat{S}_T$, $\alpha^{1/2}\hat{S}_X$, respectively. To get a pure E map, we should simultaneously purify with respect to B. In this case,

$$ S_{TE}^{+} \equiv \lim_{\alpha \to \infty} S(\alpha)^{-1} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & (\hat{S}_E - \hat{S}_X^\dagger \hat{S}_T^{-1} \hat{S}_X)^{-1} & 0 \\ 0 & 0 & 0 \end{pmatrix}. \quad (32) $$

The filtered map is

$$ \vec{s}_W^{pE} = S_E S_{TE}^{+} (S_{TE}^{+} + N^{-1})^{-1} N^{-1} \vec{d}. \quad (33) $$

## VI. EXAMPLES

In this section, we illustrate some of the procedures described above. For simplicity, we begin by considering

examples that lie in a flat square domain, using Fourier transforms instead of spherical harmonic transforms. We present an example on a spherical sky at the end of this section.

For the flat examples, the “full sky” is a square with periodic boundary conditions, pixelized into a 32 × 32 grid. Because this domain has no boundary, there are no ambiguous modes. The EB decomposition can be performed mode by mode in the Fourier basis: E modes have polarization direction parallel and perpendicular to the wave vector, and B modes have polarization oriented at 45° angles. All minimizations and solutions of linear systems were performed with conjugate gradient methods.

We adopt the E and B power spectra shown in the left panel of Fig. 1. The E spectrum is $P_E(k) \propto k^{-2}e^{-k^2\sigma_b^2}$, and the B spectrum is $P_B(k) \propto k^{-1}e^{-k^2\sigma_b^2}$. The exponential terms correspond to smoothing with a Gaussian beam of width $\sigma_b = 1$ pixel (an FWHM of about 2.4 pixels). (We write $P(k)$ rather than $C_l$ to emphasize that these calculations are on a flat patch rather than the sphere.) The power spectra are normalized so that the rms signal per pixel is 1, and the ratio of rms E to rms B is 10.

Figure 2 shows a realization of Gaussian random polarization maps made with these two spectra. The left and center panels show the E and B maps, respectively. (Note that the B map has been multiplied by 10 for visibility.) The right panel shows the sum of the two, with the addition of noise at a level of 0.3 times the rms signal. The gray region shows the mask—i.e., the area that will be presumed to be unobserved in the subsequent analysis.

Figure 3 shows the results of applying the noise-free purification procedure to these maps. As found in previous work, the ambiguous modes have support primarily close to the mask (except for a small number of large-scale modes), so the pure maps are suppressed near the mask.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Power spectra used to create the data used in the examples described in Sec. VI. The left panel shows the power spectra used in the flat-domain examples, and the right panel shows the spectra for the spherical example. In both cases, the solid blue curve is the E spectrum, and the dashed red curve is the B spectrum.</figcaption></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Example polarization maps: an *E* map (left), a *B* map (center), and a combined map containing the sum of the previous two, along with white noise. The noise amplitude is 0.3 times the signal. The gray region shows the mask to be applied in subsequent filtering. Note that the *B* map has been multiplied by 10 for visibility here and in subsequent figures.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 3. Decomposition of the map in Fig. 2. The input was the combined *E* and *B* map, with no noise. The output pure *E* and *B* maps were computed using the noise-free prescription in Sec. III. The ambiguous map is the difference between the input map and the sum of the two pure maps.</figcaption></figure>

For comparison, we calculated the pure *E* and *B* maps by computing an eigenbasis and projecting as described in [17]. The rms difference in the pure maps computed by the two methods differ by about 2%.

We verified that the results of the WF procedure reduce to the original pure maps in the low-noise limit [Eq. (16)] by applying the WF purification procedure to the same maps, assuming a flat power spectrum and a low noise level $\sigma = 10^{-3}$. The rms differences between the pure *E* and *B* maps resulting from this procedure and the maps shown in Fig. 3 were $6 \times 10^{-7}$ (E) and $7 \times 10^{-8}$ (B) times the rms input signal.

Figure 4 shows the pure *E* and *B* WFs applied to the noisy, masked map of Fig. 2. The pure WF simultaneously suppresses modes with low signal-to-noise (i.e., high-frequency modes) and rejects ambiguous modes.

Figure 5 illustrates the difference between the impure and pure WFs. In both panels, the input map is the *E* map from

Fig. 2, with no *B* modes or noise. The left panel shows the result of applying a WF without purification as in Eq. (18). As expected, *E* modes with high signal-to-noise (large scales) are almost fully recovered, while low signal-to-noise, small-scale modes are suppressed. The filtered *B* map is nonzero, showing *E* to *B* leakage. The right panel shows the result of applying the pure *E* and pure *B* WFs. The *B* power is negligible in this case.

We have also implemented the pure WF procedure on the sphere in the HEALPix pixelization [38]. We created a simulated polarization map based on the power spectra shown in the right panel of Fig. 1, which were calculated by CAMB [39] using the best-fit Planck cosmological parameters [40], with a tensor-to-scalar ratio $r = 0.05$. A Gaussian beam with $\sigma_b = 2.7^\circ$ was applied, and maps of Stokes *Q* and *U* were made with HEALPix resolution $N_{side} = 32$. Gaussian noise was added with signal-to-noise of 5 in each pixel. We assumed that only the southern
<figure><img src="image_7.png" /><figcaption>FIG. 4. Wiener-filtered pure *E* and *B* maps. The input was the combined *E* and *B* map, including noise (Fig. 2).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 5. The difference between the impure and pure Wiener filters. The two panels show the results of applying the “impure” WF (18) and the pure *E* and *B* WFs (27) to an input map containing only the *E* polarization signal shown in Fig. 2, with no *B* modes or noise. The solid curve is the *E* power spectrum (see Fig. 1). The points are the mean-square Fourier coefficients of the filtered maps. In the left panel, the points are |*s*<sup>*E*</sup><sub>*W*</sub>(*k*)|<sup>2</sup> (blue) and |*s*<sup>*B*</sup><sub>*W*</sub>(*k*)|<sup>2</sup> (red). The right panel shows |*s*<sup>*p*<sup>*E*</sup><sub>*W*</sub>(*k*)|<sup>2</sup> and |*s*<sup>*p*<sup>*B*</sup><sub>*W*</sub>(*k*)|<sup>2</sup>. The nonzero signal in the pure *B* map is consistent with the numerical accuracy of the conjugate-gradient minimization procedure.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6. Input data for tests of the pure WF on the sphere. All images show orthographic projections of the southern equatorial hemisphere. The left two panels are the masked *Q* and *U* maps, including *E* signal, *B* signal, and noise. The right two panels show only the *B* component.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 7. Maps of Stokes *Q* and *U* for the pure *B* WF map derived from the data shown in Fig. 6.</figcaption></figure>

hemisphere had been observed, and additionally applied
the Planck COMMANDER polarization mask [41], result-
ing in a data set with 40% sky coverage. The resulting data
set is shown in Fig. 6.

The WF pure *B* maps are shown in Fig. 7. As in the flat
domain, the linear systems are solved via preconditioned
conjugate gradient methods, to take advantage of the fact
that the signal and noise matrices are diagonal in the
harmonic and pixel bases respectively. The preconditioner
was taken to be diagonal in the harmonic domain (as it
would be for a map with all-sky coverage). To speed
convergence, we set the *E* power spectrum to 1000 times its
true value (rather than to infinity), and we assumed noise in
the masked region of 20 times the rms signal (rather than
infinity). We verified that increasing these factors made
negligible difference to the final maps.

We verified that the resulting map was pure by applying
the WF to an input data set consisting of only *B* signal
and noise, and to an input containing only the *E* signal. As
expected, the former led to a pure *B* map virtually identical
to the one shown, and the latter was mapped nearly to zero,
with rms fluctuations approximately 1% of the pure *B* map.

VII. CONCLUSIONS

Because the *E* and *B* components of a polarization map
probe different physics, it is important to be able to cleanly
separate the components. A cleaned *B* map that is guar-
anteed to be free of any contamination from the larger *E*
component is of particular value.

The methods we have presented provide efficient ways
to calculate such maps, based on the principle of the
Wiener filter. The method generalizes the original notion
of pure *E* and *B* maps, which considered only the effects
of masking and pixelization, to account for noise simulta-
neously. The method suppresses noise-dominated modes
(as expected of a Wiener filter) while guaranteeing a strict
lack of contamination from *E* signal into the pure *B* map and
vice versa.

The method can be implemented very efficiently, taking
advantage of the sparsity of the *E*-*B* decomposition in the
spherical harmonic basis.

Because the method is based on maximizing the pos-
terior probability, the filter is determined by the assumed
power spectra, mask, and noise properties, with no arbitrary
decisions such as apodization required. The resulting
maps are filtered versions of the original polarization data,
not scalar functions derived from them via a differential
relation.

In addition to the original motivation of eliminating *E*-*B*
leakage, the method described herein can naturally be
extended to include correlations with temperature data,
producing, e.g., “pure *E*” maps that have had the temper-
ature correlation removed.

The Wiener filter depends upon the choice of an input
power spectrum. One can choose a power spectrum
estimate *a priori*, or, as in methods such as Gibbs
sampling [21], samples of the power spectrum can be
obtained based on the data itself. In this case, the
probability distribution of pure WF maps and power
spectra are obtained simultaneously.

The ideas in this paper therefore generalize directly to
power spectrum inference. Polarization auto- and cross-
power spectrum inference as presented in [21] infers *EE*
and *BB* spectra without purification of the maps: the joint
likelihood function *L*(*C*<sub>*l*</sub><sup>*E*</sup>, *C*<sub>*l*</sub><sup>*B*</sup>) of the two spectra contains
all relevant information, including *E*-*B* leakage. However, a
detection of nonzero power in the WF pure *B* map may be
regarded as more robust than a detection derived from a
joint analysis, as the WF pure *B* map “tries as hard as it
can” to hide power coming from the *E* component.

The pure filtered maps are likely to be even more useful
in contexts other than power spectrum estimation. Tests for
foreground contamination, for example, or searches for
non-Gaussianity and statistical anisotropy depend on real-
space maps. Perhaps most importantly, characterization of
the *B*-type polarization produced via gravitational lensing
depends on details of the *B* map that go beyond the power
spectrum and, hence, on real-space pure *B* maps.

ACKNOWLEDGMENTS

E. F. B. is supported by National Science Foundation (NSF) Award No. 1410133. This work has been done within the Labex Institut Lagrange de Paris (ILP) (reference ANR-10-LABX-63) part of the Idex Sorbonne Universités à Paris pour l'Enseignement et la Recherche (SUPER) and received financial support managed by the Agence Nationale de la Recherche, as part of the program Investissements d'avenir under the reference ANR-11-IDEX-0004-02.
[1] R. Adam, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro et al. (Planck Collaboration), Planck 2015 results. IX. Diffuse component separation: CMB maps, Astron. Astrophys. **594**, A9 (2016).

[2] N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, N. Bartolo et al. (Planck Collaboration), Planck 2015 results. XI. CMB power spectra, likelihoods, and robustness of parameters, Astron. Astrophys. **594**, A11 (2016).

[3] C. L. Bennett, D. Larson, J. L. Weiland, N. Jarosik, G. Hinshaw, N. Odegard, K. M. Smith, R. S. Hill, B. Gold, M. Halpern, E. Komatsu, M. R. Nolta, L. Page, D. N. Spergel, E. Wollack, J. Dunkley, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker, and E. L. Wright, Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Final maps and results, Astrophys. J. Suppl. Ser. **208**, 20 (2013).

[4] E. M. Leitch, J. M. Kovac, N. W. Halverson, J. E. Carlstrom, C. Pryke, and M. W. E. Smith, Degree angular scale interferometer 3 year cosmic microwave background polarization results, Astrophys. J. **624**, 10 (2005).

[5] A. C. S. Readhead et al., Polarization observations with the cosmic background imager, Science **306**, 836 (2004).

[6] The Polarbear Collaboration, A measurement of the cosmic microwave background B-mode polarization power spectrum at sub-degree scales with POLARBEAR, Astrophys. J. **794**, 171 (2014).

[7] QUIET Collaboration et al., Second season QUIET observations: Measurements of the cosmic microwave background polarization power spectrum at 95 GHz, Astrophys. J. **760**, 145 (2012).

[8] M. Chiara Guzzetti, N. Bartolo, M. Liguori, and S. Matarrese, Gravitational waves from inflation, Riv. Nuovo Cimento **39**, 399 (2016).

[9] M. Kamionkowski, A. Kosowsky, and A. Stebbins, A Probe of Primordial Gravity Waves and Vorticity, Phys. Rev. Lett. **78**, 2058 (1997).

[10] U. Seljak and M. Zaldarriaga, Signature of Gravity Waves in the Polarization of the Microwave Background, Phys. Rev. Lett. **78**, 2054 (1997).

[11] M. Zaldarriaga and U. Seljak, All-sky analysis of polarization in the microwave background, Phys. Rev. D **55**, 1830 (1997).

[12] M. Kamionkowski, A. Kosowsky, and A. Stebbins, Statistics of cosmic microwave background polarization, Phys. Rev. D **55**, 7368 (1997).

[13] D. Munshi, P. Valageas, L. van Waerbeke, and A. Heavens, Cosmology with weak lensing surveys, Phys. Rep. **462**, 67 (2008).

[14] E. F. Bunn, Detectability of microwave background polarization, Phys. Rev. D **65**, 043003 (2002); Erratum, Phys. Rev. D **66**, 069902 (2002).

[15] A. Lewis, A. Challinor, and N. Turok, Analysis of CMB polarization on an incomplete sky, Phys. Rev. D **65**, 023505 (2001).

[16] A. Lewis, Harmonic E/B decomposition for CMB polarization maps, Phys. Rev. D **68**, 083509 (2003).

[17] E. F. Bunn, M. Zaldarriaga, M. Tegmark, and A. de Oliveira-Costa, E/B decomposition of finite pixelized CMB maps, Phys. Rev. D **67**, 023501 (2003).

[18] E. F. Bunn, Efficient decomposition of cosmic microwave background polarization maps into pure E, pure B, and ambiguous components, Phys. Rev. D **83**, 083003 (2011).

[19] B. D. Wandelt, D. L. Larson, and A. Lakshminarayanan, Global, exact cosmic microwave background data analysis using Gibbs sampling, Phys. Rev. D **70**, 083511 (2004).

[20] J. Jewell, S. Levin, and C. H. Anderson, Application of Monte Carlo algorithms to the Bayesian analysis of the cosmic microwave background, Astrophys. J. **609**, 1 (2004).

[21] D. L. Larson, H. K. Eriksen, B. D. Wandelt, K. M. Górski, G. Huey, J. B. Jewell, and I. J. O'Dwyer, Estimation of polarized power spectra by Gibbs sampling, Astrophys. J. **656**, 653 (2007).

[22] F. Elsner and B. D. Wandelt, Efficient Wiener filtering without preconditioning, Astron. Astrophys. **549**, A111 (2013).

[23] D. Kodi Ramanah, G. Lavaux, and B. D. Wandelt, Wiener filter reloaded: fast signal reconstruction without preconditioning, Mon. Not. R. Astron. Soc. **468**, 1782 (2017).

[24] K. M. Smith, Pseudo-$C_l$ estimators which do not mix E and B modes, Phys. Rev. D **74**, 083002 (2006).

[25] K. M. Smith, Pure pseudo-$C_l$ estimators for CMB B-modes, New Astron. Rev. **50**, 1025 (2006).

[26] K. M. Smith and M. Zaldarriaga, General solution to the E-B mixing problem, Phys. Rev. D **76**, 043001 (2007).

[27] J. Grain, M. Tristram, and R. Stompor, Polarized CMB power spectrum estimation using the pure pseudo-cross-spectrum approach, Phys. Rev. D **79**, 123515 (2009).

[28] A. Challinor and G. Chon, Error analysis of quadratic power spectrum estimates for cosmic microwave background polarization: sampling covariance, Mon. Not. R. Astron. Soc. **360**, 509 (2005).

[29] J. Kim and P. Naselsky, E/B decomposition of CMB polarization pattern of incomplete sky: A pixel space approach, Astron. Astrophys. **519**, A104 (2010).

[30] J. Kim, How to make a clean separation between CMB E and B modes with proper foreground masking, Astron. Astrophys. **531**, A32 (2011).

[31] W. Zhao and D. Baskaran, Separating E and B types of polarization on an incomplete sky, Phys. Rev. D **82**, 023001 (2010).

[32] J. Bowyer, A. H. Jaffe, and D. I. Novikov, MasQU: Finite differences on masked irregular stokes Q, U grids, arXiv: 1101.0520.

[33] A. Ferté, J. Grain, M. Tristram, and R. Stompor, Efficiency of pseudospectrum methods for estimation of the cosmic microwave background B-mode power spectrum, Phys. Rev. D **88**, 023524 (2013).

[34] L. Cao and L.-Z. Fang, A wavelet-Galerkin algorithm of the E/B decomposition of cosmic microwave background polarization maps, Astrophys. J. **706**, 1545 (2009).

[35] B. Leistedt, J. D. McEwen, M. Büttner, and H. V. Peiris, Wavelet reconstruction of E and B modes for CMB polarization and cosmic shear analyses, Mon. Not. R. Astron. Soc. **466**, 3728 (2017).

[36] K. K. Rogers, H. V. Peiris, B. Leistedt, J. D. McEwen, and A. Pontzen, Spin-SILC: CMB polarization component separation with spin wavelets, Mon. Not. R. Astron. Soc. **463**, 2310 (2016).
[37] BICEP2 Collaboration et al., Detection of B-Mode Polari-
zation at Degree Angular Scales by BICEP2, Phys. Rev.
Lett. **112**, 241101 (2014).

[38] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere, Astrophys. J. **622**, 759 (2005).

[39] A. Lewis and A. Challinor, CAMB: Code for anisotropies in the microwave background, Astrophysics Source Code Library (2011).

[40] P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett et al. (Planck Collaboration), Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. **594**, A13 (2016).

[41] R. Adam, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro et al. (Planck Collaboration), Planck 2015 results. IX. Diffuse component separation: CMB maps, Astron. Astrophys. **594**, A9 (2016).