BICEP2/KECK ARRAY. VII. MATRIX BASED E/B SEPARATION APPLIED TO BICEP2 AND THE KECK ARRAY

The values of the Stokes parameters Q and U for a particular location on the sky are dependent on the choice of coordinate system. By rotating the local coordinate system, Q is rotated into U and vice versa. Under rotation by an angle ϕ, the combinations Q + iU and Q − iU transform as:

$$\begin{eqnarray}(Q+{iU})^{\prime} & = & {e}^{-2i\phi }(Q+{iU})\\ (Q-{iU})^{\prime} & = & {e}^{2i\phi }(Q-{iU}).\end{eqnarray} \tag{ 1 }$$

The T, Q, and U fields can be expressed as sums of spin weighted spherical harmonics. While the temperature anisotropies can be broken down into spin-0 harmonics, the polarization field of Q and U must be expressed in terms of spin-2 spherical harmonics (Goldberg et al. 1967):

$$\begin{eqnarray}T({\boldsymbol{r}}) & = & \displaystyle \sum _{{lm}}{a}_{{lm}}^{T}({}_{0}{Y}_{{lm}}({\boldsymbol{r}}))\\ (Q+{iU})({\boldsymbol{r}}) & = & \displaystyle \sum _{{lm}}{a}_{+2,{lm}}({}_{+2}{Y}_{{lm}}({\boldsymbol{r}}))\\ (Q-{iU})({\boldsymbol{r}}) & = & \displaystyle \sum _{{lm}}{a}_{-2,{lm}}({}_{-2}{Y}_{{lm}}({\boldsymbol{r}})),\end{eqnarray} \tag{ 2 }$$

where ${}_{\pm 2}{Y}_{{lm}}$ are the spin-2 case of spin weighted spherical harmonics, and the spin-0 case are the normal spherical harmonics, ${}_{0}{Y}_{{lm}}$. Since Q + iU and Q − iU are affected by rotations of the coordinate system, it is convenient to express the coefficients of the spin-2 spherical harmonics using a set of coordinate independent scalar a^E_lm coefficients and pseudo-scalar a^B_lm coefficients:

$$\begin{eqnarray}{a}_{{lm}}^{E} & \equiv & -({a}_{+2,{lm}}+{a}_{-2,{lm}})/2\\ {a}_{{lm}}^{B} & \equiv & -i({a}_{+2,{lm}}-{a}_{-2,{lm}})/2.\end{eqnarray} \tag{ 3 }$$

We also define two combinations of spin-2 spherical harmonics:

$$\begin{eqnarray}{X}_{1,{lm}} & \equiv & {(}_{+2}{Y}_{{lm}}+{}_{-2}{Y}_{{lm}})/2\\ {X}_{2,{lm}} & \equiv & {(}_{+2}{Y}_{{lm}}-{}_{-2}{Y}_{{lm}})/2.\end{eqnarray} \tag{ 4 }$$

We can use the coefficients in Equation (3) and the combinations in Equation (4) to construct real space forms of T, Q, and U fields, according to Equation (2):

$$\begin{eqnarray}T({\boldsymbol{r}}) & = & \displaystyle \sum _{{lm}}{a}_{{lm}}^{T}({}_{0}{Y}_{{lm}}({\boldsymbol{r}}))\\ Q({\boldsymbol{r}}) & = & -\displaystyle \sum _{{lm}}({a}_{{lm}}^{E}{X}_{1,{lm}}({\boldsymbol{r}})+{{ia}}_{{lm}}^{B}{X}_{2,{lm}}({\boldsymbol{r}}))\\ U({\boldsymbol{r}}) & = & -\displaystyle \sum _{{lm}}({a}_{{lm}}^{B}{X}_{1,{lm}}({\boldsymbol{r}})-{{ia}}_{{lm}}^{E}{X}_{2,{lm}}({\boldsymbol{r}})).\end{eqnarray} \tag{ 5 }$$

Using these relations, we can write the polarization field as a vector:

$$\begin{eqnarray}P({\boldsymbol{r}}) & \equiv & \left(\begin{array}{c}Q({\boldsymbol{r}})\\ U({\boldsymbol{r}})\end{array}\right)\\ & = & -\displaystyle \sum _{{lm}}\left[\begin{array}{c}{a}_{{lm}}^{E}{X}_{1,{lm}}({\boldsymbol{r}})+{{ia}}_{{lm}}^{B}{X}_{2,{lm}}({\boldsymbol{r}})\\ {a}_{{lm}}^{B}{X}_{1,{lm}}({\boldsymbol{r}})-{{ia}}_{{lm}}^{E}{X}_{2,{lm}}({\boldsymbol{r}})\end{array}\right]\\ & = & -\displaystyle \sum _{{lm}}\left[{a}_{{lm}}^{E}\left(\begin{array}{c}{X}_{1,{lm}}({\boldsymbol{r}})\\ -{{iX}}_{2,{lm}}({\boldsymbol{r}})\end{array}\right)+{a}_{{lm}}^{B}\left(\begin{array}{c}{{iX}}_{2,{lm}}({\boldsymbol{r}})\\ {X}_{1,{lm}}({\boldsymbol{r}})\end{array}\right)\right]\\ & = & -\displaystyle \sum _{{lm}}[{a}_{{lm}}^{E}{Y}_{{lm}}^{E}({\boldsymbol{r}})+{a}_{{lm}}^{B}{Y}_{{lm}}^{B}({\boldsymbol{r}})],\end{eqnarray} \tag{ 6 }$$

where Y^E_lm and Y^B_lm have been introduced and defined in the last step. On the full sphere, Y^E_lm and Y^B_lm are orthogonal:

$$\begin{eqnarray}&&{\int }_{{S}^{2}}{Y}_{{lm}}^{E}({\boldsymbol{r}})\cdot {Y}_{l^{\prime} m^{\prime} }^{B}({\boldsymbol{r}}){dS}=0,\end{eqnarray} \tag{ 7 }$$

for all $l,l^{\prime}$ and $m,m^{\prime}$.

The inner product of two polarization fields is defined as:

$$\begin{eqnarray}&&{\boldsymbol{P}}\cdot {\boldsymbol{P}}^{\prime} \equiv {\int }_{{\rm{\Omega }}}{\boldsymbol{P}}\cdot {\boldsymbol{P}}^{\prime} d{\rm{\Omega }},\end{eqnarray} \tag{ 8 }$$

where Ω is the manifold on which the polarization field is defined: for the full sky it is the celestial sphere. In pixelized maps, the vector space of a polarization field has a finite dimension: twice the number of pixels in the map.

As demonstrated in Equation (7), E and B-mode polarization fields on the full sky are orthogonal. However, experiments produce Q/U maps of portions of the sky, and often filter spatial modes out of these maps. We define the term "observed" maps or modes to refer to these incomplete measurements of the true sky.

The spaces of observed E-modes and B-modes are non-orthogonal. The overlapping subspace between the two is called the ambiguous space. We cannot tell whether signal in the ambiguous subspace came from full sky E-modes or full sky B-modes.

The solution is to decompose vector fields on an observed manifold into three subspaces: "pure" E-modes, "pure" B-modes, and ambiguous modes. Pure E and B-modes are subspaces of the polarization vector space of a particular manifold, defined as:

• 1.  
  A pure B-mode is orthogonal to observed E-modes.
• 2.  
  A pure E-mode is orthogonal to observed B-modes.

Therefore, a pure B-mode is one that has no E to B leakage: neither pure E-modes nor ambiguous modes contribute to it.

We choose a HEALPix pixelization scheme (Gorski et al. 2005) for the input maps, ${\boldsymbol{m}}$, because it has equal area pixels on the sphere and is widely used in the cosmology community.

A true sky signal is represented by the map ${{\boldsymbol{m}}}^{{\boldsymbol{o}}}=[{T}_{{xy}}^{o},{Q}_{{xy}}^{o},{U}_{{xy}}^{o}]$, where ($x,y$) are the (R.A., decl.) coordinates of the map. Using synfast ²¹ , the unobserved input map is convolved with the array averaged beam function, $\bar{{\boldsymbol{ \mathcal B }}}$, constructed from measurements of the beam function of all detectors in the array:

$$\begin{eqnarray}&&\left[\begin{array}{c}{T}_{{xy}}\\ {Q}_{{xy}}\\ {U}_{{xy}}\end{array}\right]=\bar{{\boldsymbol{ \mathcal B }}}\ast {{\boldsymbol{m}}}^{{\boldsymbol{o}}}={\bar{{ \mathcal B }}}_{{xy}}\ast \left[\begin{array}{c}{T}_{{xy}}^{o}\\ {Q}_{{xy}}^{o}\\ {U}_{{xy}}^{o}\end{array}\right].\end{eqnarray} \tag{ 21 }$$

The input map vector is found by reforming the beam convolved two-dimensional map into a one-dimensional vector, ${\boldsymbol{m}}$, of length $3j$, where $j=1...{n}_{p}$, for n_p pixels in the input map:

$$\begin{eqnarray}&&{\boldsymbol{m}}\equiv \left[\begin{array}{c}{T}_{j}\\ {Q}_{j}\\ {U}_{j}\end{array}\right].\end{eqnarray} \tag{ 22 }$$

The observing strategies for Bicep2 and the Keck Array are very similar and borrow heavily from Bicep1. All three experiments target a region of sky centered at a right ascension of 0° and declination of −575. A detailed description of the scan strategy is contained in BICEP2 Collaboration II (2014).

• 1.  
  Halfscans: During normal observations, the telescope scans in azimuth at a constant elevation. The scan speed of 28 s⁻¹ in azimuth places the targeted multipoles of $20\lt l\lt 200$ at temporal frequencies less than 1 Hz. Each scan covers 642 in azimuth, at the end of which the telescope stops and reverses direction in azimuth and scans back across the field center. A scan in a single direction is known as a "halfscan."
• 2.  
  Scansets: Halfscans are grouped into sets of ≈100 halfscans, which are known as "scansets." The scan pattern deliberately covers a fixed range in azimuth within each scanset, rather than a fixed range in right ascension. Over the course of the 50 minute scanset, Earth's rotation results in a relative drift of azimuthal coordinates and right ascension of about 125. At the end of each scanset the elevation is offset by 025, and a new scanset commences. The telescope steps in 025 elevation increments between each scanset. All observations take place at 20 elevation steps, with a boresight pointing ranging in elevation between 55° and 5975. The geographic location of the telescope, near the South Pole, means that elevation and declination are approximately interchangeable.
• 3.  
  Phases: Scansets are grouped together into sets known as "phases." For Bicep2 and the Keck Array, CMB phases consist of ten scansets, comprising 9 hours of observations. CMB phases are grouped into seven types and each type has a unique combination of elevation offset and azimuthal position.
• 4.  
  Schedules: The third degree of freedom in the Bicep2 and Keck Array telescope mounts is a rotation about the boresight, referred to as "deck rotation." The polarization angles relative to the cryostats are fixed, so rotating in deck angle allows detector pairs to observe at multiple polarization angles.A "schedule" typically consists of a set of phases at a particular deck angle. The deck angle is rotated between schedules. There is typically one schedule per fridge cycle, occurring every ∼3 days for Bicep2 and ∼2 days for the Keck Array.

The Bicep2 and Keck Array detectors consist of pairs of nominally co-pointed, orthogonal, polarization sensitive phased array antennas coupled to TES bolometers (BICEP2/Keck & Spider Collaborations 2015). The signal in the timestream from detector "A" is:

$$\begin{eqnarray}&&{\tau }_{t}^{A}={T}_{t}+\mathrm{cos}(2{{\rm{\Psi }}}_{t}^{A}){Q}_{t}+\mathrm{sin}(2{{\rm{\Psi }}}_{t}^{A}){U}_{t},\end{eqnarray} \tag{ 23 }$$

where $[{T}_{t},{Q}_{t},{U}_{t}]$ are the Stokes parameters of the beam convolved sky signal for timestream sample t. A timestream consists of n_t time ordered measurements of the sky, $t=1...{n}_{t}$. ${{\rm{\Psi }}}^{A}$ is the angle the "A" antenna makes with the $Q,U$ axis on the sky. For the HEALPix polarization convention, this axis is a vector pointing toward the north celestial pole.

The relative gain normalized "A" timestream is summed and differenced with the normalized timestream from its orthogonal partner "B":

$$\begin{eqnarray}{s}_{t} & = & \displaystyle \frac{1}{2}({\tau }_{t}^{A}+{\tau }_{t}^{B})={T}_{t}+{\alpha }_{t}^{+}{Q}_{t}+{\beta }_{t}^{+}{U}_{t}\\ {d}_{t} & = & \displaystyle \frac{1}{2}({\tau }_{t}^{A}-{\tau }_{t}^{B})={\alpha }_{t}^{-}{Q}_{t}+{\beta }_{t}^{-}{U}_{t}.\end{eqnarray} \tag{ 24 }$$

The variables α and β are defined by:

$$\begin{eqnarray}{\alpha }_{t}^{\pm } & \equiv & \displaystyle \frac{1}{2}[\mathrm{cos}(2{{\rm{\Psi }}}_{t}^{A})\pm \mathrm{cos}(2{{\rm{\Psi }}}_{t}^{B})]\\ {\beta }_{t}^{\pm } & \equiv & \displaystyle \frac{1}{2}[\mathrm{sin}(2{{\rm{\Psi }}}_{t}^{A})\pm \mathrm{sin}(2{{\rm{\Psi }}}_{t}^{B})],\end{eqnarray} \tag{ 25 }$$

where ${{\rm{\Psi }}}^{B}$ is the angle the "B" antenna makes with the $Q,U$ axis on the sky. Assuming that "A" and "B" are perfectly co-pointed and orthogonal, the signal portion of the timestream vectors can be described with a transformation, ${A}_{{tj}}$, from the input map pixel (with index j) to the timestream sample (with index t):

$$\begin{eqnarray}{s}_{t} & = & {A}_{{tj}}{T}_{j}\\ {d}_{t} & = & {\alpha }_{t}^{-}{A}_{{tj}}{Q}_{j}+{\beta }_{t}^{-}{A}_{{tj}}{U}_{j},\end{eqnarray} \tag{ 26 }$$

where the terms ${\alpha }^{+}$ and ${\beta }^{+}$ in the pair sum timestream cancel due to the orthogonal orientation of the "A" and "B" detectors. The signal-plus-noise timestreams, in vector notation, are:

$$\begin{eqnarray}&&{\boldsymbol{s}}=\displaystyle \frac{1}{2}({{\boldsymbol{n}}}_{{\boldsymbol{A}}}+{{\boldsymbol{n}}}_{{\boldsymbol{B}}})+{\boldsymbol{A}}[{\boldsymbol{T}}]\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}{\boldsymbol{d}}=\displaystyle \frac{1}{2}({{\boldsymbol{n}}}_{{\boldsymbol{A}}}-{{\boldsymbol{n}}}_{{\boldsymbol{B}}})+\left[\begin{array}{cc}{{\boldsymbol{\alpha }}}^{-} & {{\boldsymbol{\beta }}}^{-}\end{array}\right]\left[\begin{array}{cc}{\boldsymbol{A}} & 0\\ 0 & {\boldsymbol{A}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right].\end{eqnarray} \tag{ 28 }$$

where ${{\boldsymbol{n}}}_{{\boldsymbol{A}}}$ and ${{\boldsymbol{n}}}_{{\boldsymbol{B}}}$ are the time ordered noise components of detector "A" and detector "B," assuming the noise is uncorrelated with the pointing of the detector pair. For signal only simulations, ${{\boldsymbol{n}}}_{{\boldsymbol{A}}}$ and ${{\boldsymbol{n}}}_{{\boldsymbol{B}}}$ can be ignored.

The matrix $[{{\boldsymbol{\alpha }}}^{-}\ {{\boldsymbol{\beta }}}^{-}]$ contains the information about the orientation of a pair's antennas relative to Q and U defined on the sky. We call it the detector orientation matrix. The combination:

$$\begin{eqnarray}\left[\begin{array}{cc}{{\boldsymbol{\alpha }}}^{-} & {{\boldsymbol{\beta }}}^{-}\end{array}\right]\left[\begin{array}{cc}{\boldsymbol{A}} & 0\\ 0 & {\boldsymbol{A}}\end{array}\right]\end{eqnarray} \tag{ 29 }$$

transforms input $Q,U$ maps into a pair difference timestream. $[{{\boldsymbol{\alpha }}}^{-}\ {{\boldsymbol{\beta }}}^{-}]$ is constructed from two diagonal matrices, ${{\boldsymbol{\alpha }}}^{-}$ and ${{\boldsymbol{\beta }}}^{-}$, which are filled with the sine and cosine of the detector orientations at each time sample. A graphical representation of the detector orientation matrix is shown in Figure 1. (Additional steps accounting for polarization efficiency and pair non-orthogonality are absorbed into a normalization correction to the pair difference timestream.)

Zoom In Zoom Out Reset image size

The matrix ${\boldsymbol{A}}$ = ${A}_{{tj}}$ represents the timestream forming matrix for a detector pair. It transforms the input temperature map, T_j, into the signal component of the pair sum timestream, s_t. A graphical representation of the timestream forming matrix is shown in Figure 2.

Zoom In Zoom Out Reset image size

To create timestreams with smooth transitions at pixel boundary crossings, the input maps should have a resolution higher than the spatial band limit imposed by the beam function. For this reason, N_side = 512 HEALPix maps are used, whose pixels have a Nyquist frequency $\gtrsim 2\times$ the band limit of the Bicep2 and Keck Array 150 GHz beam function.

The current Bicep2 and Keck Array CMB observations fall within the region of sky bounded in right ascension by $-{3}^{{\rm{h}}}{40}^{{\rm{m}}}\lt \alpha \lt {3}^{{\rm{h}}}{40}^{{\rm{m}}}$ and in declination by $-70^\circ \lt \delta \lt -45^\circ$. This region contains ${n}_{p}={\rm{111,593}}$ pixels in an N_side = 512 HEALPix map. The number of samples in a scanset is typically ${n}_{t}\approx {\rm{43,000}}$.

The simplest form of ${\boldsymbol{A}}$ performs nearest neighbor interpolation of the HEALPix maps, in which case ${\boldsymbol{A}}$ is (${n}_{t}\times {n}_{p}$) and is filled with ones where the detector pair is pointed and zeros otherwise.

A more sophisticated form of ${\boldsymbol{A}}$ performs Taylor interpolation on the HEALPix map, in which case ${\boldsymbol{A}}$ is $\left({n}_{t}\times \tfrac{\lambda (\lambda +1)}{2}{n}_{p}\right)$, where λ is the order of the Taylor polynomial used in interpolation. In this case, ${\boldsymbol{A}}$ is a matrix that performs Taylor interpolation, allowing sub-pixel accuracy to be recovered from the input map, and ${\boldsymbol{m}}$ must also contain derivatives of the true sky temperature and polarization field. This matrix is used to build the deprojection templates in Section 4.10 but is not used for forming timestreams because it increases the dimensions of the observation matrix, making the computation of the observation matrix more difficult.

To remove low frequency atmospheric noise from the data, a third order polynomial is fit and subtracted from each halfscan in the timestreams. Since each halfscan traces an approximately constant elevation trajectory across the target field, the polynomial filter removes power only in the right ascension direction of the maps. In multipole, l, the third order polynomial filter typically rolls off power below $l\lt 40$. This can be represented by a "filtering matrix," ${\boldsymbol{F}}$, which is block diagonal with the block size being the temporal length of a halfscan. Each block is composed of a matrix:

$$\begin{eqnarray}&&{\boldsymbol{F}}={\boldsymbol{I}}-{\boldsymbol{V}}{({{\boldsymbol{V}}}^{\top }{\boldsymbol{V}})}^{-1}{{\boldsymbol{V}}}^{\top },\end{eqnarray} \tag{ 30 }$$

where ${\boldsymbol{I}}$ is the identity matrix and ${\boldsymbol{V}}$ is the same third order Vandermonde matrix for each halfscan of equal length. The Vandermonde matrix is defined as:

$$\begin{eqnarray}&&{V}_{{tj}}={x}_{t}^{j-1},\end{eqnarray} \tag{ 31 }$$

where j = 4 for a third order filter and x_t are the coordinate locations.

For Bicep2 and the Keck Array, x_t is a vector of the relative azimuthal location of each sample in the halfscan. A representation of the polynomial filtering matrix is shown in Figure 3.

Zoom In Zoom Out Reset image size

Scan-synchronous subtraction removes signal in the timestreams that is fixed relative to the ground rather than moving with the sky. These azimuthally fixed signals are decoupled from signals rotating with the sky by the scan strategy, which observes over a fixed range in azimuth as the sky slides by (as described in Section 4.2). A template of the mean azimuthal signal is subtracted from the timestreams for each scan direction. This procedure can be represented as a matrix operator, referred to as a "scan-synchronous signal matrix."

The mean azimuthal signal can be found using a matrix ${\boldsymbol{X}}$ = ${X}_{{tt}^{\prime} }$, for which each row is only filled for entries containing the same azimuthal pointing as the diagonal entry. The scan-synchronous signal matrix subtracts off the mean azimuthal signal:

$$\begin{eqnarray}&&{\boldsymbol{G}}={\boldsymbol{I}}-{\boldsymbol{X}},\end{eqnarray} \tag{ 32 }$$

where ${\boldsymbol{I}}$ is the identity matrix. A graphical representation of the scan-synchronous signal removal matrix is shown in Figure 4. Note that while the ${\boldsymbol{F}}$ matrix is block diagonal and sparse, and the ${\boldsymbol{G}}$ matrix is sparse, once the two are combined, the resulting filter matrix is neither sparse nor block diagonal, making matrix operations more computationally demanding.

Zoom In Zoom Out Reset image size

The timestreams are weighted based on the measured inverse variance of each scanset. Pair sum and pair difference are weighted separately from weights calculated from the two timestreams, ${{\boldsymbol{w}}}^{+}$ and ${{\boldsymbol{w}}}^{-}$. The scheme assigns lower weight to particularly noisy channels and periods of bad weather. This choice of weighting is not a fully "optimal" map maker (Tegmark 1997), but instead represents a practical solution that avoids calculating and inverting a large noise covariance matrix. The weighting is represented by a matrix whose diagonal is filled with the vector ${{\boldsymbol{w}}}^{+}$ = ${w}_{{tt}}^{+}$ for pair sum and ${{\boldsymbol{w}}}^{-}$ = ${w}_{{tt}}^{-}$ for pair difference, shown in Figure 5.

Zoom In Zoom Out Reset image size

Ignoring noise and combining all the operators of Sections 4.5, 4.6 and 4.7, the sum and difference timestreams in Equation (26) are tranformed to the filtered timestreams:

$$\begin{eqnarray}\tilde{{\boldsymbol{s}}} & = & {{\boldsymbol{w}}}^{+}{\boldsymbol{GFA}}[{\boldsymbol{T}}]\\ \tilde{{\boldsymbol{d}}} & = & {{\boldsymbol{w}}}^{-}{\boldsymbol{GF}}\left[\begin{array}{cc}{{\boldsymbol{\alpha }}}^{-} & {{\boldsymbol{\beta }}}^{-}\end{array}\right]\left[\begin{array}{cc}{\boldsymbol{A}} & 0\\ 0 & {\boldsymbol{A}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right],\end{eqnarray} \tag{ 33 }$$

the second of which is graphically represented in Figure 6.

Zoom In Zoom Out Reset image size

The timestream quantities ${\boldsymbol{s}}$ and ${\boldsymbol{d}}$ are converted to maps by the pointing matrix, ${\boldsymbol{\Lambda }}$ = ${{\rm{\Lambda }}}_{{it}}$. If the pixelization of the input maps were identical to the output maps, the pointing matrix would be the transpose of the timestream forming matrix:

$$\begin{eqnarray*}&&{\boldsymbol{\Lambda }}={{\boldsymbol{A}}}^{\top }.\end{eqnarray*}$$

As discussed in Section 4.1, the input maps are HEALPix N_side = 512. However, the Bicep maps instead use a simple rectangular grid of pixels in R.A. and decl.: the size of the pixels is 025 in decl., with the pixel size in R.A. set to be equivalent to 025 of arc at the mid-declination of the map, resulting in 236 × 100 = 23,600 pixels. If the Bicep maps are naively used as "flat maps" then projection distortions are inherent. However, note that the ${\boldsymbol{A}}$ and ${\boldsymbol{\Lambda }}$ matrices together fully encode the mapping from the underlying curved sky to the observed map pixels, allowing such distortions to be accounted for. Figure 2 shows ${\boldsymbol{A}}$ for a single detector over a scanset and Figure 7 shows ${\boldsymbol{\Lambda }}$ for a single detector over a scanset.

Zoom In Zoom Out Reset image size

The pointing matrix for a single detector pair can be used to construct a pair sum "pairmap:"

$$\begin{eqnarray}&&{{\boldsymbol{m}}}_{T}={\boldsymbol{\Lambda }}{\boldsymbol{s}}={\boldsymbol{\Lambda }}{{\boldsymbol{w}}}^{+}{\boldsymbol{GFA}}[{\boldsymbol{T}}].\end{eqnarray} \tag{ 34 }$$

The pair difference timestream is converted into pairmaps using two copies of the pointing matrix. The two pair difference pairmaps correspond to linear combinations of Stokes Q and U:

$$\begin{eqnarray}&&\left[\begin{array}{c}{{\boldsymbol{m}}}_{{{\boldsymbol{\alpha }}}^{-}}\\ {{\boldsymbol{m}}}_{{{\boldsymbol{\beta }}}^{-}}\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{\Lambda }} & 0\\ 0 & {\boldsymbol{\Lambda }}\end{array}\right]\left[\begin{array}{c}{{\boldsymbol{\alpha }}}^{-}\\ {{\boldsymbol{\beta }}}^{-}\end{array}\right]{{\boldsymbol{w}}}^{-}{\boldsymbol{GF}}\left[\begin{array}{cc}{{\boldsymbol{\alpha }}}^{-} & {{\boldsymbol{\beta }}}^{-}\end{array}\right]\left[\begin{array}{cc}{\boldsymbol{A}} & 0\\ 0 & {\boldsymbol{A}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right].\end{eqnarray} \tag{ 35 }$$

For later convenience in abbreviating this equation, we define:

$$\begin{eqnarray}{\boldsymbol{ \mathcal Q }}\equiv \left[\begin{array}{cc}{\boldsymbol{\Lambda }} & 0\\ 0 & {\boldsymbol{\Lambda }}\end{array}\right]\left[\begin{array}{c}{{\boldsymbol{\alpha }}}^{-}\\ {{\boldsymbol{\beta }}}^{-}\end{array}\right]{{\boldsymbol{w}}}^{-}{\boldsymbol{GF}}\left[\begin{array}{cc}{{\boldsymbol{\alpha }}}^{-} & {{\boldsymbol{\beta }}}^{-}\end{array}\right]\left[\begin{array}{cc}{\boldsymbol{A}} & 0\\ 0 & {\boldsymbol{A}}\end{array}\right].\end{eqnarray} \tag{ 36 }$$

A potential systematic concerning polarization measurements is the leakage of unpolarized signal into polarized signal. In the case of CMB polarization, this takes the form of the relatively bright temperature anisotropy leaking into the much fainter polarization anisotropy. The leakage is caused by imperfect differencing between the orthogonal pairs of detectors. The beam functions can be well approximated by elliptical Gaussians, the difference of which correspond to gain, pointing, width and ellipticity (Hu et al. 2003; Shimon et al. 2008).

The Bicep2 and Keck Array pipeline removes leaked temperature signal from the polarization signal using linear regression to fit leakage templates to the polarization data. This method allows the beam mismatch parameters to be fitted directly from the CMB data itself, rather than relying on external calibration data sets, and is robust to temporal variations of the beam mismatch.

The templates used in the regression are constructed from Planck 143 GHz temperature maps.²² These maps contain both CMB and foreground emission at approximately the Bicep2 band. The noise in Planck 143 GHz is significantly subdominant to the CMB temperature anisotropy. For a full description and derivation of the deprojection technique, see Aikin (2013), Sheehy (2013), BICEP1 Collaboration (2014), and BICEP2 Collaboration III (2015). In this section, the entire deprojection algorithm is re-cast as a matrix operation.

For the purposes of generating deprojection templates, we use a timestream forming matrix that performs Taylor expansion around the nearest pixel center to the detector pointing location. The Taylor interpolating matrix produces higher fidelity timestreams than a nearest neighbor matrix. This is important for the deprojection algorithm since small displacements in beam position are responsible for the systematic effect that is removed. Without Taylor interpolation, pixel boundary discontinuities introduce noise and limit the effectiveness of deprojection.

A Taylor polynomial of order λ has $\tfrac{\lambda (\lambda +1)}{2}$ terms, so the dimensions of the input map vector for second order interpolation is $1\times 6{n}_{p}$. Using Equation (21), the input maps are convolved with the array averaged beam function. The smoothing is done using synfast, which contains the ability to output derivatives of the temperature (and polarization) field. Because the beam is applied first, the output derivatives are less noisy than they would be in the raw maps.

The maps are of the form:

$$\begin{eqnarray}&&{\boldsymbol{\Theta }}=\left[\begin{array}{c}{\boldsymbol{T}}\\ {{\boldsymbol{\nabla }}}_{\theta }{\boldsymbol{T}}\\ {{\boldsymbol{\nabla }}}_{\phi }{\boldsymbol{T}}\\ {{\boldsymbol{\nabla }}}_{\theta \theta }{\boldsymbol{T}}\\ {{\boldsymbol{\nabla }}}_{\phi \phi }{\boldsymbol{T}}\\ {{\boldsymbol{\nabla }}}_{\theta \phi }{\boldsymbol{T}}\end{array}\right],\end{eqnarray} \tag{ 37 }$$

where θ and ϕ are the HEALPix map's latitude and longitude.

Using the temperature map and its derivatives, we can find the Taylor interpolated temperature timestream by replacing ${\boldsymbol{A}}$ with a Taylor interpolating matrix, ${{\boldsymbol{A}}}^{\prime }$:

$$\begin{eqnarray}{{\boldsymbol{A}}}^{\prime }=\left[\begin{array}{cccccc}{\boldsymbol{A}} & {\boldsymbol{A}}{\boldsymbol{\Delta }}{\boldsymbol{\theta }} & {\boldsymbol{A}}{\boldsymbol{\Delta }}{\boldsymbol{\phi }} & {\boldsymbol{A}}\tfrac{{\boldsymbol{\Delta }}{{\boldsymbol{\theta }}}^{2}}{2} & {\boldsymbol{A}}\tfrac{{\boldsymbol{\Delta }}{{\boldsymbol{\phi }}}^{2}}{2} & {\boldsymbol{A}}{\boldsymbol{\Delta }}{\boldsymbol{\theta }}{\boldsymbol{\Delta }}{\boldsymbol{\phi }}\end{array}\right],\end{eqnarray} \tag{ 38 }$$

where ${\boldsymbol{\Delta }}{\boldsymbol{\theta }}$ and ${\boldsymbol{\Delta }}{\boldsymbol{\phi }}$ are diagonal matrices giving the difference between the detector pair's pointing and the nearest HEALPix pixel center.

A differential beam generating operator is applied to the timestreams to create differential beam timestreams. For example, the differential gain timestream is just the beam convolved temperature field:

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{\delta g}={\delta }_{g}{{\boldsymbol{A}}}^{\prime }{\boldsymbol{\Theta }},\end{eqnarray} \tag{ 39 }$$

where the fit coefficient for the gain mismatch is ${\delta }_{g}$. The differential pointing components are found from the first derivatives of the temperature field with respect to the focal plane coordinates, x and y:

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{\delta x}={\delta }_{x}{{\boldsymbol{\nabla }}}_{x}{{\boldsymbol{A}}}^{\prime }{\boldsymbol{\Theta }}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{\delta y}={\delta }_{y}{{\boldsymbol{\nabla }}}_{y}{{\boldsymbol{A}}}^{\prime }{\boldsymbol{\Theta }},\end{eqnarray} \tag{ 41 }$$

where ${\delta }_{x}$ and ${\delta }_{y}$ are the differential beam coefficients and ${{\boldsymbol{\nabla }}}_{x}$ and ${{\boldsymbol{\nabla }}}_{y}$ are partial differential operators with respect to the focal plane coordinates. Further details of this calculation and derivations for other beam modes are discussed in Appendix C of Bicep2 Collaboration III (2015).

The differential beam timestreams are transformed into maps analogously to Equation (35), creating a pairmap template, ${\tilde{{\boldsymbol{ \mathcal T }}}}_{j}$, for each differential beam mode, $j$. The template pairmaps for each scanset, ${\mathbb{S}}$, are then coadded over phases. For instance, the template pairmap for differential gain, is ${\tilde{{\boldsymbol{ \mathcal T }}}}_{1}$:

$$\begin{eqnarray}{\tilde{{\boldsymbol{ \mathcal T }}}}_{1} & = & \displaystyle \sum _{{\mathbb{S}}}^{{\mathbb{S}}\in \mathrm{phase}}\left(\left[\begin{array}{cc}{\boldsymbol{\Lambda }} & 0\\ 0 & {\boldsymbol{\Lambda }}\end{array}\right]\left[\begin{array}{c}{{\boldsymbol{\alpha }}}^{-}\\ {{\boldsymbol{\beta }}}^{-}\end{array}\right]{{\boldsymbol{w}}}^{-}{\boldsymbol{GF}}\right.\\ & & \times {\left.\left[\begin{array}{cc}{{\boldsymbol{\alpha }}}^{-} & {{\boldsymbol{\beta }}}^{-}\end{array}\right]\left[\begin{array}{cc}{{\boldsymbol{A}}}^{\prime } & 0\\ 0 & {{\boldsymbol{A}}}^{\prime }\end{array}\right]\left[\begin{array}{c}{\boldsymbol{\Theta }}\\ {\boldsymbol{\Theta }}\end{array}\right]\right)}_{{\mathbb{S}}}.\end{eqnarray} \tag{ 42 }$$

A matrix performing weighted linear least-squares regression against real pairmaps produces the fitted coefficients for each of the differential beam modes, ${\boldsymbol{c}}\equiv [{\delta }_{g},{\delta }_{x},{\delta }_{y},\ldots ]$:

$$\begin{eqnarray}&&{\boldsymbol{c}}={({\tilde{{\boldsymbol{ \mathcal T }}}}^{\top }{({{\boldsymbol{ \mathcal W }}}^{-})}^{-1}\tilde{{\boldsymbol{ \mathcal T }}})}^{-1}{\tilde{{\boldsymbol{ \mathcal T }}}}^{\top }{({{\boldsymbol{ \mathcal W }}}^{-})}^{-1}\left[\begin{array}{c}{{\boldsymbol{m}}}_{{{\boldsymbol{\alpha }}}^{-}}\\ {{\boldsymbol{m}}}_{{{\boldsymbol{\beta }}}^{-}}\end{array}\right],\end{eqnarray} \tag{ 43 }$$

where $\left[\begin{array}{c}{{\boldsymbol{m}}}_{{{\boldsymbol{\alpha }}}^{-}}\\ {{\boldsymbol{m}}}_{{{\boldsymbol{\beta }}}^{-}}\end{array}\right]$ is the real data pairmap coadded over a phase, $\tilde{{\boldsymbol{ \mathcal T }}}$ is a vector of pairmap templates, and ${{\boldsymbol{ \mathcal W }}}^{-}$ is the pair difference weight map, created from the weight matrix according to:

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal W }}}^{-}=\displaystyle \sum _{{\mathbb{S}}}^{{\mathbb{S}}\in \mathrm{phase}}{\left(\left[\begin{array}{c}{\boldsymbol{\Lambda }}{{\boldsymbol{\alpha }}}^{-}{{\boldsymbol{\alpha }}}^{-}{{\boldsymbol{w}}}^{-}{{\boldsymbol{\Lambda }}}^{\top }\\ {\boldsymbol{\Lambda }}{{\boldsymbol{\beta }}}^{-}{{\boldsymbol{\beta }}}^{-}{{\boldsymbol{w}}}^{-}{{\boldsymbol{\Lambda }}}^{\top }\end{array}\right]\right)}_{{\mathbb{S}}}.\end{eqnarray} \tag{ 44 }$$

The pairmap templates weighted by ${\boldsymbol{c}}$ are then subtracted from the real data pairmap:

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{{\boldsymbol{m}}}}_{{{\boldsymbol{\alpha }}}^{-}}\\ {\tilde{{\boldsymbol{m}}}}_{{{\boldsymbol{\beta }}}^{-}}\end{array}\right]=\left[\begin{array}{c}{{\boldsymbol{m}}}_{{{\boldsymbol{\alpha }}}^{-}}\\ {{\boldsymbol{m}}}_{{{\boldsymbol{\beta }}}^{-}}\end{array}\right]-\tilde{{\boldsymbol{ \mathcal T }}}{\boldsymbol{c}}.\end{eqnarray} \tag{ 45 }$$

This process takes the form of a matrix operator that includes each of the beam systematics, giving the deprojection matrix:

$$\begin{eqnarray}&&{\boldsymbol{D}}\equiv {\boldsymbol{I}}-\tilde{{\boldsymbol{ \mathcal T }}}{({\tilde{{\boldsymbol{ \mathcal T }}}}^{\top }{({{\boldsymbol{ \mathcal W }}}^{-})}^{-1}\tilde{{\boldsymbol{ \mathcal T }}})}^{-1}{\tilde{{\boldsymbol{ \mathcal T }}}}^{\top }{({{\boldsymbol{ \mathcal W }}}^{-})}^{-1}.\end{eqnarray} \tag{ 46 }$$

Deprojected pairmaps are then found according to:

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{{\boldsymbol{m}}}}_{{{\boldsymbol{\alpha }}}^{-}}\\ {\tilde{{\boldsymbol{m}}}}_{{{\boldsymbol{\beta }}}^{-}}\end{array}\right]={\boldsymbol{D}}\displaystyle \sum _{{\mathbb{S}}}^{{\mathbb{S}}\in {phase}}{{\boldsymbol{ \mathcal Q }}}_{{\mathbb{S}}}\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right].\end{eqnarray} \tag{ 47 }$$

The regression in Equation (43) operates simultaneously over all the modes to be deprojected. Because the templates for different modes are not in general orthogonal, the coefficient for each mode depends on the full set of modes. Therefore, the subtraction in Equation (46) must include the same mode list used in the regression in Equation (43). If the regression included more modes than the subtraction step, the regression would have extra degrees of freedom. This could result in incomplete removal of leakage signal. We avoid this possibility by deferring the regression step until immediately before the subtraction step, and explicitly using the same mode selection for both.

As in the standard pipeline, we deproject for each detector pair, after coadding scanset to phases. To reduce the computational demands, the matrix deprojection pairmaps have additionally been coadded over scan direction, whereas the standard pipeline performs regression separately for left going and right going scans. This is the only difference between simulations run with the standard pipeline and those calculated from the observation matrix and leads to a negligible difference, see Figures 11 and 12.

The deprojection matrix made for a phase is less sparse than one made for a scanset because over the course of a phase a particular pair will observe a larger range of elevation than it would in a scanset. The filled elements of the matrix ${\boldsymbol{D}}$, for one pair across one phase, is shown in Figure 8.

Zoom In Zoom Out Reset image size

An observed temperature map, $\tilde{{\boldsymbol{T}}}^{\prime}$, can be found by summing the pair sum pairmaps temporally over scansets (${\mathbb{S}}$) and over detector pairs (${\mathbb{P}}$):

$$\begin{eqnarray}&&\tilde{{\boldsymbol{T}}}^{\prime} =\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}{\rm{\Lambda }}{{\boldsymbol{w}}}^{+}{\boldsymbol{GFA}}[{\boldsymbol{T}}].\end{eqnarray} \tag{ 48 }$$

The matrix performing this transformation is defined as ${\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{TT}}}$, where the prime indicates the apodization comes from the inverse variance of the pair sum timestream, ${{\boldsymbol{w}}}^{+}$. The final apodization is applied in Section 4.13.

The transformation from pair difference pairmaps to $Q,U$ maps depends on the detector orientations during the observations. This transformation relies on an inversion of a 2 × 2 detector orientation matrix. We will now derive the matrix that performs this transformation.

Ignoring filtering, the pair difference timestream is found using the timestream forming matrix, ${A}_{{tj}}$:

$$\begin{eqnarray}&&{d}_{t}=\displaystyle \frac{1}{2}({\tau }_{t}^{A}-{\tau }_{t}^{B})={\alpha }_{t}^{-}{A}_{{tj}}{Q}_{j}+{\beta }_{t}^{-}{A}_{{tj}}{U}_{j}.\end{eqnarray} \tag{ 49 }$$

Forming linear combinations of the pair difference timestream,

$$\begin{eqnarray}&&\left[\begin{array}{c}{\alpha }_{t}^{-}{d}_{t}\\ {\beta }_{t}^{-}{d}_{t}\end{array}\right]=\left[\begin{array}{cc}{\alpha }_{t}^{-}{\alpha }_{t}^{-}{A}_{{tj}} & {\alpha }_{t}^{-}{\beta }_{t}^{-}{A}_{{tj}}\\ {\alpha }_{t}^{-}{\beta }_{t}^{-}{A}_{{tj}} & {\beta }_{t}^{-}{\beta }_{t}^{-}{A}_{{tj}}\end{array}\right]\left[\begin{array}{c}{Q}_{j}\\ {U}_{j}\end{array}\right],\end{eqnarray} \tag{ 50 }$$

and applying the pointing matrix, ${\boldsymbol{\Lambda }}$, the vectors ${{\boldsymbol{\alpha }}}^{-}{\boldsymbol{d}}$ and ${{\boldsymbol{\beta }}}^{-}{\boldsymbol{d}}$ are binned into map pixels, i. At this point we coadd over scansets and detector pairs, and apply a weighting, ${{\boldsymbol{w}}}^{-}$, equal to the inverse of the variance of the timestreams during a scanset:

$$\begin{eqnarray}&&\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}\left[\begin{array}{c}{{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\alpha }_{t}^{-}{d}_{t}\\ {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\beta }_{t}^{-}{d}_{t}\end{array}\right] & = & \left(\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}\left[\begin{array}{ll}{{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\alpha }_{t}^{-}{\alpha }_{t}^{-}{A}_{{tj}} & {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\alpha }_{t}^{-}{\beta }_{t}^{-}{A}_{{tj}}\\ {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\alpha }_{t}^{-}{\beta }_{t}^{-}{A}_{{tj}} & {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\beta }_{t}^{-}{\beta }_{t}^{-}{A}_{{tj}}\end{array}\right]\right)\\ & & \times \left[\begin{array}{c}{Q}_{j}\\ {U}_{j}\end{array}\right].\end{eqnarray} \tag{ 51 }$$

We invert the matrix on the right hand side of Equation (51) to compute a matrix that generates Q and U maps:

$$\begin{eqnarray}\left[\begin{array}{cc}{e}_{i} & {f}_{i}\\ {f}_{i} & {g}_{i}\end{array}\right]\equiv {\left(\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}\left[\begin{array}{cc}{{\rm{\Lambda }}}_{{it}}{{w}_{t}}^{-}{\alpha }_{t}^{-}{\alpha }_{t}^{-}{{\rm{\Lambda }}}_{{ti}} & {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\alpha }_{t}^{-}{\beta }_{t}^{-}{{\rm{\Lambda }}}_{{ti}}\\ {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\alpha }_{t}^{-}{\beta }_{t}^{-}{{\rm{\Lambda }}}_{{ti}} & {{\rm{\Lambda }}}_{{it}}{w}_{t}^{-}{\beta }_{t}^{-}{\beta }_{t}^{-}{{\rm{\Lambda }}}_{{ti}}\end{array}\right]\right)}^{-1},\end{eqnarray} \tag{ 52 }$$

where the t index has been summed over to find each of the elements in the 2 × 2 matrix on the right hand side and ${A}_{{tj}}$ has been replaced by ${{\rm{\Lambda }}}_{{ti}}$ so the equation now determines the $Q,U$ values in the observed map, ${Q}_{i},{U}_{i}$. There is one 2 × 2 matrix inversion performed for each pixel, i, in the observed map. In other words, one value of ${e}_{i}$, ${f}_{i}$, and ${g}_{i}$ is computed for each pixel in the observed map, and filled into the i-th diagonal element of ${\boldsymbol{e}}$, ${\boldsymbol{f}}$ and ${\boldsymbol{g}}$.

The pairmaps ${\tilde{{\boldsymbol{m}}}}_{{{\boldsymbol{\alpha }}}^{-}}$ and ${\tilde{{\boldsymbol{m}}}}_{{{\boldsymbol{\beta }}}^{-}}$ are transformed into Stokes $Q,U$ by multiplying by $\left[\begin{array}{cc}{\boldsymbol{e}} & {\boldsymbol{f}}\\ {\boldsymbol{f}} & {\boldsymbol{g}}\end{array}\right]$. Observed $Q,U$ maps are found according to:

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{{\boldsymbol{Q}}}}^{\prime }\\ \tilde{{\boldsymbol{U}}}^{\prime} \end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{e}} & {\boldsymbol{f}}\\ {\boldsymbol{f}} & {\boldsymbol{g}}\end{array}\right]\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}{{\boldsymbol{ \mathcal Q }}}_{{\mathbb{P}},{\mathbb{S}}}\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right],\end{eqnarray} \tag{ 53 }$$

where ${\boldsymbol{ \mathcal Q }}$ was defined in Equation (36).

If the sum of ${\boldsymbol{ \mathcal Q }}$ matrices could be inverted, it would be possible to use this inverse to recover an unbiased estimate of the original Q and U. However, ${\boldsymbol{ \mathcal Q }}$ is singular because it includes polynomial filtering and scan-synchronous signal subtraction, which completely remove some modes that were present in the original maps. We therefore use instead the matrix defined by the matrix inversion in Equation (52), which does not include these filtering operations. Even the inversion in Equation (52) is singular unless the coadded data contains observations at multiple detector angles, ${{\rm{\Psi }}}_{t}$. Observations at multiple detector orientations are made through deck rotations or by coadding over receivers in different orientations. As described in Section 4.2, deck rotations occur between phases, so coadding over phases makes the matrix invertible.

Including deprojection, Equation (53) becomes:

$$\begin{eqnarray}&&\left[\begin{array}{c}\tilde{{\boldsymbol{Q}}^{\prime} }\\ \tilde{{\boldsymbol{U}}^{\prime} }\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{e}} & {\boldsymbol{f}}\\ {\boldsymbol{f}} & {\boldsymbol{g}}\end{array}\right]\displaystyle \sum _{{\mathbb{P}}}\;{\left({\boldsymbol{D}}\displaystyle \sum _{{\mathbb{S}}}^{{\mathbb{S}}\in \mathrm{phase}}{{\boldsymbol{ \mathcal Q }}}_{{\mathbb{S}}}\right)}_{{\mathbb{P}}}\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right].\end{eqnarray} \tag{ 54 }$$

This represents the entire $Q,U$ map making process for signal simulations: from input maps to observed maps, including filtering operations. It can be summarized as:

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{{\boldsymbol{Q}}}}^{\prime }\\ {\tilde{{\boldsymbol{U}}}}^{\prime }\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{QQ}}} & {\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{QU}}}\\ {\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{UQ}}} & {\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{UU}}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right].\end{eqnarray} \tag{ 55 }$$

As constructed, the matrix, ${\boldsymbol{R}}^{\prime}$, contains an apodization based on the inverse variance of the timestreams, ${{\boldsymbol{w}}}^{+}$ and ${{\boldsymbol{w}}}^{-}$. We can, however, choose to remove this apodization, producing maps with equal weight across the field in units of  μK. We construct the quantities:

$$\begin{eqnarray}&&{{\boldsymbol{W}}}^{+}=\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}{({\boldsymbol{\Lambda }}{{\boldsymbol{w}}}^{+}{{\boldsymbol{\Lambda }}}^{\top })}_{{\mathbb{P}},{\mathbb{S}}},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}{{\boldsymbol{W}}}^{-}=\left[\begin{array}{cc}{\boldsymbol{e}} & {\boldsymbol{f}}\\ {\boldsymbol{f}} & {\boldsymbol{g}}\end{array}\right]\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}\;{\left(\left[\begin{array}{c}{\boldsymbol{\Lambda }}{{\boldsymbol{\alpha }}}^{-}{{\boldsymbol{\alpha }}}^{-}{{\boldsymbol{w}}}^{-}{{\boldsymbol{\Lambda }}}^{\top }\\ {\boldsymbol{\Lambda }}{{\boldsymbol{\beta }}}^{-}{{\boldsymbol{\beta }}}^{-}{{\boldsymbol{w}}}^{-}{{\boldsymbol{\Lambda }}}^{\top }\end{array}\right]\right)}_{{\mathbb{P}},{\mathbb{S}}},\end{eqnarray} \tag{ 57 }$$

and use these to remove the apodization from the observation matrix, solving for the non-apodized observation matrix, ${\boldsymbol{ \mathcal R }}$:

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal R }}}_{{TT}}={({{\boldsymbol{W}}}^{+})}^{-1}{\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{TT}}}\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}\left[\begin{array}{cc}{{\boldsymbol{ \mathcal R }}}_{{QQ}} & {{\boldsymbol{ \mathcal R }}}_{{QU}}\\ {{\boldsymbol{ \mathcal R }}}_{{UQ}} & {{\boldsymbol{ \mathcal R }}}_{{UU}}\end{array}\right]={({{\boldsymbol{W}}}^{-})}^{-1}\left[\begin{array}{cc}{\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{QQ}}} & {\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{QU}}}\\ {\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{UQ}}} & {\boldsymbol{R}}{^{\prime} }_{{\boldsymbol{UU}}}\end{array}\right]\end{eqnarray} \tag{ 59 }$$

Using the non-apodized observation matrix of Section 4.12, we can create an observation matix with an arbitrary apodization, ${\boldsymbol{Z}}$. The matrix ${\boldsymbol{R}}$ is constructed as follows:

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{{\boldsymbol{TT}}}={\boldsymbol{Z}}{{\boldsymbol{ \mathcal R }}}_{{TT}}\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}\left[\begin{array}{cc}{{\boldsymbol{R}}}_{{\boldsymbol{QQ}}} & {{\boldsymbol{R}}}_{{\boldsymbol{QU}}}\\ {{\boldsymbol{R}}}_{{\boldsymbol{UQ}}} & {{\boldsymbol{R}}}_{{\boldsymbol{UU}}}\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{Z}} & 0\\ 0 & {\boldsymbol{Z}}\end{array}\right]\left[\begin{array}{cc}{{\boldsymbol{ \mathcal R }}}_{{QQ}} & {{\boldsymbol{ \mathcal R }}}_{{QU}}\\ {{\boldsymbol{ \mathcal R }}}_{{UQ}} & {{\boldsymbol{ \mathcal R }}}_{{UU}}\end{array}\right]\end{eqnarray} \tag{ 61 }$$

A sensible choice for the apodization, ${\boldsymbol{Z}}$, may be the inverse variance mask we removed in Section 4.12, or a smoothed version thereof. However, there is freedom to choose any apodization at this point, and this may prove useful in joint analyses with other experiments, where the analysis combines maps with low noise regions in slightly different regions of the sky.

We have constructed a matrix, ${\boldsymbol{R}}$, which performs the linear operations of polynomial filtering, scan-synchronous signal subtraction, deprojection, weighting and pointing. ${\boldsymbol{R}}$ has dimensions $(3{\tilde{n}}_{p},3{n}_{p})$ where ${\tilde{n}}_{p}$ is the number of pixels in the Bicep map and n_p is the number of pixels in the input HEALPix map.

Using the observation matrix, the entire process of generating a signal simulation from an input map is:

$$\begin{eqnarray}&&\left[\begin{array}{c}\tilde{{\boldsymbol{T}}}\\ \tilde{{\boldsymbol{Q}}}\\ \tilde{{\boldsymbol{U}}}\end{array}\right]=\left[\begin{array}{ccc}{{\boldsymbol{R}}}_{{\boldsymbol{TT}}} & 0 & 0\\ 0 & {{\boldsymbol{R}}}_{{\boldsymbol{QQ}}} & {{\boldsymbol{R}}}_{{\boldsymbol{QU}}}\\ 0 & {{\boldsymbol{R}}}_{{\boldsymbol{UQ}}} & {{\boldsymbol{R}}}_{{\boldsymbol{UU}}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{T}}\\ {\boldsymbol{Q}}\\ {\boldsymbol{U}}\end{array}\right].\end{eqnarray} \tag{ 62 }$$

Here the off diagonal terms, ${{\boldsymbol{R}}}_{{\boldsymbol{QU}}}$ and ${{\boldsymbol{R}}}_{{\boldsymbol{UQ}}}$, exist because the filtering operations are performed on pair difference timestreams, which are a combination of Q and U.

The deprojection operator contains regression against a temperature map that must be chosen before constructing the observation matrix. The deprojection operator is a linear filtering operation, and it only removes beam systematics arising from one particular temperature field. One could in principle apply the deprojection operator to $Q,U$ maps corresponding to a different temperature field. The operator would remove the same modes from the polarization field, but these modes would not correspond to those which had been mixed between T and $Q,U$ by beam systematics (or TE correlation).

When constructing an ensemble of E-mode realizations for use in Monte Carlo power spectrum analysis, the TE correlation and the fixed temperature sky force us to build constrained realizations. The ensemble of simulations all contain identical temperature fields, so we cannot use them for analysis of temperature, which is acceptable because the focus of our analysis is polarization. The ensemble does contain different realizations of $Q,U$, constrained for the given temperature field, and these can be used in Monte Carlo analysis of polarization. The details of constructing these constrained input maps is the subject of the Appendix.

Because the construction of ${{\boldsymbol{R}}}_{{\boldsymbol{QQ}}}$, ${{\boldsymbol{R}}}_{{\boldsymbol{QU}}}$, and ${{\boldsymbol{R}}}_{{\boldsymbol{UU}}}$ depends on a fixed temperature field, the deprojection templates can be thought of as numerical constants. ${{\boldsymbol{R}}}_{{\boldsymbol{TT}}}$ performs a separate filtering on the temperature field that is largely decoupled from the filtering of Q and U. To include systematics that leak temperature to polarization, the terms ${{\boldsymbol{R}}}_{{\boldsymbol{TQ}}}$ and ${{\boldsymbol{R}}}_{{\boldsymbol{TU}}}$ would in principle need to be non zero. However, so long as the leakage corresponded to modes being removed by the deprojection matrix ${\boldsymbol{D}}$ the the deprojection elements in the ${{\boldsymbol{R}}}_{{\boldsymbol{QQ}}}$, ${{\boldsymbol{R}}}_{{\boldsymbol{QU}}}$, and ${{\boldsymbol{R}}}_{{\boldsymbol{UU}}}$ blocks would ensure that the output $Q,U$ maps were identical.

Although the nominal dimensions of ${\boldsymbol{R}}$ are large, our constant elevation scan strategy means that ${\boldsymbol{R}}$ is only filled for pixels at roughly the same declination. This means that ${\boldsymbol{R}}$ is largely sparse, as shown in Figure 9.

Zoom In Zoom Out Reset image size

Some intuition about the operations the observation matrix performs can be gained by plotting a column of the matrix reshaped as maps—see Figure 10. The column chosen in this case corresponds to a central pixel in the observed field. It shows how Q and U values in the observed map are sourced from a Q pixel in the HEALPix map. The bright pixel in the Q observed map corresponds to the location of the input Q. The effects of polynomial and scan-synchronous signal subtraction are visible to the left and right of the bright Q pixel. These two types of filtering are performed on scansets and are therefore confined to a row of pixels. Deprojection operates on phases, creating the effects seen at other declinations. Because all of these filtering operations are performed on pair difference data, which contains linear combinations of Q and U, signal in the observed U map can be created by signal in the input HEALPix Q map. This is why the U map in Figure 10 is non-zero.

Zoom In Zoom Out Reset image size

We can form observed maps from the real timestreams using the matrices constructed above:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{T}}}}_{\mathrm{real}}={\boldsymbol{Z}}{({{\boldsymbol{W}}}^{+})}^{-1}\displaystyle \sum _{{\mathbb{P}},{\mathbb{S}}}{({\boldsymbol{\Lambda }}{{\boldsymbol{w}}}^{+}{\boldsymbol{GFs}})}_{{\mathbb{P}},{\mathbb{S}}}\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{{\boldsymbol{Q}}}}_{\mathrm{real}}\\ {\tilde{{\boldsymbol{U}}}}_{\mathrm{real}}\end{array}\right] & = & \left[\begin{array}{cc}{\boldsymbol{Z}} & 0\\ 0 & {\boldsymbol{Z}}\end{array}\right]{({{\boldsymbol{W}}}^{-})}^{-1}\left[\begin{array}{cc}{\boldsymbol{e}} & {\boldsymbol{f}}\\ {\boldsymbol{f}} & {\boldsymbol{g}}\end{array}\right]\\ & & \times \displaystyle \sum _{{\mathbb{P}}}{\left({\boldsymbol{D}}\displaystyle \sum _{{\mathbb{S}}}^{{\mathbb{S}}\in \mathrm{phase}}{\left(\left[\begin{array}{cc}{\boldsymbol{\Lambda }} & 0\\ 0 & {\boldsymbol{\Lambda }}\end{array}\right]\left[\begin{array}{c}{{\boldsymbol{\alpha }}}^{-}\\ {{\boldsymbol{\beta }}}^{-}\end{array}\right]{{\boldsymbol{w}}}^{-}{\boldsymbol{GFd}}\right)}_{{\mathbb{S}}}\right)}_{{\mathbb{P}}}.\end{eqnarray} \tag{ 64 }$$

It is important to note that the exact same matrices are used to process the real data in Equation (64) as are used to construct the simulated maps in Equation (35).

The matrix formalism described above is self-contained and complete in the sense that it contains the tools necessary to create real data maps and simulated maps.

We demand that the map making and filtering operations be identical between the standard pipeline and the observation matrix. It is straightforward to test this equivalence: simulated maps run through the standard pipeline must be identical to the maps found with the observation matrix. Figures 11 and 12 show that the two match quite well, within a few percent over the multipoles of $50\lt l\lt 350$. The lack of a perfect match is due to the difference in deprojection timescale and because the standard pipeline uses N_side = 2048 HEALPix input maps that are Taylor interpolated, whereas the observation matrix uses N_side = 512 input HEALPix maps with nearest neighbor interpolation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A pixel on the sky at location i, has values of the Stokes parameters:

$$\begin{eqnarray}&&{x}_{i}\equiv \left(\begin{array}{c}{T}_{i}\\ {Q}_{i}\\ {U}_{i}\end{array}\right).\end{eqnarray} \tag{ 65 }$$

The 3 × 3 pixel-pixel covariance between two locations on the sky, i and j, is given by:

$$\begin{eqnarray}&&{{\boldsymbol{C}}}_{i,j}\equiv \langle {x}_{i}{x}_{j}\rangle ={\mathbb{R}}{(\alpha ){\boldsymbol{M}}({r}_{i}\cdot {r}_{j}){\mathbb{R}}(\alpha )}^{\top }.\end{eqnarray} \tag{ 66 }$$

The covariance matrix, ${\boldsymbol{M}}$, is defined with the $Q,U$ convention referenced to the great circle connecting the two points, $i,j$. For a particular spectrum ${\boldsymbol{M}}$ depends only on the dot product between the pixels, ${r}_{i}\cdot {r}_{j}$. ${\boldsymbol{M}}$ contains nine symmetric sub-matrices:

$$\begin{eqnarray}{\boldsymbol{M}}({r}_{i}\cdot {r}_{j})=\left(\begin{array}{ccc}\langle {T}_{i}{T}_{j}\rangle & \langle {T}_{i}{Q}_{j}\rangle & \langle {T}_{i}{U}_{j}\rangle \\ \langle {Q}_{i}{T}_{j}\rangle & \langle {Q}_{i}{Q}_{j}\rangle & \langle {Q}_{i}{U}_{j}\rangle \\ \langle {U}_{i}{T}_{j}\rangle & \langle {U}_{i}{Q}_{j}\rangle & \langle {U}_{i}{U}_{j}\rangle \end{array}\right).\end{eqnarray} \tag{ 67 }$$

The 3 × 3 matrix, ${\mathbb{R}}$, is applied to rotate from this local reference frame to a global frame where $Q,U$ are referenced to the North–South meridians. The angle between the great circle connecting any two points and the global frame is given by the parameter $\alpha$

$$\begin{eqnarray}{\mathbb{R}}(\alpha )=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & \mathrm{cos}2\alpha & \mathrm{sin}2\alpha \\ 0 & -\mathrm{sin}2\alpha & \mathrm{cos}2\alpha \end{array}\right).\end{eqnarray} \tag{ 68 }$$

Changing the sign of α allows us to change the polarization convention from IAU to HEALPix (U to $-U$), see Hamaker & Bregman (1996). We have chosen to use the IAU convention for Bicep2 and Keck Array covariance matrices.

The true sky pixel-pixel signal covariance matrix for the Stokes $Q,U$ parameters is derived in Kamionkowski et al. (1997), Zaldarriaga (1998). To calculate the covariances, we follow some of the suggestions in Appendix of Tegmark & de Oliveira-Costa (2001).

We use the HEALPix ring pixelization, which allows the covariance to be calculated simultaneously for all pixels at a particular latitude that are separated by the same distance.²³ This shortcut is exploited by simultaneously calculating all equidistant pixels for rows in the map that have the same latitude to within 1 × 10⁻⁴ degrees. This approximation is much smaller than the ∼7 arcminute pixels in N_side = 512 maps, and the rounding error has been found to be insignificant.

The observed signal covariance matrix contains the pixel-pixel covariance in the observed Bicep pixelized maps. Theoretically, modifying the true sky signal covariance matrix is simple: take the unobserved signal covariance of ${\boldsymbol{C}}$ of Section 5.1 and the observation matrix ${\boldsymbol{R}}$ from Section 4 and form the product:

$$\begin{eqnarray}&&\tilde{{\boldsymbol{C}}}={{\boldsymbol{RCR}}}^{\top }.\end{eqnarray} \tag{ 69 }$$

This equation results in a symmetric positive definite matrix, which is rank deficient because of the filtering steps in the observing process.

Unfortunately, performing this multiplication is computationally demanding: The input ${\boldsymbol{C}}$ is a square matrix, with 3 × 111,593 elements on a side, corresponding to the elements of T, Q, and U. To reduce the memory requirements of the calculation, we divide the covariance matrix, ${\boldsymbol{C}}$, into row subsets and calculate in parallel. Once a row subset is calculated, the observation matrix is immediately applied to transform the HEALPix covariance to the observed map covariance, which reduces the dimensions of the covariance to the 23,600 pixels of the observed maps.

The covariance and observed covariance should both be symmetric, which provides a good check on our math. Usually the output is slightly (fractionally, $\sim 1/{10}^{7}$) non-symmetric due to rounding errors in the multiplication, and we force the final matrix to be symmetric to numerical precision by averaging across the diagonal before moving to the next steps, since symmetric matrices often allow the use of faster algorithms.

A row of the observed covariance matrix can be reshaped into a map, which reveals the structure of the covariance for a particular pixel, see Figure 13.

Zoom In Zoom Out Reset image size

As written, Equation (17) is not solvable: ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{B}}}$ has a null space that is the set of pure E-modes. Similarly, the space of pure B-modes is the null space of ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{E}}}$. By adding the identity matrix multiplied by a constant, ${\sigma }^{2}{\boldsymbol{I}}$, to the covariance matrices we regularize the problem to find approximate solutions and eliminate the null spaces:

$$\begin{eqnarray}&&({\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{B}}}+{\sigma }^{2}{\boldsymbol{I}}){{\boldsymbol{x}}}_{{\boldsymbol{i}}}={\lambda }_{i}({\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{E}}}+{\sigma }^{2}{\boldsymbol{I}}){{\boldsymbol{x}}}_{{\boldsymbol{i}}}.\end{eqnarray} \tag{ 70 }$$

The amplitude of ${\sigma }^{2}$ sets the relative magnitude of the ambiguous mode eigenvalues versus the pure E and pure B-mode eigenvalues. The eigenvalues are shown in Figure 14. In our analysis, we choose ${\sigma }^{2}$ to be 1/100 the mean of the diagonal elements of the covariance matrices, ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{E}}}$ and ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{B}}}$.

Zoom In Zoom Out Reset image size

In Equation (70), ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{E}}}$ is an observed covariance matrix for $[Q,U]$ that is constructed according to Equations (11) and (69). The input spectrum is set to a steeply red E-mode spectrum, ${C}_{l}^{{EE}}=1/{l}^{2}$, ${C}_{l}^{{BB}}=0$. ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{B}}}$ is the same except for an input spectrum with only B-modes, ${C}_{l}^{{BB}}=1/{l}^{2}$, ${C}_{l}^{{EE}}=0$. The eigenmodes in ${{\boldsymbol{x}}}_{{\boldsymbol{i}}}$ with the largest eigenvalues comprise a set of vectors that span a space of pure B-modes, ${{\boldsymbol{b}}}_{{\boldsymbol{i}}}$. The pure B quality of these vectors can be seen in the fact that the product ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{E}}}{\boldsymbol{x}}$ is much smaller than the product ${\tilde{{\boldsymbol{C}}}}_{{\boldsymbol{B}}}{\boldsymbol{x}}$ . The eigenmodes in ${{\boldsymbol{x}}}_{{\boldsymbol{i}}}$ with the smallest eigenvalues comprise a set of vectors that span a space of pure E-modes, ${{\boldsymbol{e}}}_{{\boldsymbol{i}}}$.

Using the reddened input spectrum $1/{l}^{2}$ causes the magnitude of the eigenvalues to be proportional to the band of multipole, l, that each mode contains. The steepness of the spectrum ensures each mode contains power from a limited range of l. The particular choice of $1/{l}^{2}$ is arbitrary.

A basis constructed from a subset of eigenmodes with large eigenvalues spans a subspace of pure B-modes. We arbitrarily choose the pure B subspace to consist of eigenmodes whose corresponding eigenvalues are the largest 1/4 of the set. We find modes contained in this subset adequate for preserving power up to $l\sim 700$. The pure E subspace is similarly constructed with eigenmodes corresponding to the 1/4 smallest eigenvalues. Figure 14 shows the sorted eigenvalues, and Figure 15 shows four eigenmodes of the Bicep2 observed covariance.

Zoom In Zoom Out Reset image size

Using the set of pure E and pure B-mode basis vectors, two projection matrices are constructed from the outer products:

$$\begin{eqnarray}{{\boldsymbol{\Pi }}}_{E} & = & \displaystyle \sum _{i}{{\boldsymbol{e}}}_{{\boldsymbol{i}}}{{\boldsymbol{e}}}_{{\boldsymbol{i}}}^{\top }\\ {{\boldsymbol{\Pi }}}_{B} & = & \displaystyle \sum _{i}{{\boldsymbol{b}}}_{{\boldsymbol{i}}}{{\boldsymbol{b}}}_{{\boldsymbol{i}}}^{\top },\end{eqnarray} \tag{ 71 }$$

which we call the purification matrices for pure E and pure B. Operating the purification matrices on an input map projects onto the space of pure E and pure B:

$$\begin{eqnarray}{\tilde{{\boldsymbol{m}}}}_{\mathrm{pureE}} & = & {{\boldsymbol{\Pi }}}_{E}\tilde{{\boldsymbol{m}}}\\ {\tilde{{\boldsymbol{m}}}}_{\mathrm{pureB}} & = & {{\boldsymbol{\Pi }}}_{B}\tilde{{\boldsymbol{m}}}.\end{eqnarray} \tag{ 72 }$$

From the construction of the pure B basis, ${{\boldsymbol{b}}}_{{\boldsymbol{i}}}$, one can see that the vector ${\tilde{{\boldsymbol{m}}}}_{\mathrm{pureB}}$ vanishes for arbitrary input containing only E-modes: $\tilde{{\boldsymbol{m}}}={\boldsymbol{R}}\cdot \sum {a}_{{lm}}^{E}{{\boldsymbol{Y}}}_{{lm}}^{E}$, as desired for a purified B map.

The Bicep2 and Keck Array analysis pipeline contains the following attributes that can leak E-modes to B-modes: partial sky coverage, timestream filtering (including deprojection), and choice of map projection+estimator. A simulation demonstrating the leaked B-mode maps for each of these effects is shown in Figure 16. These maps are created by applying the standard E and B estimator in Fourier space and then an inversion back to map space. The observations and data reduction produce three classes of E/B leakage:

• 1.  
  Apodization: The first obvious deviation from the ideal full sky map is the partial sky coverage of Bicep2 and Keck Array maps. Once a boundary is imposed E/B leakage is created. Map boundary effects are reduced by an apodization window which tapers the maps smoothly to zero near the edges. Use of apodization windows is common practice in any Fourier transform analysis of finite regions to prevent "ringing" near boundaries. For small regions of sky, effects of the map boundary dominate the leakage even after apodization is applied. The apodization window used in Bicep2 and the Keck Array are modified inverse variance maps. We apply a smoothing Gaussian with a width of $\sigma =0\buildrel{\circ}\over{.} 5$ to the inverse variance map, and, for the combined analyses of Bicep2 and the Keck Array, we use a geometric mean of the individual experiments' apodization maps.
• 2.  
  Map projection: The total extent of the Bicep2 and Keck Array maps is about 50° on the sky in the direction of R.A., over a declination of $-70\lt \delta \lt 45$. The chosen projection is a simple rectangular grid of pixels in R.A. and decl.. Taking standard discrete Fourier transforms of such maps results in significant E/B leakage. While we note that other map projections will have significantly lower distortion, such effects will be present for all projections when subjected to Fourier transform. However, note that the ${\boldsymbol{A}}$ and ${\boldsymbol{\Lambda }}$ matrices together fully encode the mapping from the underlying curved sky to the flat sky of the observed map pixels, and therefore so does the observing matrix ${\boldsymbol{R}}$ derived from these.
• 3.  
  Linear filtering effects: There are three main analytic filters applied in the standard pipeline: polynomial filter, scan-synchronous subtraction, and deprojection. With respect to E to B leakage, all three filters are similar: by removing modes in the $Q,U$ maps, an E-mode can be turned into an observed B-mode. Polynomial filtering and scan-synchronous signal subtraction create comparable leakage, both in amplitude and morphology. Deprojection creates more power in the leaked B maps, and at smaller angular scales, than either polynomial filtering or scan-synchronous signal subtraction.

Zoom In Zoom Out Reset image size

If matrix purification is not used, the sample variance of E/B leakage in the Bicep2 BB power spectrum is comparable to the uncertainty due to instrumental noise. This is clearly highly undesirable, and led to the development of the purification matrix described in this paper. The purification matrix "knows about" all of the E/B mixing effects and how to deal with them.

The effectiveness of the purification matrix given by Equation (72) can be immediately tested by applying the operator to a vector of $[Q,U]$ maps simulated with the standard pipeline. The upper left map of Figure 17 shows an observed map whose input is unlensed-ΛCDM E-modes, and the next two rows in that column show the resulting E and B maps after projection onto the pure E and B spaces. The right column of Figure 17 shows an observed input map with both unlensed-ΛCDM and r = 0.1 projected onto pure E and B-modes.

Zoom In Zoom Out Reset image size

Matrix purification is integrated with the existing Bicep analysis code by applying the purification operator to maps before calculating the power spectra. Since the observation matrix is only used for this purification step, the purification matrix need only work well enough to result in E/B leakage less than the noise level of the experiment. Therefore, it is acceptable to use an approximate observation matrix constructed from a subset of the full observation list, as long as it is representative of the full scan strategy. This shortcut was employed in Keck Array & BICEP2 Collaborations V (2015), but the results shown in this section from Bicep2 are from the full set of observations.

Figure 18 compares the spectra for purified maps to the spectra from maps without purification, and to the spectra found using the improved estimator suggested in Smith (2006). Both the purified maps and unpurified maps use the standard E and B estimator in Fourier space, and we refer to the unpurified maps processed this way as the "normal method." Figure 18 shows the spectra for 200 noiseless unlensed-ΛCDM simulations passed through the three estimators. The leaked power is roughly three orders of magnitudes smaller when using the matrix purification than when using the normal method or Smith estimator. While the Smith estimator improves over the normal estimator by eliminating E/B leakage from apodization, it does not account for spatial filtering, which is a significant source of E/B leakage in the analysis pipeline. The mean of the leaked power is de-biased from the final power spectra, so what matters is the variance of the leaked spectra. Computing the 95% confidence limits based on the variance in each of the three methods, we find that in the absence of B-mode signal or instrumental noise, the matrix estimator achieves a limit on the tensor-to-scalar ratio of $r\lt 8.3\times {10}^{-5}$, while the normal method and Smith estimator achieve limits of $r\lt 0.17$ and $r\lt 0.074$ respectively.

Zoom In Zoom Out Reset image size

Figure 19 shows the spectra from the three estimators for input maps containing only input B-modes at the level of r = 0.1. The spectra for all three estimators show beam roll off at high l. At the lowest l, the filtering prevents large angular scale modes from being measured. For multipoles around $l\sim 100$, the matrix estimator recovers slightly less signal than the other two methods. However, the extra power measured by the other methods near $l\sim 100$ largely comes from the ambiguous modes. On the left of Figure 18, these ambiguous modes are seen as the bump in the normal and Smith method near $l\sim 100$.

Zoom In Zoom Out Reset image size

The total number of degrees of freedom in each band power can be estimated according to the formula:

$$\begin{eqnarray}&&{N}_{l^{\prime} }=2\displaystyle \frac{{({m}_{l^{\prime} })}^{2}}{{\sigma }_{l^{\prime} }^{2}},\end{eqnarray} \tag{ 73 }$$

where ${m}_{l^{\prime} }$ is the mean of the simulations in band power $l^{\prime}$ and ${\sigma }_{l^{\prime} }^{2}$ is the variance of the simulations in the band power. Table 1 shows the number of degrees of freedom for the three estimators for a tensor B-mode. The fewer degrees of freedom at low l for the matrix estimator are consistent with the decrease in recovered power on the left side of Figure 18. The highest l bins in Table 1 also show fewer degrees of freedom in the matrix estimator because the purification matrix includes a limited number of pure eigenmodes, as shown in Figure 14.

Figure 20 shows signal plus noise spectra for a set of 200 unlensed-ΛCDM+noise spectra. The noise simulations are the standard Bicep2 sign flip realizations discussed in BICEP2 Collaboration I (2014). In Figure 20, the mean noise level and leaked BB power are de-biased. The resulting ensemble of simulations is used to construct the errorbars in the final spectra. The tighter distribution of the matrix estimator simulations is a result of the decrease in E/B leakage. Using the matrix estimator results in an improvement in the r limit, for Bicep2 noise level and filtering in the absence of B-mode signal, of about a factor of two over the Smith method. The remaining variance in the BB spectrum of the matrix estimator is instrumental noise.

Zoom In Zoom Out Reset image size

The observation matrix transforms an input HEALPix map into an observed map with a simple matrix multiplication. The speed of the operation facilitates the calculation of analysis transfer functions, which are a necessary component of the pseudo-C_l MASTER algorithm (Hivon et al. 2002).

We start with input maps, ${{\boldsymbol{m}}}_{{\boldsymbol{l}}}$, which are delta functions in a particular multipole. These maps are then observed using the matrix ${\boldsymbol{R}}$. The spectra calculated from these maps represent the response in our analysis pipeline to the input delta function, in a manner conceptually analogous to Green's functions.

Our procedure uses two sets of HEALPix maps, one set corresponding to ${TT}={TE}={EE}=1$ and one set with ${TT}={BB}=1$. The observation matrix is used to create maps for l = 1 through 700, with 100 random realizations for each l. Processing the 140,000 maps would be infeasible without the observation matrix, but using the observation matrix it can be accomplished in a few hours.

7.3.1. Band Power Window Functions

The power spectra of the output maps for a particular l are averaged over the N = 100 realizations. The averaged spectra are used to form a band power window function, ${{\boldsymbol{ \mathcal M }}}_{{ll}^{\prime} }^{{XX}}$, for a particular band power, $l^{\prime}$, which is a function of the input multipole of the delta function, $l$:

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal M }}}_{{ll}^{\prime} }^{{XX}}=\displaystyle \frac{{\displaystyle \sum }_{r=1}^{N}{{\boldsymbol{ \mathcal F }}}_{{ll}^{\prime} }({{\boldsymbol{Rm}}}_{{\boldsymbol{l}}})}{N},\end{eqnarray} \tag{ 74 }$$

where ${{\boldsymbol{ \mathcal F }}}_{{ll}^{\prime} }$ is the analysis pipeline's transformation from map to power spectra²⁴ and ${XX}=\{{TT}\to {TT},{TE}\to {TE}$, ${EE}\to {EE},{EE}\to {BB},{BB}\to {BB}$, ${BB}\to {EE}\}$. When the input maps contain E-modes, and we measure the BB spectra, the result is the ${EE}\to {BB}$ band power window. The use of the purification matrix prevents leakage and makes these band power windows have much lower amplitude than the ${EE}\to {EE}$ or ${BB}\to {BB}$ ones. In Bicep2, it has been standard procedure not to use the E-mode purification matrix, hence the ${BB}\to {EE}$ band power window contains (irrelevant) leakage from B into E.

Calculating the band power window functions in this way accounts for all aspects of our instrument and analysis: beam convolution, sky cut, map projection, polynomial filtering, scan-synchronous signal subtraction, deprojection, and power spectrum estimation.

Figure 21 shows the results of the calculation. Although we bin in annular rings of the two-dimensional power spectra, the filtering operations move power from other l into those bins. This means the band powers are sensitive to a broader range in l than the nominal range of multipoles. The broad shelf in the band power window functions at lower l is due to filtering.

Zoom In Zoom Out Reset image size

Table 2 shows the "nominal" and "measured" centers and edges of the band power bins, where "nominal" refers to the defined range of annular rings in the two-dimensional power spectra and "measured" refers to the center and $\pm 1\sigma$ range found in the end to end calculation discussed in this section. The center is the mean of the bandpower window function and the $\pm 1\sigma$ interval corresponds to the percentiles of the bandpower window function between 16% and 84%.

Table 2.  BB Band Power Widths

	Nominal			Measured
Bin Number	Low	Center	High	Low	Center	High
1	20.0	37.5	55.0	37.0	46.4	54.0
2	55.0	72.5	90.0	59.0	73.4	86.0
3	90.0	107.5	125.0	92.0	107.2	121.0
4	125.0	142.5	160.0	125.0	140.7	157.0
5	160.0	177.5	195.0	158.0	173.7	192.0
6	195.0	212.5	230.0	189.0	205.4	227.0
7	230.0	247.5	265.0	220.0	237.2	262.0
8	265.0	282.5	300.0	253.0	270.2	298.0
9	300.0	317.5	335.0	285.0	302.9	333.0

Note. Nominal and measured centers and edges of the band power bins. The measured values are extracted from the band power window functions shown in Figure 21. The low/high values for the latter are the $\pm 1\sigma$ points.

Download table as:  ASCIITypeset image

7.3.2. Suppression Factor

The integrated area under the curve of each band power window function represents the response of that band power measurement to an input spectrum. We call this set of values the suppression factor, ${S}_{l^{\prime} }$, since they approximate how our analysis pipeline suppresses power.

The suppression factor is plotted in Figure 22. At small angular scales, the suppression factor is dominated by the roll-off of Bicep2's 31 arcminute beam. The measured array average 150 GHz beam function is shown as the dotted line in Figure 22. The "map window function," which includes the finite size of the map and the pixel window function for the ∼025 pixels, is shown as the dashed line.²⁵ At high l, the pixel window function is sub-dominant to the beam window function. At low l the suppression factor is dominated by the timestream filtering effects.

Zoom In Zoom Out Reset image size

Building the observation matrix requires constructing, multiplying, and finally averaging a large number of sparse matrices. Applying the observation matrix to the true sky signal covariance matrix requires a large matrix multiplication to find the product, ${{\boldsymbol{RCR}}}^{\top }$.

Matrix multiplication is the dominant contributor to computation time. We use matrix multiplication routines built in to MATLAB. These routines incorporate a number of optimized algorithms for computing matrix products, including BLAS (Lawson et al. 1979). At this time we have not compared the run time on GPUs with that on CPUs, but are aware of this as a possible avenue for reducing the compute time.

Solving the eigenvalue problem using the MATLAB function eig() takes about 48 hr and 80 GB of RAM for the full Bicep2 observed covariance. For Bicep2, this is small fraction of the total computation time. However, for experiments whose maps contain more pixels, the difficulty of the eigenvalue problem increases. In these cases, the use of distributed memory parallel code may be necessary.

For the Bicep2 results, we used computing resources provided by the Odyssey cluster at Harvard.²⁶ The Odyssey cluster contains 54K CPUs with 190 terabytes of RAM and 10 petabytes of storage. High memory nodes have access to 256 GB of RAM, which is useful for large matrix multiplications. Odyssey uses SLURM²⁷ as its queue manager, allowing our analysis to utilize the large number of cores available.

Although the raw Bicep2 data set comprises roughly 3 TB of data, the data products from the steps in the matrix analysis chain use 17 TB of storage. Processing during the Bicep2 matrix analysis steps required roughly 1 million CPU hours. This represents a significant portion of the total computing demand of the entire Bicep2 analysis effort, which comprised roughly 6 million CPU hours.

We start by fixing the coefficients of the spherical harmonics of the temperature sky, a^T_lm. The a_lm for the constrained E-modes is found in Dvorkin et al. (2008) and can be derived following the steps of a simple Cholesky decomposition. The 2 × 2 covariance between T and E for each mode ($l,m$) is:

$$\begin{eqnarray}\left[\begin{array}{cc}{C}_{l}^{{TT}} & {C}_{l}^{{TE}}\\ {C}_{l}^{{TE}} & {C}_{l}^{{EE}}\end{array}\right].\end{eqnarray} \tag{ 75 }$$

The lower triangle Cholesky decomposition of this 2 × 2 matrix determines the variance and covariance in each mode:

$$\begin{eqnarray}\left[\begin{array}{cc}\sqrt{{C}_{l}^{{TT}}} & 0\\ {C}_{l}^{{TE}}/\sqrt{{C}_{l}^{{TT}}} & \sqrt{{C}_{l}^{{EE}}-{({C}_{l}^{{TE}})}^{2}/{C}_{l}^{{TT}}}\end{array}\right]\left[\begin{array}{c}{n}_{{lm}}^{T}\\ {n}_{{lm}}^{E}\end{array}\right]=\left[\begin{array}{c}{a}_{{lm}}^{T}\\ {a}_{{lm}}^{E},\end{array}\right]\end{eqnarray} \tag{ 76 }$$

where ${n}_{{lm}}^{T}$ and ${n}_{{lm}}^{E}$ are unit norm, complex random numbers. Since a^T_lm are known constraints, we can solve the system of equations:

$$\begin{eqnarray}&&\sqrt{{C}_{l}^{{TT}}}{n}_{{lm}}^{T}={a}_{{lm}}^{T}\\ &&{C}_{l}^{{TE}}/\sqrt{{C}_{l}^{{TT}}}{n}_{{lm}}^{T}+\sqrt{{C}_{l}^{{EE}}-{({C}_{l}^{{TE}})}^{2}/{C}_{l}^{{TT}}}{n}_{{lm}}^{E}={a}_{{lm}}^{E},\end{eqnarray} \tag{ 77 }$$

substituting the first equation into the second to arrive at an expression for a^E_lm in terms of a^T_lm and the power spectra, C_l. To ensure the maps are real, the condition ${a}_{{lm}}^{* }=(-{1}^{m}){a}_{l-m}$ is demanded.

Naively, we might expect the constrained a^E_lm to have lower variance in its power spectrum than the unconstrained a^E_lm since there is significant TE correlation and the temperature component is fixed. However, this is not necessarily the case. When the a^T_lm fluctuate high, the mean of the constrained ensemble of a^E_lm's is also high, resulting in increased variance in the ensemble of simulations. For particularly high a^T_lm, this can result in larger variance for the constrained simulations than the unconstrained simulations. This happens to be the case in the Bicep field near l = 150, as seen in Figure 23.

Zoom In Zoom Out Reset image size

In practice, we take the a^T_lm from the Planck temperature Needlet Internal Linear Combination (NILC) map (Planck Collaboration XII 2014). This map uses the multi-frequency coverage of Planck to remove the galactic contribution to microwave emission, leaving a high signal to noise map of the CMB temperature field. The map has some contamination near the galactic plane. However, we have found that the impact of this is very local and does not affect the higher galactic latitudes where the Bicep field is located. The noise level in the Planck temperature map is fractionally small compared to the temperature signal, and we have found that this noise contributes a similarly small fraction to the constrained realization of E-modes.

The temperature anisotropies in the Planck NILC map have been lensed by the intervening structure between us and the surface of last scattering. This means the a^T_lm calculated from the Planck NILC map contain the effects of lensing and when used in Equation (77), the lensing distortion propagates through to the constrained a^E_lm. Ideally, the a^T_lm in Equation (77) would be from the unlensed sky, however, in the absence of an accurate map of the lensing deflection field, we have no way of de-lensing the a^T_lm.

Because our power spectrum analysis is insensitive to the off-diagonal correlations among modes with $l\;=l^{\prime}$ that are produced by lensing, a reasonable workaround for this problem is to use the Planck NILC but also to use the lensed C^TT_l spectrum. The result is an ensemble of a^E_lm for the lensed a^T_lm, but which have the correct covariance given by C^TE_l. For the multipole range of interest in Bicep2, lensing has a small impact on a^T_lm. We use the unlensed power spectra for C^TE_l and C^EE_l.

Another subtlety incorporating lensing into constrained realizations is the question of how to simulate lensing of the polarization sky. This can be accomplished using LensPix(Lewis 2011) to numerically lens the primordial E-modes, which creates B-modes with the correct statistics. For the constrained realizations, this is impossible because the NILC map is lensed by the true sky lensing field. The true sky lensing field is not known with high signal to noise, and therefore we cannot lens the E-modes by the same field.

Our solution to this problem is to lens the E-modes with a random realization of the lensing field, but use the Planck 143 GHz map as the temperature field. This procedure ignores lensing correlations in TE and TB because the deflection field is different for the polarization and temperature. However, the mean of the lensing TE and TB is zero, and the additional variance in TE and TB caused by lensing is small enough that it can be ignored on degree scales.