The spectral index of polarized diffuse Galactic emission between 30 and 44 GHz

TL;DR

This paper estimates the polarized spectral index of diffuse Galactic emission between 30 and 44 GHz using Planck data, revealing spatial variations and differences between the northern and southern Fermi bubbles.

Contribution

It introduces a new method combining objective and maximum entropy priors to map the polarized spectral index across the entire sky.

Findings

Average spectral index across the sky is -2.99 with statistical and intrinsic uncertainties.

Detected significant differences in spectral index between the northern and southern Fermi bubbles.

Spectral index varies with location, especially near the Galactic plane and the Fermi bubbles.

Abstract

We present an estimate of the polarized spectral index between the Planck 30 and 44 GHz surveys in $3.7^\circ$ pixels across the entire sky. We use an objective reference prior that maximises the impact of the data on the posterior and multiply this by a maximum entropy prior that includes information from observations in total intensity by assuming a polarization fraction. Our parametrization of the problem allows the reference prior to be easily determined and also provides a natural method of including prior information. The spectral index map is consistent with those found by others between surveys at similar frequencies. Across the entire sky we find an average temperature spectral index of $-2.99\pm0.03(\pm1.12)$ where the first error term is the statistical uncertainty on the mean and the second error term (in parentheses) is the extra intrinsic scatter in the data. We use a clustering algorithm to identify pixels with actual detections of the spectral index. The average spectral index in these pixels is $-3.12\pm0.03(\pm0.64)$ and then when also excluding pixels within $10^\circ$ of the Galactic plane we find $-2.92(\pm0.03)$. We find a statistically significant difference between the average spectral indices in the North and South Fermi bubbles. Only including pixels identified by the clustering algorithm, the average spectral index in the southern bubble is $-3.00\pm0.05(\pm0.35)$, which is similar to the average across the whole sky. In the northern bubble we find a much harder average spectral index of $-2.36\pm0.09(\pm0.63)$. Therefore, if the bubbles are features in microwave polarization they are not symmetric about the Galactic plane.

Tables

Table 1. Table 1: Estimates of the spectral index of polarized diffuse Galactic synchrotron emission at frequencies above ∼ 5 ​ GHz similar-to absent 5 GHz \sim 5\,\mathrm{GHz} .

Frequency	Geometric mean	Spectral index	Notes	Comments	Reference
range [GHz]	frequency [GHz]	$\beta$
All WMAP	-	$-3.2$		All-sky	Kogut et al. (2007)
—"—	-	$-2.98\pm 0.01$	a	Galactic plane	—"—
—"—	-	$-3.12\pm 0.04$	a	High latitude	—"—
—"—	-	$-3.03\pm 0.04$	a	Where	Dunkley et al. (2009a)
—"—	-	$-3.00\pm 0.04$	a		—"—
—"—	-	$-3.08\pm 0.06$	a	$\left|b\right|>{20}^{\circ }$	—"—
All Planck	-	$-3.5$ – $-3.0$	b	All-sky, Commander2	Planck Collaboration et al. (2018b)
—"—	-	$-3.1\pm 0.06$	a	All-sky, SMICA	—"—
1.4–22.8	5.6	$-3.21\pm 0.15$	c	$\left|b\right|>{30}^{\circ }$	Carretti et al. (2010)
2.3–33	8.7	$-3.25\pm 0.15$	c	, $\left|b\right|>{20}^{\circ }$	Krachmalnicoff et al. (2018)
4.78–28.4	11.7	$-3.1658\pm 0.0002$	a	$\delta >-{12}^{\circ }$	Jew (2017)
22.5–32.6	27.1	$-2.99\pm 0.01$	a	All-sky	Fuskeland et al. (2014)
—"—	—"—	$-2.98\pm 0.01$	a	Galactic plane	—"—
—"—	—"—	$-3.12\pm 0.04$	a	High latitude	—"—
23, 33 and 41	31.5	$-3.06\pm 0.02$	a	Average over 18 regions	Vidal et al. (2015)
28.4–44.1	35.4	$-2.99\left(\pm 0.03\right)\pm 1.12$	d	Whole sky	This work (Section 5)
—"—	—"—	$-3.12\left(\pm 0.03\right)\pm 0.64$	d	Clustered points	—"—
—"—	—"—	$-2.92\left(\pm 0.03\right)\pm 0.48$	d	Clustered points, $\left|b\right|>{10}^{\circ }$	—"—

Table 2. Table 2: Summary statistics of the Planck surveys used in this work. The centre frequency is for emission with a temperature spectral index of β = 0.0 𝛽 0.0 \beta=0.0 . The colour corrections are for emission with a temperature spectral index of β = − 3.0 𝛽 3.0 \beta=-3.0 . The noise refers to the noise levels in maps smoothed to 220 arcmin and downgraded to HEALPix N side = 16 subscript 𝑁 side 16 N_{\mathrm{side}}=16 pixels.

	Planck 30 GHz	Planck 44 GHz
Centre frequency [GHz]	28.4	44.1
Colour correction	0.999	0.987
$K_{\mathrm{RJ}}/K_{\mathrm{CMB}}$	0.979	0.951
Median $P$ [$\mu$KRJ]	7.0	2.2
Median ${\sigma }_P$ [$\mu$KRJ]	0.6	0.7

Table 3. Table 3: The average polarized spectral indices and intrinsic scatters that we find between the Planck 30 and 44 GHz surveys in the eighteen regions defined by Vidal et al. ( 2015 ) , along with the spectral indices found by Vidal et al. ( 2015 ) between the WMAP low-frequency polarization surveys. The regions are rectangles in Galactic longitude and latitude, defined by their centre coordinates ( l 0 subscript 𝑙 0 l_{0} , b 0 subscript 𝑏 0 b_{0} ) and their widths.

	$l_0$	$b_0$	$l_{\mathrm{width}}$	$b_{\mathrm{width}}$	$\beta$ from Vidal et al. (2015)			This work
Region	[deg]				$K$–$Ka$	$K$–$Q$	$Ka$–$Q$		${\sigma }_{\mathrm{intrinsic}}$	Notes
1	5.0	75.0	30.0	10.0	$-2.29\pm 0.33$	$-3.13\pm 0.47$	$-4.74\pm 1.44$	$-3.20\pm 0.29$	–	a
2	25.0	57.0	10.0	14.0	$-2.12\pm 0.48$	$-3.49\pm 0.92$	$-7.43\pm 3.34$	$-2.95\pm 0.12$	$0.02\pm 0.06$
3	35.0	26.5	10.0	27.0	$-3.30\pm 0.15$	$-3.15\pm 0.16$	$-2.80\pm 0.50$	$-3.15\pm 0.05$	$0.04\pm 0.07$
4	20.0	27.5	10.0	5.0	$-2.70\pm 0.52$	$-2.26\pm 0.45$	$-1.10\pm 1.40$	$-2.38\pm 0.15$	$0.09\pm 0.23$
5	20.0	18.5	10.0	7.0	$-3.61\pm 0.39$	$-2.61\pm 0.30$	$-0.75\pm 1.06$	$-2.81\pm 0.11$	$0.02\pm 0.08$
6	-6.5	21.5	9.0	7.0	$-2.62\pm 0.44$	$-3.04\pm 0.55$	$-3.03\pm 1.52$	–	–	b
7	0.0	13.5	10.0	7.0	$-3.52\pm 0.51$	$-2.34\pm 0.32$	$0.34\pm 1.27$	$-2.87\pm 0.75$	$1.00\pm 0.97$
8	12.0	7.5	10.0	5.0	$-2.89\pm 0.20$	$-3.05\pm 0.24$	$-3.38\pm 0.73$	–	–	b
9	-6.5	7.5	7.0	5.0	$-3.07\pm 0.20$	$-3.13\pm 0.23$	$-3.10\pm 0.68$	–	–	b
10	-14.5	11.0	7.0	8.0	$-3.22\pm 0.36$	$-3.03\pm 0.36$	$-2.75\pm 1.16$	–	–	b
11	32.0	0.0	24.0	8.0	$-3.08\pm 0.07$	$-3.06\pm 0.07$	$-3.06\pm 0.24$	$-3.1\pm 0.13$	–	a
12	10.0	0.0	10.0	6.0	$-3.07\pm 0.04$	$-3.06\pm 0.04$	$-3.07\pm 0.12$	$-3.35\pm 0.53$	$1.15\pm 0.54$
13	345.0	0.0	20.0	6.0	$-3.03\pm 0.02$	$-3.09\pm 0.03$	$-3.20\pm 0.09$	$-3.25\pm 0.25$	$0.70\pm 0.22$
14	354.0	-11.0	12.0	8.0	$-2.84\pm 0.20$	$-3.03\pm 0.25$	$-3.17\pm 0.76$	$-3.28\pm 0.05$	$0.03\pm 0.06$
15	-2.0	-26.5	16.0	13.0	$-3.05\pm 0.28$	$-4.00\pm 0.52$	$-5.89\pm 1.63$	$-3.20\pm 0.11$	$0.05\pm 0.12$
16	147.5	11.5	15.0	7.0	$-2.63\pm 0.38$	$-2.66\pm 0.41$	$-2.43\pm 1.22$	$-3.15\pm 0.07$	$0.04\pm 0.08$
17	135.0	0.0	30.0	6.0	$-2.64\pm 0.16$	$-2.67\pm 0.16$	$-2.79\pm 0.49$	$-3.57\pm 0.04$	$0.12\pm 0.05$
18	138.5	-19.0	9.0	12.0	$-2.10\pm 0.39$	$-2.37\pm 0.41$	$-3.00\pm 1.29$	$-2.98\pm 0.44$	–	a

Table 4. Table 4: The mean spectral indices in the north and south bubbles and the evidence that they are different. In the first column, when not subtracting an estimate of the CMB from the maps and instead marginalizing over the CMB amplitude. In the second column, when subtracting an 8.9 K monopole from the Haslam map that forms the weakly informative prior. In the remaining columns, after scaling the σ P subscript 𝜎 𝑃 \sigma_{P} noise levels of the Planck maps by factors of 1.0, 2.0 and 10.0. The main results in this paper correspond to the third column, i.e. with a σ P subscript 𝜎 𝑃 \sigma_{P} scaling factor of unity and with the SMICA CMB estimate subtracted from the Planck polarization maps and no Haslam monopole subtracted from the prior.

	No CMB	Remove Haslam	Main result	${\sigma }_P$ scaling factor
	subtraction	monopole	from this work	2.0	10.0
	$-2.47\pm 0.10$	$-2.37\pm 0.09$	$-2.36\pm 0.09$	$-2.52\pm 0.10$	$-3.05\pm 0.23$
	$-3.09\pm 0.05$	$-3.03\pm 0.05$	$-3.00\pm 0.05$	$-3.18\pm 0.06$	$-3.60\pm 0.22$
Bayes factor	$4\times {10}^5$	$8\times {10}^7$	$3\times {10}^7$	$3\times {10}^6$	4
	$-2.82\pm 0.06$	$-2.60\pm 0.08$	$-2.72\pm 0.06$	$-2.88\pm 0.06$	$-3.33\pm 0.16$

Equations

$$p(A_1,A_2|P_1,P_2)\propto p(P_1,P_2|A_1,A_2)p(A_1,A_2)$$

$$\beta (n)=\frac{\log \left(A_1(n)/A_2(n)\right)}{\log ({\nu }_1/{\nu }_2)}.$$

$$p(P_i|A_i,{\sigma }_P)$$

$$=\text{Rician}(P_i,A_i,{\sigma }_P),$$

$$p(P_1,P_2|A_1,A_2)=\prod _{i=1}^2\text{Rician}(P_i,A_i,{\sigma }_{Pi}).$$

$$p(P_i$$

$${\int }_0^{2\pi }\exp \left[-\frac{(P_i\cos \theta -{\mu }_{Qi})}{2{\sigma }_{Qi}}-\frac{(P_i\sin \theta -{\mu }_{Ui})}{2{\sigma }_{Ui}}\right]\mathrm{d}\theta .$$

$$p_{ref.}(A_1,A_2)\propto A_1{A}_2.$$

$$p_{max.ent.}(A_1,A_2)=\prod _{i=1}^2{\lambda }_i\exp (-{\lambda }_i{A}_i).$$

$$p(A_1,A_2)=\prod _{i=1}^2{A}_i{\lambda }_i^2\exp (-{\lambda }_i{A}_i),$$

$$p(\beta )=-\frac{6{\lambda }_1^2{\lambda }_2^2\exp (2\beta /r)}{r({\lambda }_1\exp (\beta /r)+{\lambda }_2{)}^4},$$

$$D_{\mathrm{KL}}(p_{\mathrm{posterior}}||p_{\mathrm{prior}})={\int }_{-\infty }^{\infty }{p}_{\mathrm{posterior}}(\beta )\log \left(\frac{p_{\mathrm{posterior}}(\beta )}{p_{\mathrm{prior}}(\beta )}\right)\mathrm{d}\beta .$$

$$y(\theta ,\varphi )=h\ast x(\theta ,\varphi )+n(\theta ,\varphi ),$$

$$\hat{x}(\theta ,\varphi )=(g\ast y)(\theta ,\varphi ),$$

$$G(\ell )=\frac{1}{H(\ell )}\frac{H(\ell {)}^2}{H(\ell {)}^2+N(\ell {)}^2/X(\ell {)}^2},$$

$$B(\ell )=\exp \left(-\frac{1}{2}\ell (\ell +1){\sigma }^2\right),$$

$$p(A_{\mathrm{CMB}}^{\prime })$$

$$=\frac{A_{\mathrm{CMB}}^{\prime }}{{{\sigma }_{\mathrm{CMB}}^{\prime }}^2}{\mathrm{e}}^{-{A_{\mathrm{CMB}}^{\prime }}^2/2{{\sigma }_{\mathrm{CMB}}^{\prime }}^2},$$

Full text

††pubyear: 2015††pagerange: The spectral index of polarized diffuse Galactic emission between 30 and 44 GHz–The spectral index of polarized diffuse Galactic emission between 30 and 44 GHz

1 Introduction

Diffuse Galactic synchrotron radiation is emitted by cosmic ray electrons and positrons spiralling through the Galactic magnetic field (Rybicki & Lightman, 1985). In polarization, this radiation is the dominant foreground to the cosmic microwave background (CMB) below around 100 GHz (Page et al., 2007). Unlike the already detected CMB temperature fluctuations or CMB $E$-mode signal, the primordial $B$-mode signal is possibly weaker than the polarized synchrotron foreground at all frequencies, angular scales and areas of the sky (Dunkley et al., 2009b). Better constraints on the frequency spectrum of diffuse Galactic synchrotron radiation are therefore vital if the primordial $B$-mode signal is to be unambiguously detected.

The frequency spectrum of synchrotron radiation is determined by the energy distribution of the radiating cosmic rays, which is approximately a power law above a few GeV for around two decades in energy (Casadei & Bindi, 2004; Aguilar et al., 2014). Diffuse Galactic synchrotron radiation therefore has a power-law frequency spectrum over approximately four decades in frequency.

Except in regions with continuous injection of new relativistic cosmic-ray electrons, the synchrotron spectrum has a cut-off at high frequencies caused by radiative energy loss. At low-frequency ($\lesssim 100\,\mathrm{MHz}$) the spectrum turns over due to self absorption (Condon & Ransom, 2016). Deviations from a power law can be introduced by localised variations in the distribution of cosmic-ray energies, by line-of-sight effects (Chluba et al., 2017), and in polarization by Faraday rotation (for example Lazarian & Pogosyan, 2016).

For diffuse Galactic synchrotron radiation, observations from tens of megahertz to hundreds of gigahertz find an average total-intensity temperature spectral index between 2.5-3.5 (for example Lawson et al., 1987; Reich & Reich, 1988; Platania et al., 2003; Davies et al., 2006; Guzmán et al., 2011; Planck Collaboration et al., 2015c, d). These measurements agree with simulations such as Orlando & Strong (2013).

Synchrotron radiation can be highly polarized, in principle up to 75 per cent (Rybicki & Lightman, 1985). At high Galactic latitude diffuse Galactic synchrotron is measured to be up to $\sim 40$ per cent polarized but is more typically 20 per cent polarized with much lower polarization fractions close to the Galactic plane (for example Kogut et al., 2007; Vidal et al., 2015; Planck Collaboration et al., 2015d).

A summary of measurements of the average spectral index of polarized diffuse Galactic synchrotron radiation is presented in Table 1. The tabulated results were obtained using a variety of methods and datasets. The methods include template fitting, $T-T$ plots and parametric fitting. The datasets include those from space-based CMB satellites WMAP (Bennett et al., 2013; Jarosik et al., 2011; Hinshaw et al., 2009; Bennett et al., 2003) and Planck (Planck Collaboration et al., 2018a, 2015a), and ground-based surveys like S-PASS (Carretti et al., 2019), C-BASS (Jones et al., 2018) and PGMS (Carretti et al., 2010). Many of the analyses listed in Table 1 assumed a single spectral index across the sky.

Variations in the spectral index typically trace structures in our Galaxy. At radio frequencies the most prominent structures are the radio loops and spurs, the most pronounced of which is the North Polar Spur (Loop I) (Roger et al., 1999). These loops are most likely caused either by the shock-waves from supernova explosions accelerating electrons in the inter-stellar medium or stellar winds from OB associations (Berkhuijsen et al., 1971; Heiles et al., 1980; Salter, 1983; Wolleben, 2007; Mertsch & Sarkar, 2013; Vidal et al., 2015). In the shock waves of recent supernova explosions the high energy electrons have not had time to radiate away their energy and this can result in hard (more shallow) spectra, older populations of cosmic-ray leptons emit softer (steeper) synchrotron emission.

One particularly interesting feature is the WMAP Haze. The Haze is a diffuse structure with a total intensity spectral index of about $\beta\simeq-2.5$ that is centred on the Galactic Centre (Finkbeiner, 2004b; Dobler & Finkbeiner, 2008; Dobler, 2012a; Planck Collaboration et al., 2012). The Haze is not symmetric about the Galactic plane, it is weaker in the south (Su et al., 2010), and this north/south asymmetry could be caused by differences in the circumgalactic medium (Sarkar, 2019). It is important to note that the Haze has many overlapping features including the Gould Belt (particularly north of the plane), and it spans a large region covering many tens of kiloparsecs including the Galactic centre.

On the sky, the extent of the WMAP Haze approximately coincides with that of the Fermi Bubbles, a double-lobed feature observed in Fermi-LAT data that extends to Galactic latitudes of $b=\pm 55^{\circ}$ (Dobler et al., 2010; Su et al., 2010; Dobler, 2012b). Unlike the Haze, the Fermi Bubbles are approximately symmetric about the Galactic plane and have well defined edges. A number of mechanisms for producing both the Fermi Bubbles and the WMAP Haze have been proposed including annihilation of Dark Matter particles (Finkbeiner, 2004a) and AGN-like activity at the Galactic centre (Su & Finkbeiner, 2012; Yang et al., 2013).

Detection of the haze in polarization has proven more difficult, Gold et al. (2011) found no indication for a haze of polarized hard spectrum synchrotron radiation. Polarized outflows from the Galactic centre have been observed that spatially correlate with the bubbles/haze (Jones et al., 2012; Carretti et al., 2013) and could be related (Crocker et al., 2015).

In this work we use Bayesian methods to estimate the spectral index from maps of polarized intensity. Our parametrization of the problem allows the objective reference prior for the problem to be easily determined and information from total intensity observations to be easily included. Much of the publicly available low-frequency polarization surveys (here defined to be below 50 GHz) are either low signal-to-noise or corrupted by Faraday rotation. We use the Planck 30  and 44 GHz surveys as these are two of the highest signal-to-noise polarized all-sky maps publicly available, that are also not so low in frequency for Faraday effects to be significant (except along the Galactic plane). We do not use the WMAP polarization maps as they have poorly constrained large-scale modes (Jarosik et al., 2011). These modes effectively introduce varying zero-levels across the $Q$ and $U$ maps that corrupt the polarized intensities.

We make a careful choice of prior that maximises the impact of the data on the posterior and that also includes weakly informative prior information from total intensity data. We use the relative entropy between the posterior and the prior along with a clustering algorithm to identify pixels with good detections of the spectral index.

The paper is laid out as follows. In Section 2 we introduce our method and justify our choice of prior and posterior distributions. In Section 3 we explain the processing that we applied to the data used in this work. In Section 4 we demonstrate that the method works well on simulated data. In Section 5 we present our results using the method on real data and show that our results are generally consistent with the work of others. We find that the North and South Fermi bubbles have different average spectral indices. The spectral index in the northern bubble is shallower than in the southern bubble and we discuss possible causes for this observation. We conclude in Section 6.

2 Model and method

We want to estimate the spectral index from two noisy measurements of the polarized intensity at different frequencies, and we also want to incorporate our prior knowledge into the analysis. Instead of working with maps of the Stokes $Q$ and $U$ parameters, we use the polarized intensities. The polarized intensity is insensitive to Faraday rotation (although not depolarization) and also allows us to more easily encode our prior knowledge from total intensity data. Calculating the analytic form of the posterior distribution of the spectral index is beyond the scope of this work and so instead we estimate it numerically.

We do this by sampling from the posterior probability distribution of the polarized intensities, which is constructed by the multiplication of a likelihood function by a prior (and a normalizing constant, Bayes & Price, 1763),

[TABLE]

where $A_{i}$ and $P_{i}$ are the true and measured polarized intensities at frequency $\nu_{i}$ respectively.

The spectral index between polarized emission at two different frequencies, $\nu_{1}$ and $\nu_{2}$, in the $n^{\textrm{th}}$ pixel of a map is then

[TABLE]

In Section 2.1 we describe how we approximate the polarized intensities as being Rician random variables. In Section 2.2 we describe how our prior is formed by multiplying the objective reference prior by the maximum entropy distribution for the testable information that we wish to include. In Section 2.3 we explain how we explore the resulting posterior distribution for the spectral index with MCMC methods and how we compare the posteriors to the prior.

2.1 Likelihood function

The polarized intensity, $P$, is formed from the quadrature sum of the measured Stokes $Q$ and $U$ values, which can themselves be modelled as Gaussian random variables. If the errors on $Q$ and $U$ are equal then the measured values of $P$ are drawn from a Rician distribution

[TABLE]

where $I_{0}$ is the zeroth order modified Bessel function of the first kind and $\sigma_{P}$ is the standard deviation of noise on $Q$ and $U$. The likelihood function for the true intensities in a pixel $n$ in maps at frequencies $\nu_{1}$ and $\nu_{2}$ is then the product of two Ricians,

[TABLE]

If the errors on $Q$ and $U$ are not equal then $P$ is distributed according to a Beckmann distribution (also known as a generalized Rician distribution, Beckmann, 1959; Yongjun Xie & Yuguang Fang, 2000; Zhu et al., 2018),

[TABLE]

The posterior distributions of the true polarized intensities, and therefore spectral indices, could then be determined from the posterior distributions of $\mu_{Qi}$ and $\mu_{Ui}$. There are two main reasons that we do not adopt this likelihood function or parametrization. Firstly, the integral in Equation 2.1 can not be solved analytically and so this likelihood function is expensive to compute. Secondly, to construct a weakly informative prior in this instance becomes difficult. The prior must allow for negative values of $\mu_{Q}$ and $\mu_{U}$ but there is no maximum entropy distribution for a variable with known expectation value that extends over the entire real line (Section 2.2) and furthermore it requires prior information on the expected polarization angles across the sky.

2.2 Prior

Away from the Galactic plane and brightest radio loops, the Planck polarization data are low signal-to-noise and as such the prior is important. We want as non-informative a prior on the true polarized intensities as possible, whilst still including weakly informative constraints from other observations. Our prior is formed by multiplying the objective reference prior for the problem by the maximum entropy prior for the information we want to include.

The reference prior is the prior distribution that maximizes the Kullback-Leibler divergence (or relative entropy, Kullback & Leibler, 1951) between itself and the posterior distribution, i.e. it allows the data to have maximal impact on the posterior distribution (Bernardo, 1979, 1981; Berger & Bernardo, 1992; Bernardo, 2005). Reference priors have several properties in their favour over other proposed non-informative priors; the process to generate them is completely general, they are invariant to parameter transformations, and they have consistent marginalization and sampling properties. They have been used in cosmology by Heavens & Sellentin (2018) to ensure that prior assumptions are not dominating estimates of the neutrino mass hierarchy.

Because the likelihood function in Equation 4 can be neatly factorized into the form $p(\vec{P}|\vec{A})=\prod_{i}p(P_{i}|A_{i})$, the reference prior is simply the product of the one dimensional Jeffreys priors for each term in the product. The Jeffreys prior for the Rician disitrbution is $p(A_{i})\propto A_{i}$, which corresponds to a uniform distribution in the $Q-U$ plane (Jeffreys, 1946; Lauwers et al., 2009). The reference prior for this problem is therefore

[TABLE]

We are not totally ignorant of the expected polarized intensity of diffuse Galactic emission from 1–100 GHz and therefore should include this information in our prior. As stated earlier: synchrotron radiation is typically $\sim 20\,\textrm{per cent}$ polarized, its frequency spectrum approximately obeys a power law, and low-frequency total intensity surveys are dominated by synchrotron radiation. To calculate our prior expected polarized intensity as a function of frequency we extrapolate the 408 MHz Haslam survey (Haslam et al. (1982) as reprocessed by Remazeilles et al. (2015)) to higher frequencies with temperature spectral index of -3.11 and assume a uniform polarization fraction of 0.2 across the sky. The prior expected polarized intensity could be improved in future work using a spatially varying total intensity spectral index and polarization fraction (for exmaple those produced by Miville-Deschênes et al., 2008, note that strictly the prior should not be formed from data included in the likelihood).

We do not have strong prior knowledge on the expected spread of the polarized intensities from this value. We therefore multiply the reference prior by the maximum entropy prior for a positive valued variable with known expectation value; i.e. an exponential distribution,

[TABLE]

The total prior on parameters $A_{1}$ and $A_{2}$ is therefore the product of Equations 6 and 7,

[TABLE]

where $1/\lambda_{i}$ is the prior expected value of $A_{i}$ and it is raised to the second power to ensure the prior integrates to unity.

Both the likelihood and prior for $A_{i}$ are unimodal distributions with positive skew. The prior is only more peaked than the likelihood when $1/\lambda_{i}\lesssim\sigma_{Pi}$. When this condition is met then the peak of the posterior is shifted from that of the likelihood to much closer to the prior. When it is not met the effect of the prior on the posterior is small.

From the prior distribution in Equation 8, we can analytically derive the prior distribution on the spectral index $\beta$,

[TABLE]

where $r=1/\log{(\nu_{1}/\nu_{2})}$.

2.3 Posterior distribution of the spectral index parameter

A calculation of the analytic form of the posterior distribution for $\beta$ is beyond the scope of this work. Instead we sample the posterior using Markov Chain Monte Carlo (MCMC) methods.

The posterior distribution for $A_{1}$ and $A_{2}$ is the product of Equations 4 and 8. We explore the posterior distribution for $A_{1}$ and $A_{2}$ in each pixel using the Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970), implemented using PyMC (Patil et al., 2010). For each step in the converged chain we calculate the spectral index using Equation 2, these are then samples drawn from the spectral index posterior distribution.

We can quantify the increase in information on the spectral index in each pixel by calculating the Kullback-Liebler divergence, $D_{\mathrm{KL}}$, between the analytic equation for the prior (Equation 9) and our numerical estimate of the posterior distribution (Kullback & Leibler, 1951),

[TABLE]

3 Data

In this section we describe the Planck 30 and 44 GHz polarized intensity maps (Planck Collaboration et al., 2018c) that we estimate the spectral index from along with the ancillary Haslam total intensity survey, which we use to form a weakly informative prior (Section 2.2), and the SMICA CMB map (Planck Collaboration et al., 2018b), which we use to subtract an estimate of the CMB from the Planck maps. The Planck maps used in this analysis are summarised in Table 2. In this work all maps use the HEALPix pixelization scheme (Gorski et al., 2005).

In each pixel, if the noise on $Q$ and $U$ is Gaussian with $\sigma_{Q}=\sigma_{U}$ then the polarized intensity, $P=\sqrt{Q^{2}+U^{2}}$, is distributed like a Rician random variable. In the Planck 30 GHz map the Median ratio of $\sigma_{Q}$ to $\sigma_{U}$ is $\mathrm{Median}(\sigma_{Q}/\sigma_{U})=1.00$. In the Planck 44 GHz map $\mathrm{Median}(\sigma_{Q}/\sigma_{U})=1.07$ and therefore $\sigma_{Q}\not\approx\sigma_{U}$. To ensure that $P$ is Rician distributed, we add small amounts of extra noise to the $Q$ and $U$ maps. In each pixel, if $\sigma_{Q}>\sigma_{U}$ then we add a number drawn from a Gaussian distribution with zero mean and variance $\sigma_{Q}^{2}-\sigma_{U}^{2}$, and vice versa for pixels with $\sigma_{U}>\sigma_{Q}$. The $\sigma_{P}$ map is then the maximum of $\sigma_{Q}$ and $\sigma_{U}$ in each pixel.

Joint analysis of multiple datasets is easiest if all maps are at a common resolution. Therefore, after adding small noise realizations and subtracting the SMICA estimate of the CMB from the $Q$ and $U$ maps, we deconvolve the beams from the maps and smooth to a common resolution of $3.7^{\circ}$ by reweighting the maps in $a_{\ell m}$ space. $3.7^{\circ}$ is approximately the width of an $N_{\mathrm{side}}=16$ HEALPix pixel, which we downgrade the resolution of the maps to.111The pixel resolution of a HEALPix map is given by the $N_{\textrm{side}}$ parameter, the sky is broken into $12\times N_{\textrm{side}}^{2}$ equal-area pixels. The smoothed-and-downgraded Haslam, SMICA, and CMB-subtracted, smoothed-and-downgraded Planck maps with their uncertainties are shown in Figure 1.

We approximate the Haslam beam as a Gaussian with $\theta_{\mathrm{FWHM}}=56^{\prime}$. To smooth this map to $3.7^{\circ}$ resolution therefore requires a convolution with a Gaussian with $\theta_{\mathrm{FWHM}}=\sqrt{220^{\prime 2}-56^{\prime 2}}=212.8^{\prime}$. We approximate the beam in the SMICA map as a Gaussian with $\theta_{\mathrm{FWHM}}=10^{\prime}$ and therefore smooth it to our target resolution of $220^{\prime}$ by a convolution with a Gaussian with $\theta_{\mathrm{FWHM}}=219.8^{\prime}$

For the Planck surveys, we model the sky maps as being the true sky signal, convolved by an azimuthally symmetric beam and corrupted by additive white noise,

[TABLE]

where $y$ is the measured signal, $(\theta,\phi)$ are the angular coordinates, $h$ is the beam profile, $x$ is the true un-convolved sky signal and $n$ is Gaussian white noise. We construct an estimator of $x$,

[TABLE]

where $g$ is an azimuthally symmetric function.

The Legendre transform of the function $g$ that provides the optimal least-squares-error estimator of $x$ is

[TABLE]

where capital letters denote harmonic transforms.

In the limit of infinite signal-to-noise $G$ tends towards the simple inverse filter, $1/H$. For $H$, we use the beam profiles in Planck Collaboration et al. (2015b). Both $X$ and $N$ are not known a priori but must be estimated from the data themselves. We approximate $X$ by fitting a power law to the angular power spectrum of the map between $\ell=10\--100$ and approximate $N^{2}$ as white noise from the high-$\ell$ part of the spectrum.

The deconvolution and smoothing operations are commutative and independent and so after deconvolution we smoothed the maps to a Gaussian using the transfer function

[TABLE]

where $\sigma$ is the standard deviation of the Gaussian.

We downgrade the maps to lower $N_{\mathrm{side}}$ in $E-B$ space using the HEALPix ud_grade subroutine before transforming back to $Q-U$.222The ud_grade subroutine does not include parallel transport.

After smoothing and downgrading, the variance maps need to be scaled by a correction factor. This factor is difficult to calculate analytically and so instead we estimate it using simulations.

From the smoothed and downgraded maps, in each pixel we construct $P(n)=\sqrt{Q(n)^{2}+U(n)^{2}}$ and approximate $\sigma_{P}(n)=\mathrm{max}(\sigma_{Q}(n),\sigma_{U}(n))$, where $n$ is the pixel number.

We also generate a simulated dataset on which to test the method. We use the PySM s1 model (Thorne et al., 2017) to simulate the Stokes $I$, $Q$ and $U$ diffuse Galactic synchrotron emission at 28.4 GHz and $3.7^{\circ}$ resolution in an $N_{\textrm{side}}=16$ map. Briefly this method uses the reprocessed, degree-scale smoothed Haslam $408\,\mathrm{MHz}$ map as an $I$ template (Remazeilles et al., 2015), and the 9-year WMAP 23 GHz maps smoothed to $3^{\circ}$ as $Q/U$ templates (Bennett et al., 2013). Small scale features are added to these templates in PySM by generating Gaussian realisations from extrapolations of the angular power spectra to high $\ell$. The synchrotron templates are extrapolated to other frequencies using a power-law frequency spectrum with a spatially varying spectral index taken from Miville-Deschênes et al. (2008) assumed to be the same in total intensity and polarization.

We then extrapolate the $Q$ and $U$ maps from 23 GHz to 44.1 GHz using a spatially constant temperature spectral index of $\beta=-3$. We add noise realisations to all the maps and create $P$ maps from the noisy maps of $Q$ and $U$. The noise realisation for each pixel was generated by multiplying a number drawn from a standard normal distribution by $\sigma_{P}$ in that pixel from the Planck noise maps.

We also need a simulated total intensity map to use as a prior. For this we used the PySM total intensity synchrotron map at 28.4 GHz, corrupted by a 5 per cent multiplicative error in each pixel, again scaled to 44.2 GHz with a temperature spectral index of -3, and apply a polarization fraction of 20 per cent.

4 Verification on simulated data

In this section we test our method on the simulated Planck data. We show in high signal-to-noise pixels the posterior estimates of the spectral index are approximately Gaussian and not significantly biased. In lower signal pixels the method recovers posterior estimates that encode our lack of knowledge about the spectral index. We describe how we use the mean shift clustering algorithm (Comaniciu & Meer, 2002) implemented in scikit-learn (Pedregosa et al., 2012) to group the pixels based on the posterior distributions of their spectral indices. This allows us to identify where the spectral index has been well constrained.

The posterior distributions of the amplitude and spectral index parameters for four pixels that represent high, intermediate, low and very low signal-to-noise signal-to-noise ratios are shown in Figure 2. The Kullback-Liebler divergences between the posterior and prior distributions on the spectral indices ($D_{\mathrm{KL}}$) are 0.66, 0.31, 0.22 and 0.87 in decreasing order of signal-to-noise, this pattern of large Kullback-Liebler divergences in pixels with either strong detections or such low signal-to-noise that the data can not constrain the spectral index is observed across the entire sky. Above each panel we label the mean ($\mu_{\beta}$) and standard deviation ($\sigma_{\beta}$) of the posterior distribution of the spectral index.

Histograms of $\mu_{\beta}$, applying three cuts on $\sigma_{\beta}$, are shown in Figure 3a. All three histograms peak at the true synchrotron spectral index (dashed red line). The weighted average spectral index across the sky is $\left<\beta\right>=-2.966\pm 0.006$. There is a slight excess of pixels with shallower spectral indices caused by pixels with very low signal-to-noise. Given that neither the priors nor the likelihoods are Gaussians, there is no reason for the posteriors to be Gaussians. However, we can still make that approximation and the normalized deviations of the estimates from the true value of -3.0 are shown in Figure 3b. The distribution is close to the standard normal but slightly skewed towards positive deviations (shallower indices).

The skewed distribution of the normalized deviation comes from the lowest signal-to-noise pixels. Maps of the mean ($\mu_{\beta}$), standard deviation ($\sigma_{\beta}$) and Kullback-Liebler divergence from the prior ($D_{\mathrm{KL}}$), of the posterior distribution of the spectral index are shown in Figure 4. In low signal-to-noise regions, the posterior estimates of the spectral index have high uncertainties and in the faintest pixels (top right region of the maps) the means of the posterior distributions are biased towards shallow spectral indices.

2-dimensional histograms of $\mu_{\beta}$, $\sigma_{\beta}$ and $D_{\mathrm{KL}}$ are shown in Figure 5. There are at least two populations of pixels, in the $D_{\mathrm{KL}}/\beta$ histogram there is an arc of pixels whose relative entropy increases with distance from the mean prior value and a second cluster with spectral index close to the true value and larger relative entropy ($D_{\mathrm{KL}}\simeq 0.75$). In the $D_{\mathrm{KL}}/\sigma_{\beta}$ histogram, most pixels lie along the line that is horizontal before a break at $\sigma_{\beta}\sim 0.6$ where the relative entropy begins to decrease with increasing $\sigma_{\beta}$.

We use the mean shift clustering algorithm to cluster the pixels into groups based on where they are in $\mu_{\beta}$, $\sigma_{\beta}$ and $D_{\mathrm{KL}}$ space. The mean shift algorithm starts walkers at the positions of the pixels and the walkers iteratively make steps towards regions of greater density (Comaniciu & Meer, 2002). When the positions of the walkers have converged, those that end up close together are grouped. The tunable parameter in this algorithm is the bandwidth that is used to calculate the local density. Before clustering we scale all variables so that they span a range of 1. We then cluster using a bandwidth of 0.075. The pixels belonging to the largest group (approximately half of the sky) have been identified in red on the 2-dimensional histograms. This cluster of pixels correspond to those where the spectral index has been detected and well constrained by the data.

Instead of simply calculating a weighted average spectral index, by approximating the spectral index posterior distributions as Gaussians and assuming that the true spectral indices are themselves drawn from a Gaussian distribution, we can fit for a mean value and the extra intrinsic scatter in the data that can not be accounted for by the uncertainties. Across the entire sky the mean spectral index is $\left<\beta\right>=-2.963\pm 0.006$ and the extra intrinsic scatter is $\sigma_{\mathrm{intrinsic}}=0.036\pm 0.013$, i.e. consistent with zero. In the pixels identified by the clustering algorithm we find $\left<\beta\right>=-2.970\pm 0.006$ and $\sigma_{\mathrm{intrinsic}}=0.03\pm 0.01$. These averages indicate that any bias on the mean of the posterior is small compared to the uncertainties. The low intrinsic scatters, which are consistent with zero, show that the uncertainties are representative of the spread of the estimates.

5 Results on real data

Maps of the estimated spectral index between the real Planck 30 and 44 GHz $P$ maps, the uncertainty of those estimates, and the Kullback-Liebler Divergence between the posteriors and the priors are shown in Figure 6.

Once again, we use the mean shift clustering algorithm to group the pixels into clusters. The largest cluster, which corresponds to those with good detections of the spectral index, contains 747 pixels (around one quarter of the sky). 2-dimensional histograms of $\mu_{\beta}$, $\sigma_{\beta}$ and $D_{\textrm{KL}}$ are shown in Figure 7, the locations of the pixels identified by the clustering algorithm are shown in red. The histograms of real data show similar structures as the simulated data (cf. Figure 5).

5.1 Average spectral index and intrinsic scatter

The weighted average spectral index across the whole sky is -3.38. This average is dominated by pixels along the Galactic plane with typically steep spectra and low uncertainty. Along the Galactic plane we expect complicated frequency spectra due to the superposition of radiation from multiple synchrotron sources along the line of sight experiencing different levels of Faraday rotation as well as turbulent magnetic fields introducing frequency dependent fluctuations (for example Lazarian & Pogosyan, 2016). The weighted average just including pixels that belong to the cluster is -3.09, and if then also ignoring pixels within $10^{\circ}$ of the Galactic plane the average is -2.96.

Again, by approximating the errors as Gaussian, we fit for both a mean value and the extra intrinsic scatter in the data that can not be accounted for by the errors alone. The means and intrinsic scatters, along with values found by others, for large sky areas are listed in Table 1. Direct comparison between the averages found in this work and by others is difficult due to the different sky areas considered. Estimates from data at lower frequencies typically find steeper spectral indices, perhaps indicating a shallowing of the spectrum with frequency.

Across the whole sky, we find $\left<\beta\right>=-2.99\pm 0.03$ and $\sigma_{\mathrm{intrinsic}}=1.12\pm 0.02$. The mean and intrinsic scatter of just the pixels identified by the clustering algorithm are $\left<\beta\right>=-3.12\pm 0.03$ and $\sigma_{\mathrm{intrinsic}}=0.64\pm 0.02$. If we also neglect pixels within $10^{\circ}$ of the Galactic plane we find $\left<\beta\right>=-2.92\pm 0.03$ and $\sigma_{\mathrm{intrinsic}}=0.48\pm 0.02$.

These intrinsic scatters are significantly larger than the value found from simulated data, which had no intrinsic scatter. When interpreting the intrinsic scatters, it is important to remember that we have approximated the extra spread of spectral indices as Gaussian. This is likely not the case and the intrinsic scatter values are therefore just indicative of the variation, $\sim 20$ per cent from the mean.

Vidal et al. (2015) used $T-T$ plots to estimate the polarized spectral index between the WMAP $K$-, $Ka$- and $Q$-band maps in eighteen regions. For each region in our spectral index map, we fitted for an average value and extra intrinsic scatter. The indices that we found along with those found by Vidal et al. (2015) are plotted in Figure 8 and listed in Table 3. We could not determine the average spectral index in four regions because three of them do not enclose the central coordinates of any $N_{\mathrm{side}}=16$ pixels (regions 6, 9 and 10) and one (region 8) contained no pixels identified by our clustering algorithm. We could not estimate the intrinsic scatter in three of the regions (regions 1, 11 and 18) as they only contained a single pixel identified by our clustering algorithm. The spectral indices found by us between Planck 30 and 44 GHz are generally consistent with those found by Vidal et al. (2015) between the WMAP maps.

5.2 Angular power spectrum

We use PolSPICE333http://www2.iap.fr/users/hivon/software/PolSpice/ (Chon et al., 2004) to estimate the angular power spectrum of the polarized spectral index map, masking $25^{\circ}$ either side of the Galactic plane. The spectrum is shown in Figure 9 and is consistent with both the default spectral index map used by PySM (from Miville-Deschênes et al., 2008) and with the power spectrum of the spectral index between 1.4 and 28.4 GHz found by Krachmalnicoff et al. (2018).

In Figure 9 the angular power spectrum of the mean spectral index map found in this work is shown in black. To estimate the impact of the noise on the spectrum we; draw 100 samples of the spectral index from the converged MCMC chains for each pixel; subtract the mean estimate of the spectral index in that pixel from the 100 samples; form 100 maps from these mean-subtracted samples and estimate their power spectra. These power spectra are shown in grey. At each multipole, the mean of these noise spectra is the noise bias and the standard deviation is an estimate of the uncertainty. The noise de-biased spectrum and $1-\sigma$ error bars are plotted incyan.

For comparison we have also plotted the power spectrum found by Krachmalnicoff et al. (2018) and the power spectrum of the spectral index map used by PySM by default. Both of these comparison spectra are consistent with the noise de-biased spectrum of the spectral indices found in this work.

5.3 Fermi bubbles/WMAP haze

There is a region north of the Galactic centre with high-significance detections of shallow spectral indices. A gnomonic projection of this region is shown in Figure 10.

Spatially the shallow spectrum emission correlates with the location of the northern Fermi bubble (the coordinates defining the Fermi bubbles are from Su et al. (2010) and outlined in red). The average spectral indices are different in the north and south bubbles. In the north bubble, only including pixels identified by the clustering algorithm, we find $\left<\beta\right>=-2.36\pm 0.09$ with an intrinsic scatter of $\sigma_{\mathrm{intrinsic}}=0.63\pm 0.07$. In the south bubble, again only including pixels identified by the clustering algorithm, we find $\left<\beta\right>=-3.00\pm 0.05$ and $\sigma_{\mathrm{intrinsic}}=0.35\pm 0.04$. Combining both regions together we find $\left<\beta\right>=-2.72\pm 0.06$ and $\sigma_{\mathrm{intrinsic}}=0.57\pm 0.04$.

For comparison, Planck Collaboration et al. (2012) included a diffuse 2-dimensional Gaussian centred on the Galactic centre as a template for this emission during template fitting. They found a spectra index for the haze component in total intensity of $\beta=-2.56\pm 0.05$, which is remarkably consistent with the average spectral index we find across both bubbles including all pixels (not just those identified by the clustering algorithm), $\left<\beta\right>=-2.57\pm 0.08$ and $\sigma_{\mathrm{intrinsic}}=1.06\pm 0.06$.

To determine whether the average spectral indices in the north and south bubbles are different we compute the Bayes factor. The Bayes factor is the ratio of the evidence (or marginalized likelihood) for the case where the true averages are different to evidence for the case where the true averages are the same. Bayes factors greater than unity are evidence in favour of the averages being different, factors less than unity are evidence for the averages being the same. Jeffreys (1983) suggests the following heuristic for interpreting Bayes factors; Bayes factors between 1–3 are barely worth mentioning, Bayes factors between 3–10 are substantial, factors between 10–30 are strong evidence, factors between 30–100 indicate very strong evidence, and factors greater than 100 are decisive evidence. The Bayes factor in this instance is $3\times 10^{7}$.

The shallow spectral indices that we have found in the northern bubble could be caused by several effects and not all of them are astrophysical. Inaccurate $\sigma_{P}$ values, incorrect CMB subtraction, $I$ to $P$ leakage at 44 GHz, depolarization at 30 GHz, zero-level errors in the prior, anomalously highly-polarized free-free emission or shallow spectrum synchrotron emission could be responsible. We discuss each of these in turn below.

5.3.1 Inaccurate $\sigma_{P}$ values

The $\sigma_{P}$ values for the Planck maps could underestimate the true noise levels in the maps. To determine the impact of using higher noise levels on the estimated spectral indices, we scaled the $\sigma_{P}$ values for both Planck maps by a factor of two and a factor of ten, estimated the spectral index in each pixel and found the average in the Fermi bubbles. The values are listed in Table 4.

Without scaling the noise levels, the evidence of the difference is decisive with a Bayes factor of $3\times 10^{7}$, scaling the noise levels by a factor of two only marginally decreases the evidence to $3\times 10^{6}$, which is still decisive evidence for the averages being different. The noise has to be scaled by very large factors to make the significance of the difference drop substantially. When the noise is scaled by a factor of 10 the Bayes factor falls to around 4, just considered substantial by Jeffreys (1983). We do not expect the noise levels on the maps to be wrong by such a large factor, but even if they are we still find substantial evidence for the average spectral index in each bubble being different.

Unsurprisingly, we find that as the noise scaling increases the combined average spectral index across both bubbles moves closer to the prior ($-3.11$).

5.3.2 CMB subtraction

The SMICA estimate of the CMB will be contaminated by residual foregrounds, particularly close to the Galactic plane. We re-ran the fitting on Planck maps without the CMB subtracted and instead, to account for the fluctuations in the CMB, marginalized over the CMB amplitude in each pixel with a Rayleigh prior

[TABLE]

where $A_{\mathrm{CMB}}^{\prime}$ is the amplitude of the CMB fluctuation (in thermodynamic units) and $\sigma_{\mathrm{CMB}}^{\prime}=0.27,\mu\mathrm{K}_{\mathrm{thermo.}}$ is the standard deviation of the CMB polarization fluctuations (again in thermodynamic units). 444The quadrature sum of two Gaussian random variables with zero mean and common variance is a new random variable with a Rayleigh distribution.

Of all the tests we made to the robustness of our estimated spectral index map, this had the largest impact on the estimated spectral indices. A scatter plot of the estimated spectral indices when subtracting the SMICA CMB estimate to the indices when instead marginalizing over the CMB amplitude is shown in Figure 11. The points lie along the line $y=0.99x-0.12$, i.e. a line with slope very close to unity but a non-zero offset. By marginalizing over the CMB amplitude we shift the spectral index posterior distributions to steeper values and also widen them by around 13 per cent.

None the less, we still find significant differences between the average spectral index in the region of the north and south fermi bubbles. The mean spectral indices in the bubbles and the Bayes factor that quantifies the evidence for the means being different, calculated without the SMICA CMB estimate subtracted and instead marginalizing over the CMB amplitude, are listed in the first column of Table 4. Subtracting the SMICA CMB estimate shifts the averages to steeper values and reduces the Bayes factor to $4\times 10^{5}$, which is still decisive evidence that the averages are different.

5.3.3 $I$ to $P$ leakage in the 44 GHz map

$I$ to $P$ leakage in the Planck 44 GHz map would result in shallower spectral indices, particularly in the north where the emission is brighter. There are two reasons we believe this to be an unlikely cause; firstly the bandpass mismatch $a$-factors for the 44 GHz horns are less than those for the 30 GHz horns (Planck Collaboration et al., 2018a), secondly typical foreground spectra are falling over this frequency range and so would also require the total intensity foreground emission in this region to be anomalously bright at 44 GHz compared to 30 GHz.

5.3.4 Depolarization in the 30 GHz map

Depolarization, caused by Faraday rotation in the Planck 30 GHz survey, could lead to the measurement of spuriously shallow synchrotron spectral indices. Faraday rotation effects are expected to be more significant close to the Galactic plane, however these same effects are expected to be small at frequencies $\sim 30\,\mathrm{GHz}$. In the worst affected regions directly in the Galactic plane, we may see Faraday depths $\phi\approx 500\,\mathrm{rad/m^{2}}$ (Wolleben et al., 2019; Oppermann et al., 2012). Assuming the polarization angle varies linearly with the square of the wavelength, $\lambda^{2}$, we can estimate the worst Faraday rotation effects to induce position angle rotations $\Delta\psi\sim 3^{\circ}$. It is worth noting here that, in reality the position angle dependence for synchrotron emission will not vary linearly with $\lambda^{2}$. However this estimate is at least indicative of the low level of expected depolarization even in the most affected parts of the sky. Further, observations of the Faraday depth across the sky show that the level of Faraday rotation is even lower over the full extent of the Fermi bubbles. For these reasons it seems unlikely that the consistent estimation of shallow spectral indices in the northern bubble is the result of depolarization at $30\,\mathrm{GHz}$.

5.3.5 Zero level of prior

A systematic offset in the expected polarized intensity could result in shallower spectral indices when the signal-to-noise in the 44 GHz map is less than that in the 30 GHz map. To test the impact of this we subtracted an 8.9 K offset from the Haslam total intensity map that formed the prior. This is the monopole found in the Haslam map by Wehus et al. (2017) using $T-T$ plots.

The average spectral indices are hardly affected by this change and actually the evidence for the average in the bubbles being different increases by around a factor of two.

5.3.6 Astrophysical source

Finally, the measured spectral indices may in fact be caused by some astrophysical source. The hard-spectrum emission in the northern bubble is unlikely to be polarized free-free emission. Firstly, free-free emission is typically observed to have a spectral index of $\beta\approx-2.1$, whilst we find an average spectral index in the northern bubble of $\beta=-2.36\pm 0.09$. Secondly, free-free emission is expected to be intrinsically un-polarized due to the random scattering directions giving rise to the emission. Free-free emission can be polarized at the edges of optically dense H ii clouds, where the polarization fraction could reach $10\%$ (Rybicki & Lightman, 1985; Keating et al., 1998). However, generally the polarization fraction will be close to zero.

In previous analyses with WMAP and Planck total intensity data, it was argued that the Galactic Haze emission likely arose from a hard synchrotron component. In our analysis we find the polarized spectral indices in the southern bubble to be similar to the soft synchrotron emission that dominates over most of the sky. The spectral indices in the northern bubble are harder. It is important to note that these regions have many overlapping features and along the line-of-sight these features can sum to give unusual spectra, particularly in polarization (for example Lazarian & Pogosyan, 2016). However, we have significant detections of shallower spectral indices in the region defined by the northern Fermi bubble.

6 Conclusions

We have used Bayesian methods to estimate the spectral index between the Planck 30 and 44 GHz polarized intensity maps. We modelled the polarized intensities as Rician random variables and used a prior formed from the multiplication of the objective reference prior for the problem with a maximum entropy prior that incorporates total intensity information. We tested our method on simulated data and found it gave unbiased estimates with uncertainties representative of the scatter on the data.

In each pixel we calculated the mean, standard deviation and Kullback-Liebler divergence from the prior of the spectral index posterior distribution. We used a clustering algorithm on these parameters to identify pixels with good detections of the spectral index, and found good detections in 747 pixels (which corresponds to approximately one quarter of the sky).

We fitted for the mean and intrinsic scatter of the spectral index, assuming that the true spectral index across the sky is a Gaussian random variable and that the posterior distributions are Gaussian. The global average spectral index that we found both across the whole sky and in pixels identified by the clustering algorithm are consistent with values found by others (these are summarised in Table 1). Across all pixels we found an average $\left<\beta\right>=-2.99\pm 0.03$ and in pixels identified by the clustering algorithm we found an average of $\left<\beta\right>=-3.12\pm 0.04$. The angular power spectrum of our spectral index map is consistent with the spectrum found by Krachmalnicoff et al. (2018) (Figure 9).

We found a statistically significant difference between the average spectral index in the North and South Fermi bubbles. Only including pixels identified by the clustering algorithm, the average spectral index in the Northern bubble is $\left<\beta\right>=-2.36\pm 0.09$, the average spectral index in the Southern bubble is $\left<\beta\right>=-3.00\pm 0.05$, and the average spectral index across both bubbles is $\left<\beta\right>=-2.72\pm 0.06$. The average across both bubbles including all pixels, not just those identified by the clustering algorithm, is $\left<\beta\right>=-2.57\pm 0.08$ and this is consistent with the spectral index found by Planck Collaboration et al. (2012) who treated both bubbles as a single component and found a spectral index of $-2.56\pm 0.05$.

The largest potential systematic error in this analysis is the CMB subtraction step. Wen we marginalize over the CMB amplitude (instead of subtracting an estimate) the indices are typically shifted to more negative values, $\Delta\beta=-0.12$, and the uncertainties are increased by around 13 per cent.

In future work, repeating this analysis with higher signal-to-noise maps will increase the fraction of the sky with detections and the weakly informative part of the prior could be improved with a spatially varying temperature spectral index and polarization fraction.

Acknowledgements

We would like to thank Clive Dickinson and Tim Pearson for comments on a draft of this paper. We would also like to thank Clive Dickinson further for useful discussions about diffuse Galactic emission at radio frequencies.