A PERTURBATIVE ANALYSIS OF SYNCHROTRON SPECTRAL INDEX VARIATION OVER THE MICROWAVE SKY

From a mixture of n physical components present in different frequency bands, our aim is to isolate a best-fit map of each one of them by removing the rest. This can be achieved by reconstructing any one component first and then repeating the method for all other components. Let us denote the component to be reconstructed by a Greek letter subscript, e.g., α and other components to be removed by any Roman letter subscript taking values from the difference set $\{1,2,\ldots ,n\}\mbox{--}\{\alpha \}$. We denote a shape vector for the component to be reconstructed by ${[{\boldsymbol{s}}]}_{\alpha }$, while the ith component from remaining $n-1$ components has a shape vector of ${[{\boldsymbol{s}}]}_{i}$. To reconstruct the cleaned template, ${\hat{T}}_{\alpha }(p)$, for component α, we form a linear combination of available maps, following

$$\begin{eqnarray}&&{\hat{T}}_{\alpha }(p)=\displaystyle \sum _{i=1}^{{n}_{b}}{w}_{\alpha }^{i}{Q}_{i}(p),\end{eqnarray} \tag{ 9 }$$

where the linear weights, ${w}_{\alpha }^{i}$, for α component is obtained by minimizing the total pixel–pixel variance, ${\sigma }_{\alpha }^{2}$, of the cleaned template, ${\hat{T}}_{\alpha }(p)$. Pixel variance of this template is defined as

$$\begin{eqnarray}&&{\sigma }_{\alpha }^{2}=\displaystyle \sum _{p=1}^{{{ \mathcal N }}_{\mathrm{pix}}}{\left({\hat{T}}_{\alpha }(p)-\displaystyle \sum _{p^{\prime} }\displaystyle \frac{{\hat{T}}_{\alpha }(p^{\prime} )}{{{ \mathcal N }}_{\mathrm{pix}}})\right)}^{2},\end{eqnarray} \tag{ 10 }$$

where ${{ \mathcal N }}_{\mathrm{pix}}$ denotes the total number of pixels available for analysis in a single frequency map after any suitable mask has been applied. Using the matrix notation above, the equation can easily be written as

$$\begin{eqnarray}&&{\sigma }_{\alpha }^{2}={[{\boldsymbol{W}}]}_{\alpha }[{\boldsymbol{V}}][{\boldsymbol{W}}{]}_{\alpha }^{T},\end{eqnarray} \tag{ 11 }$$

where, ${[{\boldsymbol{W}}]}_{\alpha }$ is a $1\times {n}_{b}$ weight matrix, $({w}_{\alpha }^{1},{w}_{\alpha }^{2},\ldots .,{w}_{\alpha }^{{n}_{b}})$, for component α. $[{\boldsymbol{V}}]$ represents the ${n}_{b}\times {n}_{b}$ covariance matrix, with elements ${V}_{{ii}^{\prime} }$ representing the covariance between frequency maps Q_i(p) and ${Q}_{i^{\prime} }(p^{\prime} )$. As described in detail in Saha et al. (2008), the choice of weights that minimizes ${\sigma }_{\alpha }^{2}$ is given by

$$\begin{eqnarray}&&{[{\boldsymbol{W}}]}_{\alpha }=\displaystyle \frac{{[{\boldsymbol{s}}]}_{\alpha }{[{\boldsymbol{V}}]}^{+}}{{[{\boldsymbol{s}}]}_{\alpha }{[{\boldsymbol{V}}]}^{+}[{\boldsymbol{s}}{]}_{\alpha }^{T}},\end{eqnarray} \tag{ 12 }$$

where ${[{\boldsymbol{V}}]}^{+}$ represents the Moore–Penrose generalized inverse (Moore 1920; Penrose 1955) of $[{\boldsymbol{V}}]$. We note that, in practice, ${[{\boldsymbol{s}}]}_{\alpha }$ is not known a priori since the spectral indices for the components to be cleaned (e.g., synchrotron, thermal dust etc.) are also, at the same time, to be determined from the data. However, as described in Section 5.2, we find these weights for different components by varying their spectral indices in different pre-defined intervals. Each set of weights thus obtained then can be used to obtain a cleaned template for a given component.

Below, we derive an analytical expression for weights (i.e., Equation (15)) that minimizes variance due to physical components by giving an outline of the detailed mathematical procedure. The basic idea behind this equation is that under the assumption of negligible detector noise⁶ , reconstruction of any desired component (with known spectral index) can be achieved if the weight vector becomes "parallel" to the shape vector of this component, while simultaneously being orthogonal to the shape vector of each one of the rest of the components that are to be removed.

Using variance and covariance of different components, we can write

$$\begin{eqnarray}&&{[{\boldsymbol{V}}]={\sigma }_{\alpha }^{2}{[{\boldsymbol{s}}{]}_{\alpha }^{T}{[{\boldsymbol{s}}]}_{\alpha }+[{\boldsymbol{r}}{]}_{\alpha }^{T}[{\boldsymbol{s}}]}_{\alpha }+[{\boldsymbol{s}}{]}_{\alpha }^{T}[{\boldsymbol{r}}]}_{\alpha }+[{\boldsymbol{F}}],\end{eqnarray} \tag{ 13 }$$

where ${[{\boldsymbol{r}}]}_{\alpha }$ is a $1\times {n}_{b}$ matrix describing random chance correlation between the component to be reconstructed and all of the other components at each of the n_b frequency bands. $[{\boldsymbol{F}}]$ describes the ${n}_{b}\times {n}_{b}$ covariance matrix of n–1 components, hence it has a rank of n–1. The above equation can be written as

$$\begin{eqnarray}&&{[{\boldsymbol{V}}]={\sigma }_{\alpha }^{2}{[{\boldsymbol{s}}{]}_{\alpha }^{T}[{\boldsymbol{s}}]}_{\alpha }+[{\boldsymbol{A}}]}_{1},\end{eqnarray} \tag{ 14 }$$

where ${[{\boldsymbol{A}}]}_{1}={[{\boldsymbol{r}}{]}_{\alpha }^{T}{[{\boldsymbol{s}}]}_{\alpha }+[{\boldsymbol{A}}]}_{2}$ and ${[{\boldsymbol{A}}]}_{2}={[{\boldsymbol{s}}{]}_{\alpha }^{T}[{\boldsymbol{r}}]}_{\alpha }+[{\boldsymbol{F}}]$. Using the above equation, we can apply the generalized Shermann–Morisson formula (Meyer 1973; Baksalary et al. 2003) successively to obtain both the numerator and the denominator of Equation (12) (e.g., as in Saha et al. 2008). After a long set of algebra and neglecting all terms containing ${[{\boldsymbol{r}}]}_{\alpha }$, Equation (12) reduces to

$$\begin{eqnarray}&&{[{\boldsymbol{W}}]}_{\alpha }=\displaystyle \frac{{[{\boldsymbol{s}}]}_{\alpha }([{\boldsymbol{I}}]-[{\boldsymbol{F}}]{[{\boldsymbol{F}}]}^{+})}{{[{\boldsymbol{s}}]}_{\alpha }([{\boldsymbol{I}}]-[{\boldsymbol{F}}]{[{\boldsymbol{F}}]}^{+}){[{\boldsymbol{s}}]}_{\alpha }^{T}}.\end{eqnarray} \tag{ 15 }$$

The term ${[{\boldsymbol{F}}][{\boldsymbol{F}}]}^{+}$ can be interpreted as the projector on the column space of $[{\boldsymbol{F}}]$. Hence $({[{\boldsymbol{I}}]\mbox{--}[{\boldsymbol{F}}][{\boldsymbol{F}}]}^{+})$ is the orthogonal projector on the null space of $[{\boldsymbol{F}}]$. We note that ${[{\boldsymbol{s}}]}_{i}\in { \mathcal C }([{\boldsymbol{F}}])$, for $i=\{1,2,3,\ldots ,n\}-\{\alpha \}$. Hence, ${[{\boldsymbol{F}}][{\boldsymbol{F}}]}^{+}[{\boldsymbol{s}}{]}_{i}^{T}=[{\boldsymbol{s}}{]}_{i}^{T}$ i.e., ${[{\boldsymbol{W}}]}_{\alpha }[{\boldsymbol{s}}{]}_{i}^{T}=0$; however, since ${[{\boldsymbol{s}}]}_{0}\notin { \mathcal C }([{\boldsymbol{F}}])$, we have ${[{\boldsymbol{W}}]}_{\alpha }[{\boldsymbol{s}}{]}_{\alpha }^{T}=1$, which is the condition for the reconstruction of the component under consideration. Equation (15) is not directly usable by us since a priori we do not have any knowledge about the templates to be reconstructed from data and hence about the covariance matrix $[{\boldsymbol{F}}]$. We, however, see that the crucial properties of ${[{\boldsymbol{W}}]}_{\alpha }$, viz., ${[{\boldsymbol{W}}]}_{\alpha }[{\boldsymbol{s}}{]}_{i}^{T}=0$ and ${[{\boldsymbol{W}}]}_{\alpha }[{\boldsymbol{s}}{]}_{\alpha }^{T}=1$ that must be satisfied for reconstruction of the α component remain unchanged if we replace $[{\boldsymbol{F}}]$ in Equation (15) by a new ${n}_{b}\times n-1$ matrix, $[{\boldsymbol{G}}]$, with column vectors that are given by each of the shape vectors of the $n-1$ components, which are to be removed from the data. We therefore compute weights replacing $[{\boldsymbol{F}}]$ by $[{\boldsymbol{G}}]$ in Equation (15).

Once all components have been reconstructed, using weights following Equation (15), for the complete set of possible shape vectors of all components, we can form a reconstructed frequency map, ${\hat{Q}}_{i}(p)$, at each frequency, ${\nu }_{i}$, following

$$\begin{eqnarray}&&{\hat{Q}}_{i}(p)=\displaystyle \sum _{j=1}^{n}{\hat{T}}_{j}(p){s}_{i}^{j}.\end{eqnarray} \tag{ 16 }$$

Due to the underlying linearity in our method, one can clearly see that noise, ${\tilde{Q}}_{k}^{n}(p,{\nu }_{i})$, due to the kth reconstructed template at frequency ${\nu }_{i}$, has a contribution from noise of all input frequency maps. Specifically,

$$\begin{eqnarray}&&{\tilde{Q}}_{k}^{n}(p,{\nu }_{i})={s}_{i}^{k}\displaystyle \sum _{j=1}^{{n}_{b}}{w}_{k}^{j}{Q}_{j}^{n}(p).\end{eqnarray} \tag{ 17 }$$

Hence the total noise, ${\tilde{Q}}_{i}^{n}(p)$, due to all recovered components at frequency ${\nu }_{i}$, is given by,

$$\begin{eqnarray}{\tilde{Q}}_{i}^{n}(p) & = & \displaystyle \sum _{k=1}^{n}{\tilde{Q}}_{k}^{n}(p,{\nu }_{i})=\displaystyle \sum _{k=1}^{n}{s}_{i}^{k}\displaystyle \sum _{j=1}^{{n}_{b}}{w}_{k}^{j}{Q}_{j}^{n}(p)\\ & = & \displaystyle \sum _{j=1}^{{n}_{b}}{\tilde{w}}_{i}^{j}{Q}_{j}^{n}(p),\end{eqnarray} \tag{ 18 }$$

where we have defined

$$\begin{eqnarray}&&{\tilde{w}}_{i}^{j}=\displaystyle \sum _{k=1}^{n}{s}_{i}^{k}{w}_{k}^{j}.\end{eqnarray} \tag{ 19 }$$

We form a difference map at each frequency, ${\nu }_{i}$, by subtracting the recovered and input frequency map at the same frequency, ${\nu }_{i}$. Using Equation (18), noise in the ith difference map is given by

$$\begin{eqnarray}&&{Q}_{i}^{n}(p)-{\tilde{Q}}_{i}^{n}(p)=(1-{\tilde{w}}_{i}^{i}){Q}_{i}^{n}(p)-\displaystyle \sum _{j\ne i}^{{n}_{b}}{\tilde{w}}_{i}^{j}{Q}_{j}^{n}(p).\end{eqnarray} \tag{ 20 }$$

Assuming that the noise of different frequency bands has zero mean and that noise properties of different detectors are uncorrelated with one another, noise covariance of the ith difference map is given by

$$\begin{eqnarray}{N}_{{pp}^{\prime} }^{i,D} & = & {(1-{\tilde{w}}_{i}^{i})}^{2}\langle {Q}_{i}^{n}(p){Q}_{i}^{n}(p^{\prime} )\rangle \\ & & +\displaystyle \sum _{j\ne i}^{{n}_{b}}{({\tilde{w}}_{i}^{j})}^{2}\langle {Q}_{j}^{n}(p){Q}_{j}^{n}(p^{\prime} )\rangle \\ & = & {(1-{\tilde{w}}_{i}^{i})}^{2}{N}_{{pp}^{\prime} }^{i}+\displaystyle \sum _{j\ne i}^{{n}_{b}}{({\tilde{w}}_{i}^{j})}^{2}{N}_{{pp}^{\prime} }^{j},\end{eqnarray} \tag{ 21 }$$

where ${N}_{{pp}^{\prime} }^{i(j)}$ is the noise covariance matrix of the $i\mathrm{th}(j\mathrm{th})$ input frequency map.

Using the noise covariance matrix obtained in the previous section, we define the mismatch between the data and the model by means of the usual ${\chi }^{2}$ statistic, following

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{N}\displaystyle \sum _{i,p,p^{\prime} }({Q}_{i}(p)-{\tilde{Q}}_{i}(p)){({N}^{i,D})}_{{pp}^{\prime} }^{+}({Q}_{i}(p^{\prime} )-{\tilde{Q}}_{i}(p^{\prime} )),\end{eqnarray} \tag{ 22 }$$

where ${({N}^{i,D})}_{{pp}^{\prime} }^{+}$ represents the Moore–Penrose generalized inverse of ${N}_{{pp}^{\prime} }^{i,D}$ and $N={n}_{b}\times {{ \mathcal N }}_{\mathrm{pix}}$ denotes the total number of surviving pixels counting all frequency bands after any suitable mask has been applied on the maps. Since ${\tilde{Q}}_{i}(p)$ depends upon ${\beta }_{0s},{\beta }_{d}$ through Equation (16) and ${N}_{{pp}^{\prime} }^{i,D}$ depends upon ${\beta }_{0s},{\beta }_{d}$ and frequency band noise covariance matrices through Equation (21), we have ${\chi }^{2}\equiv {\chi }^{2}({\beta }_{0s},{\beta }_{d})$. Clearly, ${\beta }_{0s}$ and ${\beta }_{d}$ enter as the minimizing variables for the ${\chi }^{2}$ statistic. Our best-fit values for these spectral indices are the ones that together minimize the ${\chi }^{2}$ given by the above equation.

Although Equation (22) may resemble an ML estimation of spectral indices (assuming a uniform prior) for a given noise model, we note that our method is not an ML approach, since we solve only the first order perturbative expansion, e.g., as given in Equation (4), which is not an exact equation, instead of using a full and completely accurate power-law model as given in Equation (2). However, the perturbative expansion will tend to an ML approach asymptotically. Moreover, our method can be extended to more than merely first order in perturbation expansion, whenever it is allowed by the data, (i.e., when the instrument's noise level is relatively small and the total number of free parameters to be estimated is less than or equal to the total number of frequency maps available). At some maximum order of perturbative expansion (given the small variation of the synchrotron spectral index parameter), our results will practically be indistinguishable from the true ML results within the noise limit of instruments. So sufficiently higher order perturbative analysis will possess the desirable property that its results will be identical in terms of both parameter values and error bars derived on them with the ML results, within the instrument's noise uncertainty.

5.1.1. Frequency bands

We include simulated Stokes Q polarization maps of WMAP K, Ka, Q, V bands, all three Planck LFI bands, and three lowest frequency Planck HFI bands. We do not include the WMAP W band in our analysis since it contains a somewhat larger noise level compared to other WMAP bands. We do not include the two highest frequency Planck HFI maps since they only measure temperature. Moreover, with the inclusion of the HFI 353 GHz map in our analysis, ${\chi }^{2}$ values become increasingly more sensitive to the thermal dust model, giving lesser weightage to the synchrotron component. Hence we also exclude this frequency band from our analysis. Detailed specifications of the frequency bands used in this work are given in Table 1.

Table 1.  Frequency Maps

Frequency, ν	WMAP/Planck	$a(\nu )$	${\sigma }_{0}(\mu {\rm{K}})$	${\sigma }_{w}(\mu {\rm{K}})$
(GHz)	Map Name		Ref.a	Ref.a
23	WMAP K	1.014	40.86	1435.0
30	Planck, LFI	1.023	36.56	⋯
33	WMAP Ka	1.029	41.19	1472.0
41	WMAP Q1	1.044	56.30	2254.0
41	WMAP Q2	1.044	53.41	2140.0
44	Planck, LFI	1.051	37.03	⋯
61	WMAP V1	1.099	70.86	3324.0
61	WMAP V2	1.099	63.05	2958.0
71	Planck, LFI	1.136	37.11	⋯
100	Planck, HFI	1.283	15.83	⋯
143	Planck, HFI	1.646	11.80	⋯
217	Planck, HFI	2.961	19.39	⋯

Note. List of WMAP and Planck frequency maps used in the analysis. The third column denotes antenna to thermodynamic conversion factors at different frequencies. In the last column, ${\sigma }_{w}$ values represent the polarized noise level of WMAP detectors for a single observation. In the fourth column, ${\sigma }_{0}$ values describing Planck polarized noise levels are consistent with the Planck blue book (Consortium 2005) values and scaled to ${N}_{\mathrm{side}}=512$. For WMAP ${\sigma }_{0}$, values represent ${\sigma }_{w}/\sqrt{\langle {N}_{\mathrm{obs}}(p)\rangle }$ where $\langle \rangle$ denotes the all-sky average and ${N}_{\mathrm{obs}}(p)$ represents the effective number of observations at a given pixel p for a given ${N}_{\mathrm{side}}=512$ map. The fourth column clearly shows a lower noise level of Planck compared to WMAP. Wherever multiple differencing assemblies are available for WMAP, a polarized noise level for each differencing assembly is also mentioned in order. For WMAP, the actual noise level per pixel at ${N}_{\mathrm{side}}=512$ is given by $\sigma (p)={\sigma }_{w}/\sqrt{{N}_{\mathrm{obs}}(p)}$. For a Planck frequency, the band noise level per pixel at ${N}_{\mathrm{side}}=512$ is assumed to be uniform with $\sigma (p)={\sigma }_{0}$ as the specified values at ${N}_{\mathrm{side}}=512$.

^aStokes (Q, U) polarization noise level in $\mu {\rm{K}}$ thermodynamic temperature unit at ${N}_{\mathrm{side}}=512$.

Download table as:  ASCIITypeset image

Since pixel values of spectral index variation map, ${\rm{\Delta }}{\beta }_{s}(p)$, are small compared to the foreground emission, any significant residual noise in the reconstructed template maps results in its potentially noisy reconstruction following Equation (8). Owing to the delicate noise sensitivity of the reconstructed spectral index map, we chose to minimize the detector noise at the beginning of the analysis. For this purpose, we use Stokes Q frequency maps at ${N}_{\mathrm{side}}=8$, smoothed by a Gaussian beam function of FWHM = 20°.

5.1.2. Mask

5.1.3. Physical Components

5.1.4. Spectral Index Map

5.1.5. Detector Noise Maps and Covariance Matrix

Generating noise maps and thus modeling the noise properties in the low-resolution maps are somewhat more involved than the simple smoothing operations that needed to be performed on foreground templates and CMB maps at ${N}_{\mathrm{side}}=8$ in Section 5.1.3. We take into account WMAP's scan-modulated detector noise pattern by using ${N}_{\mathrm{obs}}(p)$ information associated with each nine-year Differencing Assembly (DA) map available on the LAMBDA website at ${N}_{\mathrm{side}}=512$. We first compute noise variance at each pixel following ${\sigma }^{2}(p)={\sigma }_{w}^{2}/{N}_{\mathrm{obs}}(p)$ at ${N}_{\mathrm{side}}=512$ for all DA maps between the K to V bands. Since for each of the Q and V bands two DA maps are available, we simply weigh down the average variance of these two DA maps at each pixel, by a factor of two, to obtain noise variance maps corresponding to an average DA map within the frequency bands. For Planck LFI and HFI maps, we use a simple uniform white noise model. Pixel white noise levels, $\sigma (p)$, for Planck maps are obtained from Planck blue book specifications, ${\sigma }_{\mathrm{bb}}$, which represents noise in the beam area of the corresponding map in units of $\mu {\rm{K}}/{\rm{K}}$. Pixel noise levels for Planck maps at ${N}_{\mathrm{side}}=512$ are then given by, $\sigma (p)={\sigma }_{\mathrm{bb}}{T}_{\mathrm{CMB}}{B}_{\mathrm{fwhm}}/{\rm{\Delta }}\theta$, where B_fwhm represents the FWHM of Planck beams in arcmin, ${\rm{\Delta }}\theta =6\buildrel{\,\prime}\over{.} 9$, represents the pixel size of the ${N}_{\mathrm{side}}=512$ map and ${T}_{\mathrm{CMB}}=2.73\,{\rm{K}}$.

In principle, one can use all the noise variance maps at ${N}_{\mathrm{side}}=512$ corresponding to all frequency bands of Table 1 to obtain representations of noise maps at ${N}_{\mathrm{nside}}=512$ and then downgrade them to ${N}_{\mathrm{side}}=8$ for use in our simulation. However, one needs to simulate many such maps at ${N}_{\mathrm{side}}=512$, since we need to obtain pixel noise covariance at ${N}_{\mathrm{side}}=8$ by downgrading the high-resolution noise maps. Since generating noise maps at ${N}_{\mathrm{side}}=512$ and subsequently downgrading them many times is a computationally expensive process, which requires $\sim {N}_{\mathrm{pix}}\,$ $=\,12\times 512\times 512$ of floating point operations per simulation, we devise a completely equivalent approach that is computationally significantly fast compared to the former method. Instead of downgrading ${N}_{\mathrm{side}}=512$ noise maps, we obtain noise variance maps at ${N}_{\mathrm{side}}=8$ by simply multiplying the total noise variance of all ${N}_{\mathrm{side}}=512$ pixels lying inside each ${N}_{\mathrm{side}}=8$ pixel by the ratio of pixel numbers of a HEALPix map at ${N}_{\mathrm{side}}=8$ and ${N}_{\mathrm{side}}=512$. Such a method is always feasible for our method and does not introduce any pixel–pixel correlation at ${N}_{\mathrm{side}}=8$ since noise in the original ${N}_{\mathrm{side}}=512$ maps is uncorrelated from pixel to pixel. Using these noise variance maps, we simulate noise maps at ${N}_{\mathrm{side}}=8$. These noise maps are then smoothed by a Gaussian beam of FWHM = 20°, to generate realistic noise maps at ${N}_{\mathrm{side}}=8$ for each frequency band. We obtain 10000 Monte-Carlo simulations of smoothed noise maps at ${N}_{\mathrm{side}}=8$. Since the smoothing operations introduce pixel to pixel noise covariance, we obtain a full pixel to pixel noise covariance matrix over surviving pixels (after application of pm10 mask) at ${N}_{\mathrm{side}}=8$ corresponding to each frequency band. We have shown a noise covariance matrix for the V-band masked map, in Figure 1. The non-diagonal nature of this matrix indicates non-trivial noise properties in the maps.

Zoom In Zoom Out Reset image size

Using the ${N}_{\mathrm{side}}=8$ smoothed CMB Stokes Q polarized map, Q_c(p), synchrotron and thermal dust templates (Q_0s and Q_0d respectively), spectral index map, ${\beta }_{s}(p)$, and detector noise map, ${Q}_{i}^{n}(p)$ ($i=1,2,\ldots ,{n}_{b}$) as obtained above, we generate simulated frequency maps for WMAP and Planck frequency bands, using Equation (2) and knowing the values of a_i. We use ${\beta }_{d}=2.0,{\nu }_{0s}=23\,\mathrm{GHz}$ and ${\nu }_{0d}=94\,\mathrm{GHz}$ in this work.

Since spectral indices of neither synchrotron, nor thermal dust components are known a priori, we perform a reconstruction of all components defined by Equation (5) for a suitable range of both of these indices. We vary ${\beta }_{0s}$ between −2.65 to −3.05 with a step of −0.01. For each of these values ${\beta }_{d}$ is varied between 1.80 and 2.20 with a step size of 0.01. For each pair (${\beta }_{0s},{\beta }_{d}$), such generated, we compute shape vectors for all foreground template components. Combining these with the shape vector of CMB, we estimate a set of weights using Equation (15) to reconstruct all template components including CMB. For each pair (${\beta }_{0s},{\beta }_{d}$), we form reconstructed frequency maps following Equation (16) using all the reconstructed template components. We then form the difference map at each frequency by subtracting a reconstructed frequency map from the corresponding input frequency map. The noise covariance of the ith difference map over the surviving pixels for a given pair of indices is given by the matrix, ${N}_{{pp}^{\prime} }^{i,D}$, defined by Equation (21). Once the difference map is formed, we obtain ${\chi }^{2}({\beta }_{0s},{\beta }_{d})$ for all indices following Equation (22). Our best-fit values for indices are the ones that minimize ${\chi }^{2}$ over all indices. Using these values, we obtain a spectral index variation map following Equation (8). We validate our method by performing 1000 Monte-Carlo simulations of the entire procedure described above. We note that since CMB follows a blackbody distribution its shape vector remains identical for all foreground indices and for all simulations. Since the weights, as given by Equation (15), depend only upon foreground and CMB shape vectors (without any dependence on random detector noise, which varies from simulation to simulation) for the purpose of reconstruction of all components, we calculate these weights only once for the entire range of spectral indices and store them on the disk. We then use these weights to reconstruct all components for all Monte-Carlo simulations for all indices. Moreover, we also precompute Moore–Penrose generalized inverses of detector noise covariance matrices of all difference maps for all pairs of spectral indices $({\beta }_{0s},{\beta }_{d})$ and store them on the disk. In effect, a Singular Value Decomposition algorithm with a cutoff of $1.0\times {10}^{-7}$ was used to implement Moore–Penrose inverses. An alternative approach would have been to invert these covariance matrices after adding a small regularization noise to each one of them.

For each simulation described in Section 5, we construct a spectral index variation map, ${\rm{\Delta }}{\hat{\beta }}_{s}(p)$, over the reference index, ${\beta }_{0s}$, following Equation (8). The results for reconstruction of spectral index map are shown in Figure 2. The colored area of the Mollweide projection of all the three panels shows the region on the sky that survives after applying the pm10 mask. The top panel shows the mean reconstructed synchrotron spectral index map, $\langle {\hat{\beta }}_{s}(p)\rangle$, defined as, $\langle {\hat{\beta }}_{s}(p)\rangle =\langle {\beta }_{0s,\min }\rangle +\langle {\rm{\Delta }}{\hat{\beta }}_{s}(p)\rangle$, where ${\beta }_{0s,\min }$ denotes a value of ${\beta }_{0s}$ (e.g., see Equation (4)) corresponding to ${\chi }_{\min }^{2}$ and the average is performed over all 1000 simulations. The middle panel is a simple standard deviation map for any one of the reconstructed spectral index maps. The lowermost panel shows the difference between the average reconstructed index map and the input spectral index map, $\langle {\hat{\beta }}_{s}(p)\rangle \,\mbox{--}\,{\beta }_{s}(p)$. Clearly, the difference map shows that our method can accurately reconstruct the spectral index map inside the pm10 mask. The reconstruction difference toward higher (lower) galactic coordinates becomes slightly larger due to weak polarized signal in these regions.

Zoom In Zoom Out Reset image size

For any given Monte-Carlo simulation, ${\chi }^{2}({\beta }_{0s},{\beta }_{d})$ values show a clear minimum, ${\chi }_{\min }^{2}$, at some values of $({\beta }_{0s},{\beta }_{d})=({\beta }_{0s,\min },{\beta }_{d,\min })$. Although a topographical plot of ${\chi }^{2}({\beta }_{0s},{\beta }_{d})$ depends upon the particular random simulation under consideration, the features of such plots that arise on average due to foregrounds, CMB, and detector noise, can be studied by making a plot of $\langle {\chi }^{2}({\beta }_{0s},{\beta }_{d})\rangle$ (where the average is performed over all 1000 simulations for each pair $({\beta }_{0s},{\beta }_{d})$) with respect to both indices. We have shown this at the top panel of Figure 3. The elongated contours of this plot show the relatively higher sensitivity of the data with the variation of synchrotron spectral index compared to the same variation of the dust spectral index. In the bottom panel of this figure, we show a standard deviation map of ${\chi }^{2}$ values for all indices. Lower standard deviation values around $\langle {\chi }_{\min }^{2}\rangle$ indicate that the minimum is recovered satisfactorily by our method. Using results from Monte-Carlo simulations, we estimate an average value of $({\beta }_{0s,\min },{\beta }_{d,\min })$ corresponding to ${\chi }_{\min }^{2}$ as $(-2.84\pm 0.01,2.00\pm 0.004)$. A small error on the spectral index values corresponding to ${\chi }_{\min }^{2}$ shows that our method can reliably recover global spectral indices.

Zoom In Zoom Out Reset image size

Using the Monte-Carlo simulations, we find that values of ${\chi }_{\min }^{2}$ are distributed over a narrow range of 0.82–1.05. A plot of the histogram of ${\chi }_{\min }^{2}$ is shown in the top panel of Figure 4. We estimate, $\langle {\chi }_{\min }^{2}\rangle =0.93$, averaging over 1000 Monte-Carlo simulations, with a sample standard deviation of 0.04. The probability density function of this figure follows a very symmetric distribution for ${\chi }_{\min }^{2}$. This can be expected due to a very large number of degrees of freedom (∼4260) available in the analysis. In the bottom panel of this figure, we have shown the probability distribution of 20° smoothed ${N}_{\mathrm{side}}=8$ input spectral index map, from the region that survives after the application of the pm10 mask. A bimodal structure of the distribution is readily identifiable corresponding to spectral indices of $\sim -2.85$ and −3.16 respectively. Comparing these values with the mean reconstructed spectral index −2.84 ± 0.01, it is interesting to note that our method picks out the peak corresponding to a flatter index from the available pair of doubly degenerate ones.

Zoom In Zoom Out Reset image size

One very important component of this work is the weights. We show weights for reconstruction of ${Q}_{0s}(p)$ and ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$ templates for the entire range of spectral indices, respectively, in the top and bottom panels of Figure 5. At the low-frequency side, where synchrotron emission is dominant, weights, for reconstruction of the synchrotron template, tend to be positive for all indices. We see this for the K, LFI 30 GHz, and Ka bands. For all other bands except the HFI 217 GHz band, weights are negative for all indices. HFI 217 GHz band weights become positive for all indices to cancel out signals from CMB, thermal dust, ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$, and at the same time to preserve Q_0s. For the reconstruction of ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$, template weights must completely project out its shape vector from the data, at the same time, the weight vector needs to be orthogonal to shape vectors of synchrotron and other components. Because of the close analytical similarity in the shape of vectors ${Q}_{0s}(p)$ and ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$, the orthogonality condition leads to a somewhat larger magnitude of weights for ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$ than the ${Q}_{0s}(p)$ template. The relatively larger magnitude of weights can more easily be seen from the second panel than the first one. For the second panel and for the K band, weights take the largest negative value, since, from the expression of the shape vector for the ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$ component, ($a(\nu ){(\nu /{\nu }_{0s})}^{-{\beta }_{0s}}\mathrm{ln}(\nu /{\nu }_{0s})$), the shape vector has a zero component at the K band. For low-frequency LFI 30 GHz, Ka, Q, and LFI 44 GHz frequency bands, all weights are positive, indicating a larger shape vector at lower frequency. For V, LFI 71 GHz, HFI 100 GHz, and HFI 143 GHz frequency bands weights become negative, since shape vector decreases in this frequency range. Finally, for HFI 217 GHz, frequency band weights again become positive to project out the ${Q}_{0s}(p){\rm{\Delta }}{\beta }_{s}(p)$ component completely and to cancel out the rest of the components.

Zoom In Zoom Out Reset image size

As associated results of our method, we reconstruct asymptotically ML templates for all components using the best-fit spectral index values ${\beta }_{0s,\min },{\beta }_{d,\min }$ corresponding to ${\chi }_{\min }^{2}$ obtained in Section 6.1. We show the results in Figure 6. Unlike the case for reconstructed synchrotron spectral map, we do not apply a mask on the recovered templates, since CMB and foreground signals are expected to be stronger than the spectral index variations. Top, middle, and bottom panels of this figure show results for CMB, thermal dust, and synchrotron template reconstruction, respectively, obtained from a randomly chosen Monte-Carlo simulation with seed = 100. All the color maps of this figure are in $\mu {\rm{K}}$ units (thermodynamic for CMB, antenna for foreground templates). Following the usual row, column convention of a matrix, we represent any image of this figure by the pair (i, j). (1, 1) image represents the reconstructed CMB map against the input CMB map for the same simulation, shown as the (1, 2) image. As one can easily visually identify, the large-scale features of these two maps are very nearly identical. The (1, 3) image shows the difference between the reconstructed and the input CMB map and hence is dominated by residual detector noise. The noisy nature of the difference map is also indicated by the symmetric color table plotted below the (1, 3) image. The (1, 4) image represents the standard deviation map estimated from the set of all (1, 3) maps resulting from 1000 Monte-Carlo simulations of the component separation algorithm. As expected, the standard deviation map clearly indicates a scan-modulated noise pattern as applicable for WMAP observations. There exist some foreground induced errors on the galactic plane. Such errors can be reduced by extending our method to a second-order perturbative analysis.

Zoom In Zoom Out Reset image size

The (2, 1) image represents the reconstructed thermal dust template at 217 GHz for seed = 100. The reconstructed template matches very well with the input thermal dust template at the same frequency, viz., the (2, 2) image. The difference between the reconstructed and input templates ((2, 3) image) is dominated by noise. The standard deviation map for reconstructed thermal dust template is shown as the (2, 4) image. Clearly, WMAP's scan-modulated noise pattern is visible in this figure. There also exist some residual foreground errors in this map, which again can be eliminated by performing a higher order perturbative analysis.

The (3, 1) image shows the reconstructed synchrotron template at 23 GHz for the same random seed. The input synchrotron template at 23 GHz is shown as the (3, 2) image. Both these maps match with each each other very well. The difference of the first and second maps is shown as the (3, 3) image. The standard deviation map, image (3, 4), shows the scan-modulated detector noise pattern, without any visible signature of residual foreground error.

The normalized reconstruction error in each of the above reconstructed templates can be quantified by dividing the difference of the output reconstructed template and the input one by the corresponding standard deviation map (i.e., dividing each map of column 3 by the corresponding map of column 4 of Figure 6). The dimensionless pixel values of these error maps should have zero mean and unity standard deviation for reliable template reconstruction. In Figure 7, we show these reconstruction error maps for CMB (top panel), thermal dust (middle panel), and the synchrotron template (bottom panel). We have shown the mean, $\hat{\mu }$, and standard deviation, ${\hat{s}}_{d}$, values for CMB, thermal dust, and synchrotron templates in Table 2 for the entire sky, pm10 mask, and from the complementary region of the pm10 mask. Comparing the mean values, we conclude that, among all of these three regions, the reconstruction of CMB and thermal dust template components are most accurate for the pm10 sky region. The synchrotron template also has a small mean-error on the same sky region, compared to the maximum and minimum of the input template, indicating an accurate reconstruction. Furthermore, the minimum and maximum pixel values obtained from the entire sky for each of these maps are, respectively, (−2.79, 2.96), (−2.61, 2.92) (−2.92, 3.11), indicating a faithful reconstruction of input templates.

Zoom In Zoom Out Reset image size

Table 2.  Mean $\hat{\mu }$, and Standard Deviation, ${\hat{s}}_{d}$, Statistics of Normalized Error Maps

Component Name	Fullsky		pm10 Mask		pm10c Mask
	$\hat{\mu }$	${\hat{s}}_{d}$	$\hat{\mu }$	${\hat{s}}_{d}$	$\hat{\mu }$	${\hat{s}}_{d}$
CMB	0.169	0.995	0.017	0.913	0.359	1.060
Thermal Dust	−0.087	1.016	0.036	0.953	−0.240	1.070
Synchrotron	−0.210	0.947	−0.294	0.918	−0.106	0.974

Note. Mean and standard deviation values estimated from the normalized error maps of different template components. The pm10c mask is complementary to the pm10 mask.

Download table as:  ASCIITypeset image