An improved model of diffuse galactic radio emission from 10 MHz to 5 THz

Abstract

We present an improved Global Sky Model (GSM) of diffuse Galactic radio emission from 10 MHz to 5 THz, whose uses include foreground modelling for cosmic microwave background (CMB) and 21 cm cosmology. Our model improves on past work both algorithmically and by adding new data sets such as the Planck maps and the enhanced Haslam map. Our method generalizes the principal component analysis approach to handle non-overlapping regions, enabling the inclusion of 29 sky maps with no region of the sky common to all. We also perform a blind separation of our GSM into physical components with a method that makes no assumptions about physical emission mechanisms (synchrotron, free–free, dust, etc). Remarkably, this blind method automatically finds five components that have previously only been found ‘by hand’, which we identify with synchrotron, free–free, cold dust, warm dust, and the CMB anisotropy. Computing the cross-power spectrum between these blindly extracted components and Planck component maps, we find a strong correlation at all angular scales. The improved GSM is available online at http://github.com/jeffzhen/gsm2016.

radiation mechanisms: general, methods: data analysis, techniques: interferometric, cosmic background radiation, radio continuum: general, radio lines: general

1 INTRODUCTION

Modelling diffuse Galactic radio emission has received great interest in cosmology. Cosmological signals, such as the cosmic microwave background (CMB) and redshifted 21 cm emission, have to pass through our Galaxy before reaching us, so they are inevitably mixed with foreground contamination from Galactic emission. For CMB experiments, the unpolarized CMB signal dominates over foregrounds in out-of-plane directions, though the situation for polarized CMB signal is about one to two orders of magnitude worse (for reviews, see de Oliveira-Costa & Tegmark 1999 and references therein). For the redshifted 21 cm signals from the epoch of reionization (EoR), the foregrounds are thought to be four or more orders of magnitudes higher than the 21 cm signal (see Furlanetto, Oh & Briggs 2006; Morales & Wyithe 2010 for reviews, and Dillon et al. 2014; Parsons et al. 2014; Ali et al. 2015 for the latest 21 cm power spectrum upper limits), which makes foreground modelling even more important (Datta, Bowman & Carilli 2010; Morales et al. 2012; Parsons et al. 2012; Trott, Wayth & Tingay 2012; Vedantham, Udaya Shankar & Subrahmanyan 2012; Pober et al. 2013, 2016; Thyagarajan et al. 2013, 2015; Liu, Parsons & Trott 2014a,b; Dillon et al. 2015; Pober 2015; Mozdzen et al. 2016). Among many efforts to model the foreground are the Planck Sky Model (PSM; Delabrouille et al. 2013) and the GSM.

In the original GSM (de Oliveira-Costa et al. 2008), the authors carried out an exhaustive survey of existing sky maps from 10 MHz to 100 GHz, and performed a principal component analysis (PCA) on 11 of the highest quality maps. They found that using just the top three principal components could explain between 90 and 99 per cent of the variations in all the maps, depending on frequency and sky direction. By interpolating these three components over frequency, the GSM can therefore model the diffuse Galactic emission anywhere between 10 MHz and 100 GHz.

The considerable flexibility and power of the GSM software has resulted in its widespread application. Given that foregrounds are not only ubiquitous in cosmological surveys, but also interesting probes of astrophysical phenomena in their own right, it is unsurprising that the GSM has had broad influence within the astrophysics community.

Astrophysical research using the GSM has ranged, for example, from limits on radio emission from annihilating dark matter (e.g. Spekkens et al. 2013) to searches for non-Gaussianity in Planck CMB maps (e.g. Novaes et al. 2015) and from EoR power spectrum limits (e.g. Ali et al. 2015) to the discovery of a Fast Radio Burst (e.g. Burke-Spolaor & Bannister 2014).

However, there are some notable limitations to the GSM. In terms of frequency coverage, there is a lack of direct observation in the EoR frequency range. The closest maps in frequency are a 45 MHz map with about 3| $_{.}^{\circ}$ |6 resolution (Alvarez et al. 1997; Maeda et al. 1999) and a 408 MHz map with 56 arcmin resolution (Haslam et al. 1981, 1982; Remazeilles et al. 2015), and for the high-resolution GSM version, the entire model is locked to the 408 MHz map. In terms of sky coverage, the PCA algorithm is performed on a rather small region covered by all 11 maps, mostly in out-of-plane directions, so the resulting principal components may suffer from biased frequency dependences skewed towards that of the small region covered. In terms of accuracy, each of the 11 data sets are normalized beforehand so that they are treated equally, regardless of their relative accuracy. Thus, the model may suffer unnecessarily from noise and systematics in the maps of lesser quality. Lastly, since the PCA produces orthogonal principal components and orthogonal principal maps, these components are not likely to correspond to actual physical processes such as synchrotron or dust emission. Actual emission mechanisms have similar spatial structures (such as strong emission in the Galactic plane) and thus are not mutually orthogonal.

In this work, we present a new GSM-building method that naturally extends the original PCA algorithm. We use PCA as the initial step to obtain crude estimates of the principal components and their corresponding maps, and we iterate between the components and the maps to find the best fit to all of the data available. This method allows us to include 29 sky maps in the frequency range 10 MHz to 5 THz (including the Parkes maps at 85 MHz and 150 MHz; Landecker & Wielebinski 1970), which share no completely common sky coverage. This method can be further extended in future to incorporate uncertainty information from each map into the fitting process. Furthermore, we use a blind component separation technique to recombine the orthogonal components into those that are physically interpretable as different emission mechanisms. Both our blind spectra and blind component maps agree remarkably well with existing physical models, and additionally are consistent with the Planck component maps.

The remainder of this paper is structured as follows. In Section 2, we describe our improved GSM-building method. In Section 3, we describe the 29 sky maps included in this work. In Section 4, we present the improved GSM from 10 MHz to 5 THz with six components in Section 4.1, and estimate the predictive accuracy of the new GSM in Section 4.2. In Section 5, we use a blind approach to recombine those components into physically meaningful contributions in Section 5.1, and compare our blind spectra and component maps with existing results in the literature in Sections 5.2 and 5.3, along with a quantitative comparison to Planck component maps in Section 5.4. Finally, in Section 6, we conclude by summarizing this work and discussing potential future improvements for the GSM.

2 ITERATIVE ALGORITHM FOR BUILDING A GSM

2.1 Framework

We describe the GSM in the form of two matrices, an n_pix × n_c map matrix | $\bf{M}$ | and an n_c × n_f normalized spectrum matrix | $\bf{S}$ |⁠, where n_pix is the number of pixels in each sky map, n_c is the number of components, and n_f is the number of frequencies for which we have maps. Furthermore, we encapsulate all the sky map data into an n_pix × n_f matrix | $\bf{D}$ |⁠, where all maps are normalized to the same level at each frequency.¹ The GSM then models the sky by

\begin{equation} \bf{D}\approx \bf{M}\bf{S}. \end{equation}

(1)

Thus to construct a GSM is to find the pair of | $\bf{M}$ | and | $\bf{S}$ | that minimizes the cost function

\begin{equation} \Vert \bf{D}-\bf{M}\bf{S}\Vert ^2. \end{equation}

(2)

In general, because both components in the product | $\bf{M}\bf{S}$ | are unknown, we have degeneracies in the form of an invertible n_c × n_c matrix Ψ. For any solution | $\bf{M}$ | and | $\bf{S}$ |⁠, an alternative solution | $\bf{M}^{\prime } = \bf{M}\mathbf {\Psi }^{-1}$ | and | $\bf{S}^{\prime } = \mathbf {\Psi }\bf{S}$ | will produce an identical prediction for the data | $\bf{D}$ |⁠.

To use the GSM to predict sky maps at previously unmeasured frequencies or sky regions, one interpolates at the desired frequency both the normalized spectra in | $\bf{S}$ | and the overall normalization. For the rest of this work, we choose linear interpolation for the normalized spectra over the log of frequency, and linear interpolation for the log of normalization over the log of frequency. Because normalized spectra are interpolated linearly, the choice of degeneracy matrix Ψ will have no effect on the predictions of the GSM.

2.2 PCA algorithm

We first describe the original PCA algorithm in de Oliveira-Costa et al. (2008), before generalizing it to our iterative algorithm. Given the data matrix | $\bf{D}$ |⁠, one first performs an eigen-decomposition of | $\bf{D}^t\bf{D}$ |⁠:

\begin{equation} \bf{D}^t\bf{D}= \bf{C}^t\mathbf {\Lambda }\bf{C}, \end{equation}

(3)

where | $\bf{D}^t$ | denotes the transpose of | $\bf{D}$ |⁠, | $\bf{C}$ | is an n_f × n_f orthogonal matrix with eigenvectors as its rows, and | $\mathbf {\Lambda }$ | is an n_f × n_f diagonal matrix with eigenvalues on its diagonal. If the sky can be described by n_c components where n_c < n_f, then | $\mathbf {\Lambda }$ | has only n_c non-zero eigenvalues on its diagonal,² so

\begin{equation} \bf{D}^t\bf{D}= \widetilde{\bf{C}}^t\widetilde{\mathbf {\Lambda}}\widetilde{\bf{C}}, \end{equation}

(4)

where | $\widetilde{\bf{C}}$ | and | $\widetilde{\mathbf {\Lambda }}$ | are the n_c × n_f and n_c × n_c parts of their non-tilde counterparts corresponding to the non-zero eigenvalues. One then takes the principal components | $\widetilde{\bf{C}}$ | and solves for the best | $\bf{M}$ | that satisfies

\begin{equation} \bf{M}\widetilde{\bf{C}} \approx \bf{D}. \end{equation}

(5)

Comparing equations (1) and (5), we see that the PCA algorithm obtains one solution to equation (1), with | $\bf{S} = \widetilde{\bf{C}}$ |⁠. However, in practice, the sky maps have various sky coverages with potentially very few overlapping pixels, so rather than using the full | $\bf{D}$ | for equation (3), the original GSM uses a | $\bf{D}^{\ast }$ | that consists of a subset of the rows in | $\bf{D}$ | that corresponds to the common pixels covered at all frequencies. Therefore, the process of obtaining | $\widetilde{\bf{C}}$ | is ignoring the majority of information in | $\bf{D}$ |⁠, and the solutions are thus not minimizing the cost function in equation (2). Furthermore, if one wishes to include more data sets, such as we do in this work, there are no overlapping regions that cover all of the frequency range, so it is difficult to apply the original PCA algorithm as-is.

2.3 Iterative algorithm

To overcome the challenges discussed above, we extend the PCA algorithm by iterating on the results obtained from it. We start by temporarily excluding the data sets with smallest sky coverage one by one, until the remaining data sets have more than 5 per cent common sky coverage. We then use the 5 per cent common pixels to obtain n_c, | $\bf{M}^{(0)}$ |⁠, and | $\widetilde{\bf{C}}$ | following the PCA algorithm described in the previous section, where | $\bf{M}^{(0)}$ | denotes our starting | $\bf{M}$ | at the 0th iteration. Once we obtain | $\bf{M}^{(0)}$ |⁠, we fold the temporarily excluded frequencies back into | $\bf{D}$ | and start iterating. At the ith iteration we compute

\begin{equation} \bf{S}^{(i)} = (1-\eta )\bf{S}^{(i-1)} + \eta (\bf{M}^{(i-1)t}\bf{M}^{(i-1)})^{-1}\bf{M}^{(i-1)t}\bf{D}, \end{equation}

(6)

\begin{equation} \bf{M}^{(i)} = (1-\eta )\bf{M}^{(i-1)} + \eta (\bf{S}^{(i)}\bf{S}^{(i)t})^{-1}\bf{S}^{(i)}\bf{D}^t, \end{equation}

(7)

where 0 < η ≤ 1 is a step size. Intuitively, at each iteration we are computing the least square solution | $\bf{S}^{(i)}$ | to equation (1), treating | $\bf{M}^{(i-1)}$ | as the truth, and vice versa. We keep iterating until the cost function decreases by less than 0.01 per cent.

Due to incomplete sky coverage, there are two further tweaks to the above equations. First, we do not directly compute the above equations in a matrix form. Rather, we compute | $\bf{S}^{(i)}$ | column by column, with each column representing one frequency, and we modify | $\bf{M}^{(i-1)}$ | to exclude the pixels not covered at that frequency. Similarly for equation (7), each column corresponds to a pixel, and we modify | $\bf{S}^{(i)}$ | to exclude frequencies that do not cover that pixel. Thus, at every iteration, the algorithm incorporates all pixel data in | $\bf{D}$ |⁠, regardless of the sky coverage at each frequency or the frequency coverage of each pixel. In addition, since the maps have different sky coverages, simply normalizing valid pixels in each map will overweigh the maps covering low temperature out-of-plane regions compared to those that only cover the Galactic Centre. Thus, we re-normalize the data matrix | $\bf{D}$ | at every iteration by dividing each map by the norm (root-sum-square) of its best-fitting map | $\bf{M}^{(i)}\bf{S}^{(i)}$ |⁠.

2.4 Incorporating more general data formats

Equation (1) assumes that all data sets contained in | $\bf{D}$ | are sky maps with the same pixelization and angular resolution, which is what we focus on in this work. However, this comes with the drawback of forcing us to work with the lowest common resolution for all maps, thus discarding the high-frequency information in higher quality maps. Furthermore, many low-frequency interferometers such as MWA (Tingay et al. 2013) and PAPER (Parsons et al. 2010, 2014; Ali et al. 2015) have produced high-quality data, but because they sample sparsely in the Fourier domain, it is challenging to reduce those measurements into pixelized sky maps. Fortunately, we can generalize our algorithm and allow all data sets that are linearly related to the pixelized sky maps. The data set at frequency index f can be described by

\begin{equation} D_{if} = \sum _{jn}B_{ij}^f M_{jn} S_{nf}, \end{equation}

(8)

which is a generalized version of equation (1). At any given frequency index f, | $\bf{B}^f$ | is a known matrix describing the linear relationship between the pixelized sky map and the data set. For pixelized sky-maps, | $\bf{B}^f$ | is the identity matrix and equation (8) reduces to equation (1). For low-resolution data sets, | $\bf{B}^f$ | contains the point spread function for all pixels in the sky, which is typically the antenna beam pattern. For drift-scanning interferometer data sets, | $\bf{B}^f$ | contains the primary beam pattern as well as the rotating fringe patterns for all baselines, and D_fi at f is a flattened list of visibilities over both time and baselines. Given | $\bf{B}$ | and | $\bf{D}$ |⁠, one can iteratively solve for | $\bf{S}$ | and | $\bf{M}$ |⁠. It is worth noting that while including the | $\bf{B}$ |-matrix does not significantly increase the amount of computations for obtaining | $\bf{S}$ | at each iteration, it makes the computation of | $\bf{M}$ | much more demanding. This is because one can no longer use equation (7) to compute | $\bf{M}$ | pixel by pixel, as the | $\bf{B}$ | mixes pixels in the sky. One now needs to operate matrices with n_pix × n_c rows rather than n_c rows as in equation (7). For a HEALPix map (Górski et al. 2005) with n_side = 64 corresponding to pixel size of 1° and six principal components, the matrix size is about 3 × 10⁵ on each side, which demands significant computing resources beyond the scope of this work.

3 SKY SURVEY DATA SETS

Table 1 and Fig. 1 list all sky maps we use in this work, ranging from 10 MHz to 5 THz. We have included all sky maps in this frequency range that have angular resolutions better than 5° and sky coverage larger than 20 per cent of the full sky. All data sets are publicly available online,³ except those at 45 MHz and 2.33 GHz, which are available on request from the respective data owners. We manually mask out the ecliptic plane in the 3.33 THz AKARI data and the 5 THz IRIS data, and exclude IRIS data above 5 THz, in order to remove zodiacal contamination. Unlike the original GSM, which used physical models to pre-remove the CMB anisotropy from the WMAP and Planck maps, we choose to let the data speak for themselves and not pre-remove the CMB anisotropy. These 29 maps form regions with 120 different combinations of map overlap, with no common region between all maps (see the bottom-right plot in Fig. 1), so the improved GSM-making method is necessary to combine all these data sets.

Figure 1.

29 sky maps used in this work from 10 MHz to 5 THz, plotted in Galactic coordinate under Mollweide projection centred at the Galactic Centre, on arcsinh scales, where the constant for each arcsinh is set to the overall amplitude of each map, as shown in Fig. 2. The colour scale follows the rainbow order, with red being the highest, purple the lowest, and grey signifying no data. The 11 bold and underscored frequencies are those included in the original GSM. The last panel shows the 120 different frequency coverage regions, each represented by a different colour (the progression of colour implies no particular ordering), and none of which contains all 29 frequencies.

Table 1.

List of sky maps we use in our multi-frequency modelling. F: full sky, S: southern sky, N: northern sky, E: equatorial plane, P: partial sky. CHIPASS at 1.39 GHz has a bandwidth of 64 MHz, so its frequency largely overlaps with the Stokert+Villa Elisa map at 1.42 GHz.

Project/Instrument	ν(GHz)	Area	Resolution	Reference(s)
DRAO, CAN	0.01	N	2.6 × 1| $_{.}^{\circ}$ |9	Caswell (1976)
DRAO, CAN	0.022	N	1.1 × 1| $_{.}^{\circ}$ |7	Roger et al. (1999)
MRAO+JMUAR	0.045	F	∼3| $_{.}^{\circ}$ |6	Alvarez et al. (1997); Maeda et al. (1999)
Parkes	0.085	E	3.8 × 3| $_{.}^{\circ}$ |5	Landecker & Wielebinski (1970)
Parkes	0.15	E	2| $_{.}^{\circ}$ |2	Landecker & Wielebinski (1970)
GER, AUS, ENG	0.408	F	56 arcmin	Haslam et al. (1981, 1982); Remazeilles et al. (2015)
Dwingeloo, NLD	0.82	N	1| $_{.}^{\circ}$ |2	Berkhuijsen (1971)
CHIPASS	1.39	S	14.4 arcmin	Calabretta, Staveley-Smith & Barnes (2014)
Stokert, Villa Elisa	1.42	F	36 arcmin	Reich (1982); Reich & Reich (1986), Reich, Testori & Reich (2001)
Rhodes/HartRAO	2.33	S	20 arcmin	Jonas, Baart & Nicolson (1998)
WMAP	22.8	F	49 arcmin	Hinshaw et al. (2009)
Planck	28.4	F	32 arcmin	Planck Collaboration I (2016a)
WMAP	33.0	F	37 arcmin	Hinshaw et al. (2009)
WMAP	40.6	F	29 arcmin	Hinshaw et al. (2009)
Planck	44.1	F	24 arcmin	Planck Collaboration I (2016a)
WMAP	60.8	F	20 arcmin	Hinshaw et al. (2009)
Planck	70.4	F	14 arcmin	Planck Collaboration I (2016a)
WMAP	93.5	F	13 arcmin	Hinshaw et al. (2009)
Planck	100	F	10 arcmin	Planck Collaboration I (2016a)
Planck	143	F	7 arcmin	Planck Collaboration I (2016a)
Planck	217	F	5 arcmin	Planck Collaboration I (2016a)
Planck	353	F	5 arcmin	Planck Collaboration I (2016a)
Planck	545	F	5 arcmin	Planck Collaboration I (2016a)
Planck	857	F	5 arcmin	Planck Collaboration I (2016a)
AKARI	1875	P	1.5 arcmin	Doi et al. (2015)
AKARI	2143	P	1.5 arcmin	Doi et al. (2015)
IRAS (IRIS)	3000	F	4.3 arcmin	Miville-Deschênes & Lagache (2005)
AKARI	3333	P	1.5 arcmin	Doi et al. (2015)
IRAS (IRIS)	5000	P	4 arcmin	Miville-Deschênes & Lagache (2005)

Due to the large frequency range covered by these sky maps and the variety of instruments involved in making them, there is considerable variation in their characteristics. For example, some instruments have variable beam sizes, while others have significant sidelobes. For this work, we approximate all beams as Gaussian, whose full width at half-maximum size (θ_FWHM) equals the resolution quoted in Table 1. We smooth all maps to 5° by applying Gaussian smoothing kernels of sizes | $\sqrt{{5^\circ }^2-\theta _{\rm FWHM}^2}$ |⁠, and remove another 3° from edge areas in incomplete maps. All map are pixelized on to a HEALPix grid with n_side = 64.

Additionally, each instrument faces unique calibration challenges that causes them to differ in error properties. While it is possible to quantify various error properties in the form of error covariance matrices and insert them into equations (6) and (7), we leave such effort to future work.

In terms of unit, the sky maps come in three different units. The maps below 20 GHz are in Rayleigh–Jeans temperature, K_RJ, below 500 GHz and above 20 GHz in CMB temperature, K_CMB, and above 500 GHz in MJy sr⁻¹. We choose to convert all maps into to MJy sr⁻¹ using the following conversion equations:

\begin{equation} \frac{1 {\rm K}_{\rm RJ}}{1 {\rm MJy\, sr^{-1}}} = 2 \times 10^{20} \frac {k_{\rm B}}{\lambda ^2}, \end{equation}

(9)

\begin{equation} \frac{1 {\rm K}_{\rm CMB}}{1 {\rm K}_{\rm RJ}} = \left(\frac{h\nu }{k_{\rm B}T_{\rm CMB}}\right)^2\frac{e^{\frac{h\nu }{k_{\rm B}T_{\rm CMB}}}}{\left( e^{\frac{h\nu }{k_{\rm B}T_{\rm CMB}}}-1\right)^2}, \end{equation}

(10)

where ν is the frequency at the centre of the frequency band of the map, λ = c/ν is the wavelength, k_B is Boltzmann's constant, h is Planck's constant, and T_CMB is the CMB temperature. The root mean square (rms) amplitude of all the maps after unit conversion are shown in Fig. 2.

Figure 2.

The overall amplitudes of the 29 sky maps included in this work. The sizes of the circles represent the sky coverage of each map. The amplitudes can be thought of as the rms value of the input maps, with some iterative adjustments explained in Section 2.3.

4 RESULTS I: THE IMPROVED GSM

4.1 Orthogonal components result

We apply the algorithms described in Section 2 on maps described in Section 3. To obtain the initial 5 per cent common coverage, maps at 10 MHz, 85 MHz, 150 MHz, and 5 THz are temporarily excluded. In the iteration process, it takes 12 iterations with step size η = 1 to converge. We also tried a smaller η = 0.2, and obtained the same result with four times more iterations. After the first convergence, we mask out the 1 per cent of the pixels with the highest fitting residual in order to eliminate point sources, and reiterate for another eight iterations to arrive at the final convergence. The new GSM with six components is shown in Fig. 3. We choose the number of components to be six based on the predictive accuracy of the resulting GSM, which we discuss in more detail in Section 4.2. We find that the increase in predictive accuracy from using more than six components is negligible compared to the fitting errors. We plan to perform a more rigorous model selection to determine the optimal number of components in a future work that will incorporate uncertainty information from each input map.

Figure 3.

Top: the six orthogonal components. Bottom: their spectra. The maps are plotted on arcsinh scales where the constant in arcsinh is set to the median of pixel amplitudes. The percentage on each component represents the fraction of all variations among the 29 input maps explained by that component. Here, the spectra are not directly | $\bf{S}$ |⁠, but the rows of | $\bf{S}$ | multiplied by the overall amplitudes of the input maps shown in Fig. 2. Note that the spectra here have no physical meanings attached to them, which we explore in Section 5.

4.2 Error analysis

To assess the predictive power of the GSM, we compute three types of rms errors at each frequency. The first is the fitting residual that is the difference between the complete GSM model and the input sky maps. This fitting error is typically 1 per cent for CMB frequencies between 20 and 100 GHz, and around 5 per cent everywhere else, as shown in blue in Fig. 4. This is a combination of the measurement errors in the input maps and weak emission mechanisms comparable to such errors, which we are unable to capture. The fitting error represent the lower bound of the predictive accuracy of our GSM.

Figure 4.

Three different rms error percentage estimations for our GSM, compared to the error quoted with the original GSM. Agnostic error is an overly pessimistic upper bound to the predictive accuracy of the GSM, whereas the map agnostic error is the lower bound. Fitting error is the fitting residual in each map not explained by the GSM. The sizes of the circles represent the sky coverage of each map. The improved GSM has less error compared to the original GSM across the entire frequency range.

The second type of error is the map ‘agnostic error’, where at a given frequency, we use the same normalized spectra but do not use the input map at that frequency to calculate the principal component maps. We calculate the rms of the difference between the resulting model and the actual sky map. As shown in green in Fig. 4, this error is typically 1.5–2 times larger than the fitting error. Among the three errors discussed in this section, this error is the most probable estimate of the predictive accuracy of the GSM, and is also used in the original GSM work.

The third type of error is a conservative upper-bound to the predictive accuracy of the GSM. Here, we pretend to have zero knowledge of each map and exclude it from the algorithm from the outset, and produce a new GSM that is agnostic of that map. We then compare the excluded map with the predicted map using the new agnostic GSM, and calculate the rms of the difference after an overall renormalization. The renormalization is typically less than 15 per cent. This is an overly pessimistic estimate of the predictive power of the GSM, because we are not in the regime of an overabundance of sky maps, especially at frequencies below 20 GHz. The only two full-sky maps below CMB (at 408 MHz and 1.4 GHz) serve as ‘cornerstones’ of the GSM, so the agnostic error caused by removing them completely (and thus crippling the GSM) is not a good estimate of the predictive power of the complete GSM at an unexplored frequency such as 300 MHz.

The peaks in the agnostic error curves that rise above 10 per cent fall into three categories. The first category includes those near spectral lines, where the high errors are caused by having a spike at one end of the interpolation. One such example is the strongest peak at 2.3 GHz, whose error mostly comes from interpolating between the peak at 1.4 GHz and the lowest WMAP frequency. Another example is the peak around 100 GHz, which is contaminated by the 115 GHz CO line. The second category includes those at the edge of our frequency range, namely 10 MHz and 5 THz. They have high errors due to extrapolation. The third category includes those that interpolate between incomplete and typically non-overlapping maps. One example is the 408 MHz Haslam map, which is interpolated between the 150 MHz equatorial map and the 820 MHz northern sky map. Another example is the 3 THz IRIS map that is interpolated between two AKARI maps with poor equatorial coverage.

In summary, for a map predicted by the GSM, we expect the rms accuracy to be above the fitting error, and below the completely agnostic error. This is typically between 5 per cent and 15 per cent for most frequencies, and for frequencies near 50 GHz and 500 GHz, the accuracy can be as good as 2 per cent. In addition, we expect an overall amplitude offset of less than 15 per cent. In addition to a wider and denser frequency coverage, the improved GSM is seen to be consistently more accurate, by up to a factor of 2, than the original GSM across the entire frequency range.

5 RESULTS II: RECOMBINE THE GSM TO EXTRACT PHYSICAL COMPONENTS

5.1 Blind component separation

Since the component maps we obtain (Fig. 3) are mutually orthogonal by construction, they are linear combinations of underlying emission mechanisms such as synchrotron and thermal dust. As discussed in Section 2.1, we have the freedom to apply any 6 × 6 invertible matrix Ψ and its inverse to the matrices | $\bf{S}$ | and | $\bf{M}$ | without changing our model, | $\bf{M}\bf{S}$ |⁠. Applying Ψ will not change any linear interpolation results either. In this section, we discuss our two-step automatic algorithm to determine Ψ and blindly extract the physical contributions. This algorithm focuses on manipulating the normalized spectrum matrix | $\bf{S}$ |⁠. The only two pieces of information that we use in order to determine Ψ are: (1) the absence of monopole in CMB and (2) compact frequency support for each normalized spectrum, meaning that the each component dominates a limited frequency range. Our algorithm is thus blind to any existing models for known emission mechanisms.

We start with a Ψ that equals the identity matrix. The first step tries to search for the smallest frequency range outside which there are only five components, meaning that the sixth component must be limited to said frequency range. To do this, we enumerate all possible frequency ranges, and for each frequency range, we remove the columns in | $\bf{S}$ | that correspond to that range to form | $\bf{S}^{\ast }$ |⁠. We then calculate the eigenvalues and eigenvectors for | $\bf{S}^{\ast }\bf{S}^{{\ast }t}$ |⁠. For the smallest frequency range whose smallest eigenmode explains less than 1 per cent of all variations in | $\bf{S}^{\ast }$ |⁠,

\begin{equation} \frac{\lambda _{\rm min}}{{\rm tr}(\bf{S}^{\ast }\bf{S}^{{\ast }t})} < 1{\rm per cent}, \end{equation}

(11)

we multiply Ψ by its eigenvectors | $\bf{C}^{\ast }$ |⁠. This way, the last row of | $\mathbf {\Psi }\bf{S}$ |⁠, which corresponds to the smallest eigenvalue, is our first separated component. If there is no frequency range that satisfies this criterion, we keep doubling the 1 per cent threshold until one valid range appears. In practice, the synchrotron component is separated first with 2 per cent threshold. We then repeat the above procedure, and separate the CMB component at 4 per cent. Lastly, the ‘HI’ component is separated at 8 per cent. The first step ends since we cap the threshold at 10 per cent.

Before we proceed to the second step, we attempt to clean the foreground contamination in the CMB component found in the first step. To do this we find the Ψ that minimizes the power in the CMB component map. This is the only procedure in the entire algorithm where we are working with the component map rather than normalized spectra, also the only procedure that uses physical knowledge, namely the lack of monopole in the CMB.

The automatic procedure mentioned above is not able to separate the last three modes, which suggests that they cover largely overlapping frequency ranges. For our second step, we demand that the last three components each dominate one of the three frequency ranges: 10 GHz to 100 GHz, 100 GHz to 1 THz, and above 1 THz. To do this, we find the Ψ that minimizes the normalized spectrum of each component outside its designated range. The results we obtain are shown in filled dots in Fig. 5.

Figure 5.

Top: six component maps from blind component separation. Bottom: six component spectra (filled points) and their best fits (dashed lines). For plotting purpose, each component map is normalized to have a median of 1 and plotted on arcsinh scales, where the colour range corresponds to 2 and 98 percentiles in each map. In the bottom panel, we only plot the two most dominant component spectra at each frequency to reduce visual clutter. The best-fitting parameters are listed in Table 2.

5.2 Fitting the blind components

As shown in Fig. 5, five of the six components obtained in our blind separation have spectra that closely resemble those of known physical processes: synchrotron, free–free, CMB, and two thermal dust species at different temperatures. In this section, we take the known spectrum models of these mechanisms, such as a power law for synchrotron, and find the best parameters for them that fit our blind spectra. We first compute six spectra by taking the new normalized spectrum matrix | $\mathbf {\Psi }\bf{S}$ | and multiplying back the normalization (shown in Fig. 2). We then fit each spectrum using its corresponding physical model: synchrotron and free–free as power laws with spectral indices β, thermal dust as blackbody spectra multiplied by power laws with temperatures T and spectral indices β, and CMB anisotropy as the first-order Taylor expansion of the blackbody spectrum with temperature T. Note that we do not change our blind spectra based on these models.

The fits are performed in log(frequency)–log(surface brightness) space, and input variances are empirically estimated using the rms of best-fitting residuals. We limit the fitting range for each component to the frequency range it dominates, as shown in Fig. 5. For power-law fits, we perform a standard linear regression to obtain the best fits and error bars. For the dust and CMB components, we first numerically search for the best fit, and then perform Monte Carlo on the input spectra using the empirical variances to estimate the error bars. Although our component separation is blind to any of these model spectra, the best-fitting parameters for all five of our blind spectra are seen to agree very well with existing literature values to within 2σ, as shown in Table 2.

Table 2.

List of model parameters for our blind components compared to previously published values. All parameters that fit our blind components agree with previously published values to within 2σ. The error bar from Platania et al. (2003) represents spatial variation rather than statistical uncertainty. The error bars for previously published dust model parameters are estimated using the difference in values between Finkbeiner, Davis & Schlegel (1999) and Meisner & Finkbeiner (2015).

Component	Parameter	Best fit for blind spectra	Previous value	Ref.
Synchrotron	spectral index β below Haslam	−2.519 ± 0.018	−2.5 ± 0.1	Rogers & Bowman (2008)
	spectral index β above Haslam	−2.715 ± 0.082	−2.695 ± 0.120	Lawson et al. (1987)
				Reich & Reich (1988)
				Davies, Watson & Gutierrez (1996)
				Platania et al. (2003)
Free–free	spectral index β	−2.175 ± 0.032	−2.1 ± 0.03	Dickinson, Davies & Davis (2003)
CMB	temperature T	2.748 ± 0.016 K	2.72548 ± 0.00057 K	Fixsen (2009)
Warm dust	temperature T	20.0 ± 3.3 K	15.95 ± 0.25 K	Finkbeiner et al. (1999),
	spectral index β	2.78 ± 0.70	2.76 ± 0.06	Meisner & Finkbeiner (2015)
Cold dust	temperature T	10.4 ± 1.2 K	9.58 ± 0.18 K	Finkbeiner et al. (1999),
	spectral index β	2.54 ± 0.51	1.65 ± 0.02	Meisner & Finkbeiner (2015)

For future versions of the GSM, a worthy goal would be to rediscover the anomalous microwave emission (AME; Kogut 1996; Leitch et al. 1997; Banday et al. 2003; Lagache 2003; de Oliveira-Costa et al. 2004; Finkbeiner 2004; Davies et al. 2006; Bonaldi et al. 2007; Dobler & Finkbeiner 2008; Miville-Deschênes et al. 2008; Ysard, Miville-Deschênes & Verstraete 2010; Gold et al. 2011) component in addition to the five components mentioned above. Given the currently available data sets, the AME component is likely absorbed into both fitting error and other components such as thermal dusts. The absence of a stand-alone AME component can be attributed to insufficient signal-to-noise ratio, since AME peaks around 20 GHz, which is very close to the 3–20 GHz frequency gap in existing data sets. Future GSM analysis with more data in the relevant frequency range will likely be able to separate the AME component.

5.3 High-resolution GSM

In addition to the 5° resolution GSM we have discussed, we also produce a high-resolution version using the same normalized spectra as the 5° version. Because of successful physical component separation as shown in Fig. 5, we are able to adopt two different resolutions throughout the frequency range. For frequencies above 40 GHz, we can safely neglect the synchrotron component, because it is more than two orders of magnitude below the other four components at most of these frequencies, as shown in Fig. 5. We thus take those WMAP, Planck, AKARI, and IRIS maps whose resolution is better than 24 arcmin, smooth them to 24 arcmin, and solve for the CMB, the free–free, and the two dust component maps. Then, for frequencies below 10 GHz, we take the 408 MHz and 1.42 GHz maps, smooth the latter to 56 arcmin, smooth and remove the four component maps obtained in the high-frequency band, and fit to the synchrotron and ‘HI’ maps. Overall, the high-resolution product has 56 arcmin angular resolution at frequencies below 10 GHz, and 24 arcmin angular resolution above 10 GHz, shown in Fig. 6. Note that while the high-resolution component maps share the same spectra and large-scale features as their low-resolution counterparts, directly smoothing them to 5° will not produce exactly the low-resolution versions, because the high-resolution component maps are calculated using only a subset of the 29 input maps.

Figure 6.

High-resolution version of the six component maps. The five non-CMB maps follow the same colour scales as those in Fig. 5, whereas the CMB component is plotted on a linear scale. The top two components have 56 arcmin angular resolution, and the other four 24 arcmin.

Our blind high-resolution component maps agree remarkably well with existing maps and spatial templates. Our synchrotron, CMB anisotropy, free–free, and cold dust maps share most of the features seen in their counterparts in both the WMAP 9 yr results and Planck 2015 results (see e.g. fig. 19 in Bennett et al. (2013) and fig. 16 in Planck Collaboration I 2016a). Free–free emission is known to closely trace H α emission (Finkbeiner 2003), and our free–free map shows all the key features visible in the composite H α map presented in Finkbeiner (2003). Moreover, our cold dust map share all the key spatial features of the 94 GHz dust map presented in Finkbeiner et al. (1999).

5.4 Comparison with Planck and other foreground models

The blind spectra can be fitted very well in the frequency ranges they dominate, as shown in Fig. 5 and Table 2. On the flip side, however, because our blind separation approach minimizes the normalized spectrum of each component outside the frequency range it dominates, this approach cannot recover very well the parts of the spectrum outside its dominant range. One such example is the free–free spectrum, which is well fitted by a power law with spectral index of −2.175 ± 0.032 near the CMB frequencies, but due to the minimization procedure, it largely disappears in the MHz range where synchrotron dominates. It is still contained in the GSM, but absorbed into the synchrotron and the HI+other components.

A similar problem is evident with our cold dust component. From Table 2, we see a significant difference in our best-fitting spectral index (β = 2.54 ± 0.51) and accepted literature values (e.g. β = 1.65 ± 0.02 from Meisner & Finkbeiner 2015). A visual inspection of Fig. 5 suggests that this is a result of degeneracies with other components. One sees that 100 to 1000 GHz (where cold dust is bright) is a rather ‘crowded’ region with many other components contributing non-negligibly to the sky emission. Moreover, the dominance of cold dust over these other components is weak at best. This makes our cold dust component particularly susceptible to leakage from components such as warm dust, driving our spectral index to a steeper value. We therefore see that our blind component extraction should be interpreted with some degree of caution when applied to frequency ranges where a wide range of foreground emission and mechanisms are at play. To some extent, the resulting degeneracies are captured by the large error bars associated with our dust parameters, but it is important to bear in mind that there may also be systematic errors arising from the intrinsic limitations of our method.

The sixth component, HI+other, is likely a combination of multiple mechanisms and systematic effects. Its strong peak near 1.4 GHz suggests local 21 cm emission as the dominant mechanism. In its component map, we notice two interesting features. There is a blue (low temperature) region whose shape resembles the synchrotron map, so this might indicate another weak synchrotron component with a different spectral index. We also notice a dipole component that is aligned with the equatorial poles, with striping artefacts especially visible near the northern celestial pole, which suggests contribution from scanning systematics in the low-frequency maps.

To further examine our physical foreground components, we compare our data products with those derived from Planck data. The Planck collaboration proposed a global representation of the multi-component sky over nine frequency bands between 30 and 857 GHz. Built upon a Bayesian analysis framework, Planck's foreground maps are created by fitting its observational data to parametric models (Planck Collaboration X 2016b). Since these models require an explicit modelling of foregrounds, they provide a useful contrast to our blind extractions.

To compare the foreground components extracted from the GSM and their counterparts from Planck, we calculate a normalized angular cross-power spectrum | $R_{l}^\mathrm{Planck, GSM}$ |

\begin{equation} R_{\ell }^\mathrm{Planck, GSM}= \frac{C_{\ell }^\mathrm{Planck, GSM}}{\sqrt{C_{\ell }^\mathrm{Planck} C_{\ell }^\mathrm{GSM}}} \end{equation}

(12)

where | $C_{\ell }^{X,Y} = \frac{1}{2l +1}\sum _{m} a_{\ell m}^{X} a_{\ell m}^{Y*}$ | with a_ℓm denoting a spherical harmonic expansion coefficient for spherical harmonic indices ℓ and m. We choose to normalize our cross-powers to minimize effects due to differences in instrumental properties between Planck and the wide variety of data used to construct the GSM. We also choose to split our comparison into regions in and out of the Galactic plane, defined by the map shown in Fig. 7. Note that Planck uses a one-component grey-body thermal dust model, and for our comparison we correlate both dust components with the single dust map from Planck.

Figure 7.

‘In-plane’ (red region) and ‘out-of-plane’ (blue region) masks used for the comparison of our physical foreground components with those from Planck.

Fig. 8 shows the resulting cross-power spectra. The extracted CMB maps are most strongly correlated outside the Galactic plane, where possible confusion with foregrounds (particularly with our blind extraction method) is least likely. The lower correlation within the Galactic plane suggest that there is indeed some leakage from the foreground components. For most of our foreground components, we generally find strong cross-correlations between our maps and those from Planck, with only a very weak dependence on an angular scale. These trends are relatively robust to the precise mask that we use to define the regions inside and outside the plane (for example, one could have chosen a mask that is morphologically more similar to free–free emission than our current mask that is more similar to the patterns of dust emission). There are, however, two minor exceptions to this: we find the out-of-plane curves for free–free and warm dust to shift up and down slightly depending on the exact mask used. This is perhaps not surprising, given that warm dust and free–free are two of the components that are more compactly confined to the Galactic plane, and thus the most sensitive to the definition of the ‘plane’. This does, however, suggest a Galactic latitude dependence on the accuracy of our modelling procedures. A detailed examination of such dependence is beyond the scope of this paper, and investigations are currently underway to incorporate in our modelling the fact that the Galactic plane is a more physically complex region than the Galactic poles. Finally, we find that our sixth component, HI+other, is not highly correlated with any of Planck's foreground products.

Figure 8.

The normalized cross-correlation | $R_{\ell }^\mathrm{Planck, GSM}$ | from equation (12), for the CMB, free–free, synchrotron, and dust components. Overall, there is a reasonable correlation between our extracted foreground components and those from Planck throughout all angular scales on most of components. Red lines denote the cross-correlation inside the Galactic plane, while blue lines show the cross-correlation outside of the Galactic plane. Our CMB map is seen to be only weakly correlated with the corresponding Planck map in the Galactic plane, while highly correlated outside the plane.

Ultimately, it is perhaps unsurprising that there is a reasonable correlation between our extracted components and those from Planck. After all, Planck maps were used to form the GSM in the first place. Thus, the exercise performed here is more a self-consistent demonstration of the ability of our blind method to identify independent foreground components. Our method may thus be helpful for future iterations of the GSM, where it may be used to detect the presence of new foreground components demanded by better data, prior to the availability of a physical model.

6 SUMMARY AND OUTLOOK

We have presented an improved algorithm that builds upon the original PCA algorithm, allowing us to include many more incomplete sky maps across a larger frequency range. We have presented an improved GSM, with expected predictive accuracy between 5 per cent and 15 per cent for most frequencies, and around 2 per cent near CMB frequencies, with overall amplitude offsets less than 15 per cent. It is worth pointing out that these accuracy estimations are based on rms values over the whole sky, so for a particular sky direction, there can be higher errors. We have also presented a blind component separation technique that identifies five physical components that agree very well with existing physical models. Lastly, we also create a high-resolution GSM, with 56 arcmin resolution at frequencies below 10 GHz, and 24 arcmin resolution above 10 GHz.

There are two ways to further improve the GSM in future: adding more data sets and improving the algorithm. On the data side, many high-quality maps will become available in the near future, such as CBASS (King et al. 2010; Irfan et al. 2015), SPASS (Carretti et al. 2009), GMIMS (Wolleben et al. 2009), GEM (Tello et al. 2013), and QUIJOTE (Génova-Santos et al. 2015). Some of these upcoming maps will fill in the gap between 2.3 and 23 GHz. In addition, in a separate paper (Zheng et al. 2016), we present a new imaging method that allows existing low-frequency interferometers such as the MWA to produce high-quality foreground maps that cover nearly the full sky. Last but not least, we aim to produce a polarized GSM using existing and upcoming polarized sky maps.

On the algorithmic side, we would like to include uncertainty information for all sky maps in the form of error covariance matrices of the maps into equations (6) and (7). This would allow us to properly weigh each map rather than resorting to renormalizing all maps and treating them all equally. In addition, this will allow us to perform rigorous model selection in order to decide the optimal number of principal components. Last but not least, we are interested in improving the blind physical component separation procedure to learn more about the physical mechanisms that correspond to the components found in the GSM.

With the help of upcoming high-quality sky maps, the improved GSM will not only serve as a powerful predictive model for diffuse Galactic emission, but will also have the potential to uncover new physical mechanisms that contribute to this emission.

Acknowledgments

This work was supported by NSF grants AST-1105835 and AST-1440343. AL acknowledges support from the University of California Office of the President Multicampus Research Programs and Initiatives through award MR-15-328388, and from NASA through Hubble Fellowship grant HST-HF2-51363.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. The NASA funding for the eGSM is NASA grant no. NNX16AF55G. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. This research has made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We wish to thank Jacqueline Hewitt, Gianni Bernardi, Douglas Finkbeiner, Aaron M. Meisner for helpful comments and suggestions, and Angelica de Oliveira-Costa for pioneering the GSM approach. We also wish to thank our anonymous referee for extremely helpful and detailed feedback.

1

Due to incomplete and different sky coverages between maps, a naive normalization will not weigh the maps properly. We discuss our normalization in more detail in Section 2.3.

2

In practice, no eigenvalues are perfect zeros, and we will discuss our choice of n_c in more details in Section 4.1.

3

MPIfR: http://www3.mpifr-bonn.mpg.de/survey.html, LAMBDA: lambda.gsfc.nasa.gov/product/, Planck: pla.esac.esa.int/pla.