Maximum Likelihood Foreground Cleaning for Cosmic Microwave Background Polarimeters in the Presence of Systematic Effects

# Maximum Likelihood Foreground Cleaning for Cosmic Microwave Background Polarimeters in the Presence of Systematic Effects

## Abstract

We extend a general maximum likelihood foreground estimation for cosmic microwave background polarization data to include estimation of instrumental systematic effects. We focus on two particular effects: frequency band measurement uncertainty, and instrumentally induced frequency dependent polarization rotation. We assess the bias induced on the estimation of the -mode polarization signal by these two systematic effects in the presence of instrumental noise and uncertainties in the polarization and spectral index of Galactic dust. Degeneracies between uncertainties in the band and polarization angle calibration measurements and in the dust spectral index and polarization increase the uncertainty in the extracted CMB -mode power, and may give rise to a biased estimate. We provide a quantitative assessment of the potential bias and increased uncertainty in an example experimental configuration. For example, we find that with 10% polarized dust, tensor to scalar ratio of , and the instrumental configuration of the EBEX balloon payload, the estimated CMB -mode power spectrum is recovered without bias when the frequency band measurement has 5% uncertainty or less, and the polarization angle calibration has an uncertainty of up to 4.

cosmic microwave background — instrumentation: polarimeters — methods: data analysis

## 1 Introduction

The inflationary paradigm posits that the universe underwent a period of exponential expansion within its first fraction of a second. One of the generic predictions of inflation is the production of stochastic gravitational wave background and the imprint of a -mode pattern in the polarization of the cosmic microwave background (CMB) radiation on large angular scales (Kamionkowski et al., 1997; Seljak & Zaldarriaga, 1997). The level of the inflationary gravitational wave -mode signal encodes information about the inflationary energy scale and is characterized by a tensor-to-scalar ratio  which quantifies the relative strength of gravity waves and density perturbations produced during inflation. The current upper limit from -mode observation alone is  (BICEP2/Keck and Planck Collaborations, 2015; Planck Collaboration, 2015a). On small angular scales there is more dominant source of -mode polarization: lensing of CMB photons by the large scale structure of the universe, which converts CMB -mode to -mode (Zaldarriaga & Seljak, 1998).

According to recent observations, polarized Galactic thermal dust emission is a significant source of contamination for the inflationary -mode signal (Gold et al., 2009; Fraisse et al., 2013; Planck Collaboration, 2014; BICEP2/Keck and Planck Collaborations, 2015). To monitor and subtract Galactic foreground contamination many CMB polarimeters observe at multiple frequency bands. In this paper we address foreground estimation for CMB polarimeters in the presence of two systematic effects: uncertainty in measurement of the frequency band, and frequency dependent polarization rotation.

Band measurement uncertainty is a common systematic effect. Absolute calibration uncertainty, which affects all signal components within a given frequency band by the same factor, has already been discussed in Stompor et al. (2009), Stivoli et al. (2010) and Errard & Stompor (2012). In this paper we focus on band measurement uncertainty, which affects sky signals with different emission spectra differently. In the context of foreground estimation with unknown foreground spectral parameters, the effect of band measurement uncertainties has not yet been studied much in detail. Church et al. (2003) gave an analytic calculation of the effect of error in the knowledge of bands given sky signals with fixed and known spectral shape. Bao et al. (2012) studied the impact of band measurement uncertainty in the presence of an achromatic half-wave plate (AHWP) and proposed an approximate dust estimation and subtraction process.

For CMB polarimeters using an AHWP, such as the E and B experiment (EBEX) (Reichborn-Kjennerud et al., 2010), or using a sinuous antenna multi-chroic pixel (Arnold et al., 2010), these elements induce a frequency-dependent instrumental polarization rotation. Uncertainties in foreground emission laws and in the relative intensity of the CMB and foregrounds - at high frequencies predominantly Galactic dust - induce uncertainties in correcting the instrumental polarization rotation.

The goal of this paper is to develop a foreground estimation formalism in the presence of these two systematic effects and to assess its performance in a practical instrumental configuration. Our work is based on the maximum likelihood component separation method presented in Stompor et al. (2009). For concreteness we adopt the AHWP model, frequency bands, and approximate noise information that were applicable to the design of the EBEX balloon-borne polarimeter.

In Section 2 we explain the theoretical framework of the maximum likelihood foreground cleaning algorithm. In Section 3 we briefly describe the details of our simulation. The results of the simulations are reported in Section 4. Finally in Section 5 we make the concluding remarks.

## 2 Theoretical Framework

In this section we discuss the mathematical framework of the maximum likelihood foreground estimation. This framework is developed entirely in the map domain. We start with a brief review of the basic formalism in the absence of systematic effects. We then extend it to include band measurement uncertainty and frequency dependent polarization rotation.

### 2.1 Basic Formalism

We model the sky signal observed in multiple frequency bands for a single pixel in the map as

 dp = Ap sp + np. (1)

Here the subscript denotes quantities for a single pixel;  is the data vector containing the measured signals for frequency bands and Stokes parameters;  is the underlying sky signal vector for sky signal components and Stokes parameters;  is the noise vector for frequency bands and Stokes parameters;  is the component ‘mixing matrix’, which maps the sky signals to the observed signals. At each pixel p, the mixing matrix  has blocks. Each block is an by diagonal matrix with all of its diagonal elements equal to each other. It is parameterized by a set of unknown parameters describing the spectral shape of the components. In this paper we assume that the parameters are spatially uniform across the patch. For pixels we remove the subscript and the data model becomes

 d = A s + n. (2)

The likelihood function for the data set is

 −2lnL(s,β) = const + (d−A s)t N−1 (d−A s), (3)

where  is the noise covariance matrix. When there is no correlated noise between different pixels,  is a rank square, symmetric, block diagonal noise matrix and the full data likelihood can be calculated as the sum of likelihood values calculated from each pixel. The likelihood is maximized when

 s = (At N−1 A)−1 At N−1 d, (4)

and by substituting Eq. 4 into Eq. 3 we find

 −2lnL = const − (At N−1 d)t (At N−1 A)−1 (At N−1 d). (5)

For a given set of parameters , the mixing matrix  can be calculated straightforwardly. This suggests a two step foreground estimation process: first finding the set of parameters that maximizes the likelihood function (Eq. 5), then calculating the component signals given the maximum likelihood parameters using Eq. 4.

### 2.2 Extension of the Basic Formalism

When extending the basic formalism to include systematic effects, the mixing matrix  takes a more complicated form while the main steps of the foreground estimation remain the same.

#### Frequency Band Measurement Uncertainty

A frequency band is expressed in terms of transmission as a function of frequency . We characterize an arbitrary band in terms of two parameters: the band-center and band-width. We use the same definitions as in Runyan et al. (2003). The band-center is

 νc = ∫ν T(ν) dν∫T(ν) dν, (6)

which is the mean frequency weighted by the transmission. The band-width is

 Δν = νU−νL, (7)

where the lower and upper band edges are defined as

 νU = ∫∞νcνT′(ν)dν∫∞νcT′(ν)dν , νL = ∫νc0νT′(ν)dν∫νc0T′(ν)dν, (8)

and is the derivative of the transmission curve . This definition of band edges favors frequencies where sharp transitions happen in the transmission curve. For example, for a top-hat band the lower and upper band edges are where the transmission turns on and off.

There are uncertainties in the measurements of frequency bands. To the extent that bands in this paper are characterized with two parameters, band center and width, the measured band uncertainties convert to uncertainties in these parameters. For any incoming sky signal, band measurement uncertainties propagate to uncertainties in the total measured in-band power, which also depends on the spectral shape and amplitude of the sky signal. Each sky signal component is affected differently. We parameterize the measured power uncertainties in terms of a vector of scaling coefficients  which contains components. Here we define the nominal band as the ‘true’ underlying band of the instrument and the assumed band as the band determined from the measurements. A scaling coefficient is the ratio of the band integrated power in the nominal band around frequency and the assumed band for signal component . Uncertainties in the band measurements translate to a range of probable . Because the scaling coefficient also depends on the emission spectral shape, each sky signal component has a different range of probable values. Note that does not depend on the overall amplitude of the sky components. In the case where the band is measured accurately, the scaling coefficient is unity for any sky signal. Figure 1 illustrates the scaling coefficients; it shows  for CMB and polarized Galactic thermal dust in two frequency bands centered at 150 GHz and 250 GHz. Because in this paper we assume that the spectral shape of the sky components are spatially uniform across the patch, the scaling coefficients  are also spatially uniform.

In our extended foreground estimation formalism we find the maximum likelihood solution for the foreground spectral parameters  and for the scaling coefficients . A range of probable scaling coefficients, as obtained from the likelihood function, are mappable to possible combinations of band-centers and band-widths.

The technical implementation is carried out by introducing a new mixing matrix , which has the same rank as the original mixing matrix , but explicitly includes the scaling coefficients. Each by block in  is multiplied by the scaling coefficient for the corresponding frequency channel and sky component. With the data model is

 d=AB(β,η) s + n (9)

and the likelihood function is

 −2lnLB(s,β,η) = const + (d−AB s)t N−1 (d−AB s). (10)

When the likelihood reaches maximum we have

 −2lnLB = const − (At% B N−1 d)t (AtB N−1 AB)−1 (AtB N−1 d), (11)

and

 s = (AtB N−1 AB)−1 AtB N−1 d. (12)

There exist two degeneracies in the problem. First, the signal level is degenerate with the scaling coefficients. For a particular sky component, simultaneously multiplying all the scaling coefficients by a factor is equivalent to multiplying the sky signal by the same factor. Second, the spectral parameters and the scaling coefficients are degenerate. Tilting the spectrum is equivalent to changing the scaling coefficients by corresponding factors. Figure 2 shows an illustration of this degeneracy. In the example, a mis-estimate of the polarized Galactic thermal dust spectral index by 0.2 has the same and as an 8.5 GHz shift of the band-center at the 150 band toward lower frequency and a 6 GHz reduction of the band-width at the 250 band.

The degeneracies can be broken by setting priors on the scaling coefficients based on the uncertainty of the frequency band measurements. If we assume is measured with an uncertainty around , the likelihood function becomes

 −2lnLB(β,η) = const − (AtBN−1d)t (AtBN−1AB)−1 (AtBN−1d) + ∑i(ηi−¯ηi)2σ2ηi, (13)

where the last term represents the priors on the scaling coefficients. We use Eq. 13 to find the best fit spectral parameters and the scaling coefficients. We then recover the sky component signals using the best fit parameters.

#### Frequency Dependent Polarization Rotation

Frequency dependent polarization rotation mixes the and signals. The amount of rotation within a band depends on the spectral shape of the incoming signals, the characteristic parameters of the instrument, and the frequency band shape.

To incorporate instrumental frequency dependent rotation into foreground estimation, we introduce a new mixing matrix  which has the same rank as the original mixing matrix . The mixing matrix is parameterized by the spectral parameters  and a vector  of rotation angles for each sky signal component in each frequency band. Compared to , the - blocks for each sky signal component at each frequency band in  are multiplied by rotation matrices with non-zero rotation angles . Thus each - block has non-zero off-diagonal elements and  is not block diagonal. A requirement on the uncertainty of the rotation angles translates to a requirement on the uncertainty of polarization angle calibration.

With  the data model is

 d = AR(β,θ) s + n, (14)

and the likelihood function is

 −2ln LR(s,β,θ) = const + (d − AR s)t N−1 (d − AR s). (15)

When the likelihood reaches its maximum we have

 −2ln LR = const − (AtR N−1 d)t (AtR N−1 AR)−1 (AtR N−1 d), (16)

and

 s = (AtR N−1 AR)−1 AtR N−1 d. (17)

There exists a degeneracy between the frequency dependent polarization rotation and the polarization angle of the incoming signal. Changing the band averaged rotation angles of a sky component in all bands by a given amount is equivalent to changing the incoming polarization signal by the same amount. This degeneracy can be broken by setting priors on the rotation angles based on the uncertainty of the polarization angle calibration. If we assume the band averaged rotation angle is determined within an uncertainty around given the knowledge of the instrument and the incoming signals, the likelihood function becomes

 −2lnLR(β,θ) = const − (AtRN−1d)t (AtRN−1AR)−1 (AtRN−1d) + ∑i(θi−¯θi)2σ2θi, (18)

where the last term originates from the priors on the rotation angles.

#### Combining Band Measurement Uncertainty and Frequency Dependent Polarization Rotation

To incorporate both band measurement uncertainty and frequency dependent polarization rotation at the same time we define a new mixing matrix ,,, which has the same rank as the original mixing matrix . Compared to , the mixing matrix is parameterized by the spectral parameters , the scaling coefficients  and the band averaged rotation angles . Each by block in  is multiplied by the in-band scaling coefficient  for the corresponding signal component and frequency band. The - block within each by block is multiplied by a rotation matrix with rotation angle .

With the definition of the data model is

 d = AC(β,η,θ) s + n, (19)

and the likelihood function is

 −2lnLC(s,β,η,θ) = const + (d − AC s)t N−1 (d − AC s). (20)

When the likelihood reaches its maximum we have

 −2lnLC = const − (At% C N−1 d)t (AtC N−1 AC)−1 (AtC N−1 d), (21)

and

 s = (AtC N−1 AC)−1 AtC N−1 d. (22)

Due to the degeneracies mentioned in Sects. 2.2.1 and 2.2.2, to recover the signals accurately we set priors to the scaling coefficients  and the band averaged rotation angles  based on the uncertainty of the measurements. The likelihood function becomes

 −2lnLC(β,η,θ) = const − (AtC N−1 d)t (AtC N−1 AC)−1 (AtC N−1 d)+ ∑i(ηi−¯ηi)2σ2ηi + ∑i(θi−¯θi)2σ2θi, (23)

where the last two terms come from the priors.

### 2.3 Error Propagation

Uncertainties in the spectral parameters and instrumental parameters in the mixing matrix propagate to the estimated sky component signals. We estimate the covariance matrix for the parameters using the inverse of the curvature matrix calculated at the maximum likelihood. In the basic case without systematic effects we have

 \boldmath~N\unboldmathββ = [(A,βis)t N−1 (A,βjs)−(A,βiβjs)t N−1 (d−A s)−Mt^NM]−1 (24)

for the spectral parameters , where  is the estimated uncertainty matrix of the sky component maps assuming the mixing matrix is perfectly known

 ^N = (At N−1 A)−1, (25)

and the  matrix is defined as

 M = At N−1 A,βs−At,β N−1 (d−A s). (26)

and  are the first and second order derivatives of the mixing matrix with respect to the spectral parameters, estimated at the best fit values of the spectral parameters. The noise matrix of the estimated sky signals is

 ~Nss=^N+(^NM) ~Nββ (^NM)t. (27)

When systematic effects are considered, the mixing matrix  in Eq. 24, Eq. 26 and Eq. 27, is changed to the corresponding extended mixing matrix (,  or ) and the derivative of the mixing matrix is performed with respect to all corresponding unknown parameters including the spectral parameters , the scaling coefficients  and the band averaged rotation angles .

## 3 Simulations

We test the extended foreground estimation formalism using simulations. We select a sky patch centered on RA = , DEC = and include only the CMB and Galactic dust as sky signals. The input CMB angular power spectra are generated with Code for Anisotropies in the Microwave Background (CAMB) described in Lewis et al. (2000) using the Wilkinson Microwave Anisotropy Probe (WMAP) seven-year best-fit cosmological parameters (Komatsu et al., 2011) and a tensor to scalar ratio of . With a flat-sky approximation (Kaiser, 1992) we generate random Fourier amplitudes in temperature and polarization to match the power spectra, which are then transformed to CMB maps on the chosen square patch of sky. We adopt the process to generate Galactic dust foreground detailed in Stivoli et al. (2010). The dust intensity and its frequency scaling are given by an  K blackbody with a 1.59 power-law emissivity according to the recent Planck measurement (Planck Collaboration, 2015b). The dust polarization fraction is set to 10%. Both the dust frequency scaling and the dust polarization fraction are assumed to be spatially uniform. The large scale polarization angle patterns () are derived from WMAP dust polarization template (Page et al., 2007). On small angular scales a Gaussian fluctuation power is added using the recipe described in Giardino et al. (2002). Figure 3 shows the and maps of the dust signal and one realization of CMB in the selected sky patch.

The observations are simulated in three top-hat bands centered at 150 GHz, 250 GHz and 410 GHz. The frequency ranges for the three bands and their noise levels per square pixel in the and maps are listed in Table 1. The noise is assumed to be uniform and white across the patch. The noise realization is added to the signal in the map domain in all simulations.

In the simulation we include a continuously rotating achromatic half-wave plate (AHWP) which induces a frequency dependent polarization rotation. The AHWP is composed of a stack of five single sapphire half-wave plates with an average thickness of 1.65 mm. The relative orientation angles between the optical axis of the plates and the first plate in the stack are 0, 28, 94, 28, 0. The ordinary and extraordinary refractive indices of sapphire at cryogenic temperature are used to calculate the instrumental polarization rotation angles, i.e., and  (Afsar & Chi, 1991; Loewenstein et al., 1973). We model the AHWP using Mueller matrices formalism described in Matsumura et al. (2009). Given the frequency bands, the input dust spectrum and the AHWP parameters, the rotation angles are 115.66, 103.64, 110.76 for CMB and 115.02, 103.68, 112.88 for dust at 150, 250, 410 GHz bands, respectively . Figure 4 shows the observed and maps at 150 GHz band. They include both CMB, and dust, rotation induced by the AHWP, and instrumental noise. The parameters of the AHWP and the noise levels were chosen to be close to those designed for the EBEX balloon-borne instrument.

The and power spectra are calculated simultaneously using a flat-sky approximation (Kaiser, 1992). For each set of input parameters, the simulation is run 100 times with different CMB and noise realizations while the dust signal is kept the same. The result shown for any given bin is the mean of the 100 simulations and the error bar is the standard deviation. Specifically, the error bars of the input CMB spectrum are the cosmic variances while the error bars of the output CMB spectrum also include uncertainties from the instrumental noise and the propagation of the parameter uncertainties. We also calculate residual power spectra. For each of the 100 simulations the residual spectrum is the difference between the estimated CMB spectrum, after estimation of the foreground, and the input spectrum. We plot the mean and the standard deviations of the 100 simulations. We consider the estimated CMB -mode power spectrum ‘biased’ using two criteria such that each provides somewhat different information: (1) when the estimated CMB -mode signal power is more than twice the cosmic variance away from the input in any bin, and (2) when the residual power spectrum is more than one sigma away from zero. The first criterion quantifies bias relative to irreducible cosmic variance and is particularly useful for noiseless simulations. The second is more relevant when instrument noise is included and is less stringent than the first.

For the simulations the data model in the basic formalism

 d = A s + n, (28)

becomes

 \scalebox0.78⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝I150Q150U150I250Q250U250I410Q410U410⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝A150,c00A150,d(β)000A150,c00A150,d(β)000A150,c00A150,d(β)A250,c00A250,d(β)000A250,c00A250,d(β)000A250,c00A250,d(β)A410,c00A410,d(β)000A410,c00A410,d(β)000A410,c00A410,d(β)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ICMBQCMBUCMBIdustQdustUdust⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠+⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝nI150nQ150nU150nI250nQ250nU250nI410nQ410nU410⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (29)

Here , and are the polarization stokes parameters. The subscripts denote frequency bands and components. The elements in the mixing matrix related to CMB are known since we know the CMB spectrum. For any given frequency , the element in the mixing matrix related to dust is

 Aν,d(β) = (νν0)βBν(Td,ν), (30)

where is a reference frequency typically set to the center of the highest frequency channel of the experiment, is the blackbody spectrum, is the dust temperature (which is assumed to be known) and is the spectral index of dust which is the sole unknown parameter in this case.

When band measurement uncertainties are included, the mixing matrix becomes

 Unknown environment '% (31)

There are seven unknown parameters in : the dust spectral index and six scaling coefficients.

When the frequency dependent polarization rotation is included, the mixing matrix becomes

 AR=⎛⎜ ⎜⎝AR150,cAR150,dAR250,cAR250,dAR410,cAR410,d⎞⎟ ⎟⎠, (32)

with the blocks  for CMB and dust at frequency band being

 ARν,c=⎛⎜ ⎜⎝Aν,c000Aν,ccosθν,c−Aν,csinθν,c0Aν,csinθν,cAν,ccosθν,c⎞⎟ ⎟⎠,
 ARν,d=⎛⎜ ⎜⎝Aν,d(β)000Aν,d(β)cosθν,d(β)−Aν,d(β)sinθν,d(β)0Aν,d(β)sinθν,d(β)Aν,d(β)cosθν,d(β)⎞⎟ ⎟⎠,

where the subscript runs over 150, 250, and 410 GHz. There are seven unknown parameters in : the dust spectral index and six rotation angles.

When both the band measurement uncertainty and the frequency dependent rotation effect of the AHWP are considered, the mixing matrix becomes

 AC=⎛⎜ ⎜⎝AC150,cAC150,dAC250,cAC250,dAC410,cAC410,d⎞⎟ ⎟⎠, (33)

with the blocks  for CMB and dust at frequency channel being

 ACν,c=⎛⎜ ⎜⎝ην,cAν,c000ην,cAν,ccosθν,c−ην,cAν,csinθν,c0ην,cAν,csinθν,cην,cAν,ccosθν,c⎞⎟ ⎟⎠,
 ACν,d=⎛⎜ ⎜⎝ην,dAν,d(β)000ην,dAν,d(β)cosθν,d% (β)−ην,dAν,d(β)sinθν,d(β)0ην,dAν,d(β)sinθν,d% (β)ην,dAν,d(β)cosθν,d(β)⎞⎟ ⎟⎠,

where the subscript runs over 150, 250 and 410 GHz. There are a total of 13 unknown parameters in : the dust spectral index, six scaling coefficients and six band averaged rotation angles. Note that reduce to in the case of and reduce to in the case of . In our simulations, when setting priors to parameters, we assume a Gaussian prior with quoted FWHM values.

## 4 Results

### 4.1 Consistency Check

In the absence of band measurement and polarization angle calibration uncertainties the extended formalism reproduces the input spectra with no bias; this is shown in Fig. 5. For this simulation we included the systematic effects but they were assumed to be perfectly known. In all the power spectra plots shown in the rest of this section, the data points of the estimated CMB -mode power spectrum are plotted offset along the x-axis for clarity. We also plot the theoretical underlying CMB -mode power spectrum with .

Next we assess the intrinsic bias in the estimated CMB signal induced by the formalism in the absence of instrumental noise. Due to the degeneracies discussed in Sec. 2.2, here we consider the systematic effects in only one band at a time while assuming the other two bands are perfectly known. All the parameters in the mixing matrix are estimated without any priors. We test the formalism including band measurement uncertainty only, frequency dependent polarization rotation only, or both systematic effects. The estimated CMB -mode signal is not biased in all of the cases. Figure 6 shows an example when both band measurement uncertainty and frequency dependent polarization rotation are considered at 150 GHz band. When the 150 GHz band-center or band-width is mis-estimated by 15 GHz, the estimated scaling coefficients and CMB -mode signal are without any bias. Similarly, the formalism does not bias the CMB -mode signal with a mis-estimation of 15 GHz of the band-center or band-width at 250 or 410 GHz band. We also see similar results when only one of the systematic effects is considered.

### 4.2 Band Measurement Uncertainty

When estimating the CMB -mode signal in the presence of band measurement uncertainty and instrumental noise, all the band averaged rotation angles are assumed to be perfectly known.

First we address the cases where there are uncertainties in only one of the bands. As discussed in Sec. 2.2.1 we observe a degeneracy between the dust spectral index and the dust scaling coefficient. When noise is present in the observation, changing the tilt of the dust spectrum or the scaling of the dust signal in one band give similar likelihood. The degeneracy causes bias in the recovered CMB power spectrum, particularly in the low- bins where the dust signal is high and the inflationary -mode signal resides. Figure 7 shows the 2-D likelihood plot between the dust spectral index and the dust scaling coefficient and the estimated CMB -mode power spectrum when there are uncertainties at 150 GHz band only. The points along the bright diagonal feature in the likelihood plot have similar likelihood. When the parameters are mis-estimated, the estimated -mode signal is biased by at low .

Setting priors on the scaling coefficients according to the band measurement uncertainty limits the parameter space in which the best fit values are searched for. In this study the priors are centered on the input value which is one. In the case where the uncertainty of only the 150 GHz band is considered, with 15% Gaussian priors on the scaling coefficients for both CMB and dust the estimated CMB -mode signal is not biased, as shown in Fig. 8.

When uncertainties of all three bands are included the degeneracy between the scaling coefficients and the sky signal level requires priors on all scaling coefficients. Figure 9 shows two cases where 5% or 10% Gaussian priors are set on all scaling coefficients. When all scaling coefficients have 5% Gaussian priors the -mode signal is estimated accurately. When the priors are relaxed to 10%, the estimated signal in the lowest bin is 2.5 times the cosmic variance away from the input value. The bias comes from the degeneracy between the dust scaling coefficients and the dust spectral index and therefore the bias is strongest at the range in which dust is most intense. If we factor in the effect of the instrumental noise the residual power spectrum is consistent with zero in all bins.

Specific priors on the scaling coefficients is analogous to known measurement uncertainties. Table 2 gives the mapping between prior interval on the scaling coefficients and uncertainties in band-center and band-width for the three top-hat bands used in the simulations. The values for band-center uncertainty are calculated assuming the band-width has no uncertainty and vice versa.

In practice, however, both band-center and band-width have uncertainty simultaneously. As an example, in Fig. 10 we show the contour of the scaling coefficient deviation in the 2-D parameter space of band-center and band-width for Galactic dust at 250 GHz band. Similar Figures are obtained for the CMB and other frequency bands as well. Using such Figures, which are easily calculable given any band parameters, one can map band measurement uncertainty to scaling coefficient priors.

### 4.3 Frequency Dependent Polarization Rotation

When assessing only the frequency dependent polarization rotation in the presence of instrumental noise, all the scaling coefficients are assumed to be perfectly known. We first consider the uncertainty of the band averaged rotation angles in only one band while assuming the rotation angles in the other two bands are perfectly known. Since there is no degeneracy between the incoming polarization angle and the AHWP rotation angles in this case, we can estimate the CMB -mode signal accurately without priors on the rotation angles. Figure 11 shows the result where the uncertainty of the polarization angle calibration in only 150 GHz band is considered. We get similar results for 250 GHz and 410 GHz band.

When uncertainties of polarization angle calibration in all three bands are included, the degeneracy between the frequency dependent polarization rotation and the polarization angle of the incoming signal requires priors on the band-averaged rotation angles. Here we center the priors on the input value. Figure 12 shows that with 4 Gaussian priors on all band-averaged rotation angles, the final estimated CMB -mode signal is not biased. When the priors are relaxed to 10 the estimated CMB -mode signal is more than twice the cosmic variance away from the input in a few bins at . Less stringent priors on the band averaged rotation angles result in a bigger mis-estimate of the polarization angle of the incoming signal. This mis-estimate induces leakage from CMB -mode signal to -mode signal which has a bigger effect at high . When other noise sources are included the residual power spectrum is consistent with zero in both cases.

### 4.4 Combining Band Measurement Uncertainty and Frequency Dependent Polarization Rotation

When both systematic effects are included in the presence of instrumental noise, the formalism estimates the scaling coefficients and the band averaged polarization rotation angles simultaneously. We first assess the bias in the estimated CMB -mode signal including the systematic effects at only one of the bands while assuming the other two bands are perfectly known. In Fig. 13 we show the results where there are systematic effects at 150 GHz only. The estimated CMB -mode signal is not biased with 15% Gaussian priors on the scaling coefficients for CMB and dust and no priors on band averaged rotation angles. We get similar results when there are systematic effects at 250 or 410 GHz band only.

When the systematic effects are considered for all three bands, priors are needed for all parameters due to the degeneracies. Figure 14 shows that with 5% Gaussian prior on scaling coefficients and 4 Gaussian priors on band averaged rotation angles we can estimate CMB -mode signal without any bias.

### 4.5 Foreground Estimation with Measured EBEX Bands

As a practical example, we estimate foreground using the measured EBEX bands and their uncertainties, which are shown in Fig. 15 (Zilic, 2014).

We carry out simulations of the band measurement. Each realization is a simulated measurement of the band given the measured frequency bin errors. For each realization we compute the scaling coefficients compared to the measured band, which are taken to be the nominal values of the measurement, and band averaged rotation angles for CMB and Galactic dust. The standard deviation from the simulated results, listed in Table 3, are considered as the priors on the corresponding parameters.

Figure 16 shows the power spectra of the simulation with the priors listed in Table 3. We find that with an expreimtn like EBEX can estimate the B-mode spectrum without bias at a level comparible to either cosmic variance or instrumental noise.

## 5 Summary

Recent observations suggest that Galactic dust is a significant contaminating source in all regions of the sky for CMB polarimeters targeting the inflationary -mode signal. The Galactic dust signal needs to be removed and robust foreground estimation is essential. We presented a general maximum likelihood foreground estimation formalism in the presence of two systematic effects: band measurement uncertainty and frequency dependent polarization rotation. The formalism fits the systematic effects simultaneously with the dust spectral index which allows for imperfect knowledge of the instrumental parameters.

We found several degeneracies: 1) between the dust spectral index and the dust scaling coefficients, 2) between the scaling coefficients and the signal levels, and 3) between the frequency dependent polarization rotation and the signal polarization angles. Due to the degeneracies, priors on the parameters are needed in order estimate the CMB -mode signal accurately. We quantified the degeneracies and showed, as an example, that a sub-orbital experiment like EBEX should not be limited in estimating and subtracting Galactic dust if it has 10% polarization fraction, if the tensor to scalar ratio , and with band measurement and polarization angle calibration uncertainty of 5% and 4, respectively. Such an experiment may not be limited with even higher polarization fractions or lower values. A detailed calculation is necessary for other concrete cases.

Band measurement uncertainty is a common systematic effect for all CMB instruments. Frequency dependent polarization rotation by elements in the focal plane is not unique to the use of an achromatic half-wave plate either. Sinuous antenna detectors also have frequency dependent polarization rotation (O’Brient et al., 2008). Therefore, the foreground estimation formalism developed in this paper has general applicability.

This work received major support through an NSF grant NSF ANT-094513. We are thankful for the computing resources provided by the Minnesota Supercomputing Institute. EBEX is supported through NASA grants NNX08AG40G and NNX07AP36H. C.Bao acknowledges the support from the Doctoral Dissertation Fellowship from the University of Minnesota. C.Baccigalupi acknowledges the support from the INDARK INFN Grant.

### Footnotes

1. affiliation: University of Minnesota School of Physics and Astronomy, Minneapolis, MN 55455
2. affiliation: SISSA, Astrophysics Sector, via Bonomea 265, Trieste 34136, Italy
3. affiliation: Hamline University, Saint Paul, MN 55104
4. affiliation: University of Minnesota School of Physics and Astronomy, Minneapolis, MN 55455
5. affiliation: Imperial College, London, SW72AZ, England, United Kingdom
6. affiliation: Laboratoire Astroparticule et Cosmologie (APC), 75205 Paris Cedex 13, France

### References

1. Afsar, M. N., & Chi, H. 1991, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 16th International Conference on Infrared and Millimeter Waves, iSBN 0-8194-0707-0, M. R. Siegrist, M. Q. Tran and T. M. Tran, Lausanne, Switzerland, Aug 26-30, 1991
2. Arnold, K., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7741, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
3. Bao, C., et al. 2012, ApJ, 747, 97
4. BICEP2/Keck and Planck Collaborations. 2015, Physical Review Letters, 114, 101301
5. Church, S., Knox, L., & White, M. 2003, ApJ, 582, L63
6. Errard, J., & Stompor, R. 2012, Phys. Rev. D, 85, 083006
7. Fraisse, A. A., et al. 2013, JCAP, 4, 47
8. Giardino, G., Banday, A. J., Górski, K. M., Bennett, K., Jonas, J. L., & Tauber, J. 2002, A&A, 387, 82
9. Gold, B., et al. 2009, ApJS, 180, 265
10. Kaiser, N. 1992, ApJ, 388, 272
11. Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. Lett., 78, 2058, astro-ph/9609132
12. Komatsu, E., et al. 2011, ApJS, 192, 18
13. Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
14. Loewenstein, E. V., Smith, D. R., & Morgan, R. L. 1973, Appl. Opt., 12, 398
15. Matsumura, T., Hanany, S., Ade, P., Johnson, B. R., Jones, T. J., Jonnalagadda, P., & Savini, G. 2009, Appl. Opt., 48, 3614
16. O’Brient, R., et al. 2008, in Presented at the Society of Photo-Optical Instrumentation Engineers (SPIE) Conference, Vol. 7020, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
17. Page, L., et al. 2007, ApJS, 170, 335
18. Planck Collaboration. 2014, ArXiv e-prints, 1409.5738
19. —. 2015a, ArXiv e-prints, 1502.01589
20. —. 2015b, A&A, 576, A107
21. Reichborn-Kjennerud, B., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7741, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
22. Runyan, M. C., et al. 2003, ApJS, 149, 265
23. Seljak, U., & Zaldarriaga, M. 1997, Phys. Rev. Lett., 78, 2054, astro-ph/9609169
24. Stivoli, F., Grain, J., Leach, S. M., Tristram, M., Baccigalupi, C., & Stompor, R. 2010, MNRAS, 408, 2319
25. Stompor, R., Leach, S., Stivoli, F., & Baccigalupi, C. 2009, MNRAS, 392, 216
26. Zaldarriaga, M., & Seljak, U. 1998, Phys. Rev. D, 58, 023003
27. Zilic, K. T. 2014, PhD thesis, University of Minnesota
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters