Testing for spectral index variations in polarised CMB foregrounds

We present a Bayesian parametric component separation method for polarised microwave sky maps. We solve jointly for the primary cosmic microwave background (CMB) signal and the main Galactic polarised foreground components. For the latter, we consider electron-synchrotron radiation and thermal dust emission, modelled in frequency as a power law and a modified blackbody respectively. We account for inter-pixel correlations in the noise covariance matrices of the input maps and introduce a spatial correlation length in the prior matrices for the spectral indices beta. We apply our method to low-resolution polarised Planck 2018 Low and High Frequency Instrument (LFI/HFI) data, including the SRoll2 re-processing of HFI data. We find evidence for spatial variation of the synchrotron spectral index, and no evidence for depolarisation of dust. Using the HFI SRoll2 maps, and applying wide priors on the spectral indices, we find a mean polarised synchrotron spectral index over the unmasked sky of beta-sync = -2.833 +- 0.620. For polarised dust emission, we obtain beta-dust = 1.429 +- 0.236. Our method returns correlated uncertainties for all components of the sky model. Using our recovered CMB maps and associated uncertainties, we constrain the optical depth to reionization, tau, using a cross-spectrum-based likelihood-approximation scheme (momento) to be tau = 0.0598 +- 0.0059. We confirm our findings using a pixel-based likelihood (pixlike). In both cases, we obtain a result that is consistent with, albeit a fraction of a sigma higher than, that found by subtracting spatially uniform foreground templates. While the latter method is sufficient for current polarisation data from Planck, next-generation space-borne CMB experiments will need more powerful schemes such as the one presented here.


INTRODUCTION
Measurements of the polarised cosmic microwave background (CMB) anisotropies provide a unique opportunity to constrain early Universe physics. However, Galactic polarised emission must be removed to high accuracy to extract reliable information from the primordial CMB. Disentangling the polarised CMB signal at large angular scales from its contaminants, mainly synchrotron emission and thermal dust emission, is crucial for measurement of the reionization optical depth and the detection of a signal from primordial gravitational waves. Scalar perturbations generated during inflation give a divergencelike -mode polarisation pattern in the CMB, whereas tensor perturbations would produce a distinctive curl-like -mode polarisation signature, together with an -mode pattern of approximately equal amplitude (Kamionkowski et al. 1997;Seljak & Zaldarriaga 1997). Measuring primordial -modes would fix the energy scale of inflation (Starobinskii 1983) and provide constraints on specific inflationary models (e.g. Kallosh & Linde 2019). Gravitational lensing of CMB -modes by intervening matter generates a -mode ★ E-mail: rmvd2@cam.ac.uk anisotropy contributing an additional contaminant to the primordial B-mode signal 1 (see Lewis & Challinor 2006, for a review).
A key science goal of next-generation CMB experiments from space (LiteBIRD; Sugai et al. 2020) and from the ground, such as CMB-S4 (The CMB-S4 Collaboration et al. 2022) and the Simons Observatory (SO, Ade et al. 2019) is to achieve a statistically significant detection of primordial gravitational waves if the tensor-to-scalar ratio is of order 10 −3 (Ade et al. 2019). Polarised foregrounds pose a major challenge for these ambitious projects. It is therefore important to investigate Galactic foregrounds in as much detail as possible using existing polarisation data. In this paper we use Planck data at large angular scales, which poses additional challenges because of the low signal-to-noise of the Planck polarisation maps and the presence of residual instrumental systematics.
We adopt a Bayesian approach to polarized foreground removal. The aim is to recover the polarised CMB anisotropies, together with a parametric model of polarised Galactic foregrounds, using a statistical framework that allows the foreground errors to be propagated through to cosmological parameters. We assume that lower frequencies ( < ∼ 70 GHz) are dominated by synchrotron emission while higher frequencies ( > ∼ 100 GHz) are dominated by thermal dust emission. We approximate the frequency dependence of the polarised foregrounds as a power law for synchrotron and a modified blackbody for thermal dust. Aspects of the theoretical basis for the approach used in this paper were originally developed in e.g. Eriksen et al. (2006); Dunkley et al. (2009b); Stompor et al. (2009) and Gratton (2008). We extend this framework by allowing for pixel-pixel correlations in the noise covariance matrices (NCM) and imposing priors on the variations of the spectral indices across the sky, effectively regularising noise in the maps.
For Low Frequency Instrument (LFI) Planck data (Planck Collaboration 2020a), we assume that polarised emission arises primarily from synchrotron radiation emitted by spiralling cosmic ray electrons in the Galactic magnetic field. Synchrotron emission can be highly polarised, with polarisation fractions of up to 75% (Ginzburg & Syrovatskii 1965). The specific intensity, , of synchrotron emission has power law frequency dependence, ∝ , with ∼ −3. The spectral index is related to the energy distribution of the cosmic rays of cosmic ray electrons, / ∝ − , according to = −( + 3)/2 (Rybicki & Lightman 1985). There have been many studies of angular variations of he spectral index . Radio frequency observations in the frequency range 408 MHz to 10 GHz (Davies et al. 1996;Platania et al. 2003;Bennett et al. 2003) find values of in the range -2.8 to 3.2 with shallower values towards the Galactic plane. This is consistent with models in which the steepening of the spectral index is caused by energy dependent diffusion and radiative losses as electrons propagate into the Galactic halo (see e.g. Dickinson 2018). There have also been analyses reporting steeper synchrotron spectral indices in the radio loops, perhaps reflecting an upper limit to the electron energy spectrum generated by shock acceleration within the loops (Lawson et al. 1987;Reich & Reich 1988;Borka 2007). Galactic foreground emission in total intensity over the frequency range 20 − 100 GHz is discussed in detail by Planck Collaboration (2016d) and is difficult to interpret because of the complexities associated with disentangling synchrotron from free-free and anomalous microwave emission. Polarised Galactic emission over this frequency range, which is the subject of this paper, is dominated by synchrotron emission. However, polarization maps from Planck and the Wilkinson Microwave Anisotropy Probe (WMAP; Hinshaw et al. 2013) are noisy and affected by instrumental systematics. Previous work has largely been limited to measuring the synchrotron spectral index leading to values ≈ −3 (consistent with the synchrotron intensity spectral index) with some tentative evidence of spectral index variations (Page et al. 2007;Fuskeland et al. 2014;Vidal et al. 2015;Dickinson 2018;Svalheim et al. 2020;Fuskeland et al. 2021;Martire et al. 2022).
At frequencies > ∼ 217 GHz covered by the Planck High Frequency Instrument (HFI) (Planck Collaboration 2020b; Delouis et al. 2019) thermal Galactic dust emission dominates the foreground emission in intensity over most of the sky, with the cosmic infrared background (CIB) dominating at high Galactic latitudes. Galactic dust emission has been discussed by the Planck collaboration (Planck Collaboration 2014c, 2015). Using map-based statistics, these papers find that the spectral energy distribution (SED) of Galactic dust emission in intensity is well approximated by a modified black body distribution with a spectral index = 1.59 ± 0.12, and dust temperature = 20.3 ± 1.3K. A power-spectrum based analysis, which allows accurate separation of Galactic dust emission from the CIB and primordial CMB leads to consistent results of = 1.49 ± 0.05 and = 22.7 ± 2.8K ). There is no evi-dence at the Planck frequencies for angular variations of the diffuse Galactic dust emission in intensity.
Since the highest polarised Planck channel is at 353 GHz, there is a more restricted frequency range available with which to study Galactic dust emission in polarization compared to intensity. Nevertheless, the polarised dust SED is consistent with a modified black body distribution with the same values of and as in intensity (Planck Collaboration 2015). Variations in the dust temperature and orientation of dust grains with respect to the Galactic magnetic field can lead to line-of-sight differences in the polarised dust SED (Tassis & Pavlidou 2015). An analysis of the Planck BB power spectrum suggested a decorrelation of polarised dust spectrum between 217 and 353 GHz (Planck Collaboration 2017). However a more detailed analysis of the Planck maps, using more robust statistical techniques, showed no evidence for any decorrelation (Planck Collaboration 2020f). More recently Pelgrims et al. (2021) have used a combination of Planck polarization data and measurements of neutral hydrogen (HI) velocity components to find that lines of sight with multiple HI components display greater depolarisation compared to other lines of sight.
Many component separation methods have been developed to disentangle the primary CMB signal from foregrounds (see e.g. Planck Collaboration 2020c, for examples). The methods can be roughly split into three categories: template fitting, parametric fitting methods and multi-frequency internal linear combination schemes. The analysis of Planck polarization maps has been based on subtracting a low frequency synchrotron template and a high frequency Galactic dust template (Planck Collaboration 2016f, 2020dPagano et al. 2020;de Belsunce et al. 2021). Bayesian parametric fitting methods, such as the Markov Chain Monte Carlo algorithm Commander (Eriksen et al. 2004(Eriksen et al. , 2008, have been applied successfully to disentangle multipole foreground components in intensity (see e.g. Planck Collaboration 2020c). However these methods have been much less successful when applied to large-scale polarization because of the higher noise levels (leading to convergence problems) and the contribution of instrumental systematics. For an early application of pixel-by-pixel parametric foreground fitting to WMAP polarization data see Page et al. (2007). In this paper, we develop and apply stable Bayesian parametric methods to the removal of Galactic foregrounds from the Planck polarization maps. This paper is organised as follows: in Sec. 2 we introduce our Bayesian parametric fitting method and in Sec. 2.2 we present our likelihoods. In Sec. 3 we discuss the Planck data that we use and the pixel-pixel noise covariance matrices required for the Bayesian component separation method. In Sec. 4 we test the parametric foreground fitting approach on Planck data and use two likelihoods to quantify its effect on the optical depth to reionization parameter. Sec. 5 presents our conclusions.

METHODS
Following the formalism discussed in e.g. Eriksen et al. (2008) where P ( |d, M) is the posterior, P (d| , M) the likelihood, P ( |M) the prior and P (d|M) the evidence.
The Planck satellite provides , , and maps, each containing pix pixels, in frequency bands. We assume Gaussian noise that is uncorrelated between frequency channels but may have inter-pixel correlations. Fitting parametric foreground models is computationally intensive at full Planck resolution. Since in this paper we are interested in polarized foregrounds at large angular scales, we degrade the Planck LFI and HFI maps to a resolution of side = 16 containing pix = 12 · 2 side = 3072 pixels following the procedure described in de Belsunce et al. (2021, hereafter B21).

Bayesian parametric foreground estimation
We assume that the observed sky can be modelled as where d are the observed (masked) sky maps in antenna temperature, Am is the sky signal which consists of a frequency-scaling matrix A and model amplitudes m, and n is the instrumental noise. The parametric model for the sky signal is composed of three components: the primordial CMB ( ), electron-synchrotron emission ( ), and thermal dust radiation ( ). Foregrounds are taken to be independent of the CMB and we model the contaminants, namely electron synchrotron and thermal dust, as a power law and a modified blackbody respectively, see e.g. Planck Collaboration (2016e, 2020c). Given Eq. (2), we disentangle the CMB and the foregrounds by jointly estimating the posterior distribution P (A, m|d) of the scaling matrix and the model amplitudes given the observed data.
The position-dependent scaling matrix A describes the frequency dependence of the CMB and foreground signals. The first columns of A consist of repeated block-diagonal identity matrices 1, weighting the primordial CMB, and the following columns are populated with factors weighting the foreground channels. The vector m consists of var stacked amplitudes of the model components: m with = { , , }. The sky signal can thus be written in thermodynamic temperature as the sum over the different components where 0 ( ∈ { , }) are the pivot frequencies assigned to the channels assumed to best trace synchrotron and dust, respectively. The function ( ) relates thermodynamic to brightness temperature and ( ) is a frequency-dependent function accounting for a modified blackbody. The arguments are = ℎ / and = ℎ / , with the Boltzmann, ℎ the Planck constant and ( ) the CMB (dust) temperature 2 .
We introduce a physically motivated prior term for the spectral indices to vary continuously across the sky. Following the procedure outlined in Crowe (2013), we therefore introduce a prior for the spectral indices incorporating a spatial correlation angle corr , which adds non-zero off-diagonal entries to the prior: (4) Here the angular separation 3 between the pixels across a sky map 2 For the exact derivation and numerical values of the conversion factors, see appendix A. 3 Taking the limit of a very large correlation angle reduces the scheme to a spatially uniform template cleaning procedure. is given by , and ( ) 2 is the variance per spectral index. The negative log-prior term thus reads 1 2 where is a stacked vector of the spectral indices ( , ).
To determine the best solution, we minimise the negative loglikelihood (S ≡ − ln P): where N is the noise covariance matrix (NCM) which we will introduce in Sec. 3.2. The prior term acts as a penalty function on spectral index variations and effectively regularises scatter in the recovered maps. We solve this minimisation problem partly analytically and partly numerically. Uncertainties in our solution are estimated as the inverse of the curvature matrix of Eq. (6).

Implementation details
We sample the spectral indices jointly for sp = 2 Planck polarisation maps, namely and , with pix = 3072 pixels in = 6 frequency channels: For LFI, we use 30 and 70 GHz and for HFI 100, 143, 217, and 353 GHz. We set the synchrotron pivot frequency to 30 GHz and the dust one to 353 GHz. We neglect noise correlations between frequency channels, using block-diagonal noise matrices in frequency space. The consideration of inter pixel noise covariance and correlated spectral index priors renders the minimisation problem computationally labelled collectively as x, and our scheme requires inverting square matrices of corresponding size. The minimisation is performed by a multi-dimensional pseudo-Newton-Raphson (NR) method (Press et al. 2007) following a twostep procedure. First, for given spectral indices, we analytically compute the amplitudes that minimize Eq. (6). Second, we compute first and second derivatives of Eq. (6) for all parameters and evaluate these at the assumed 's and the calculated amplitudes. Using these derivatives, we take a step in parameter space Δx = − H −1 · ∇S, where is the Hessian of the likelihood at the current trial point, and ∇S the vector of first derivatives of Eq. (6). We iterate until convergence is achieved. The parameter ∈ [0, 1] is adaptively varied during the process to speed up convergence.
For an initial starting point, we chose to use the 30 GHz maps for the synchrotron amplitudes, the 100 GHz maps for the CMB, and 353 GHz maps for the thermal dust amplitudes. For the spectral index maps, we chose to start at the peak of the assumed prior, to be discussed below.
In our analyses we choose the mask presented in Fig. 1 which retains a fraction of = 0.52 of the sky (effects of mask choice are discussed in appendix C). For some choices of broad spectral priors, the scheme did not immediately converge. In these situations, we first solved with a sufficiently tight prior to achieve convergence. Then we gradually increased the prior towards the desired value, and 4 The presented approach is still much less expensive than Gibbs sampling approaches, see e.g. (Eriksen et al. 2008;Armitage-Caplan et al. 2011). For low-resolution polarisation Planck data with pixel-pixel correlations, the foreground cleaning takes a few CPU-hours at side = 16.
using the previous solution as the starting point for the next one. Newton-Raphson-type methods are not guaranteed to converge and we saw examples of this phenomenon for a small number of pixels (∼ 10) at the boundary of our original mask. We prevented this issue from occurring in our analyses by masking these pixels together with their neighbours (for a total of ∼ 50 additionally masked pixels).

Noise and prior correlations
To explore the sensitivity of our scheme to inter pixel correlations, we test three different cases. First, we take the noise covariance matrix (NCM) N and prior matrix C to be diagonal; we denote this 'diagNCM + diagPrior'. This simplification renders the problem computationally trivial as the minimisation can be carried out pixel by pixel. Second, we allow for pixel-pixel correlations in the NCM while keeping the prior matrix diagonal; we call this 'corrNCM + diagPrior'. Third, we allow for inter pixel correlations in both the NCM and the prior matrix; we denote this 'corrNCM + corrPrior'.
The resulting spectral index maps naturally depend on the spectral index priors because the data is not very constraining. Therefore, it is critical to choose physically motivated priors 5 . We introduce wide priors following the form of Eqs. (4) and (5). These are centred around mean values of previous spectral index measurements from Planck Collaboration (2020c) and radio frequency observations. Hence, we use sync = −3.0 with a prior width of ( sync ) = 2.0, and dust = 1.5 with a prior width of ( dust ) = 1.0, for synchrotron and thermal dust respectively.
Significant variations of the spectral indices on small scales are unlikely because of the large interaction length of cosmic rays measured with gamma ray observations. Together with evidence for spatially varying spectral indices measured at 400 MHz -10 GHz from radio frequency observations we are confident that spectral index variations occur at scales no smaller than 5 − 10 • . Therefore, we choose an angular correlation length for both foreground components of corr = 5 • .

Likelihoods
In the following, we briefly introduce and compare the two likelihoods used to assess the performance of the foreground cleaning scheme for cosmological parameter estimation. We do this by constraining the optical depth to reionization parameter from Planck polarisation maps and compare our results to the templatecleaned analysis presented in B21. In Sec. 2.2.1 we present a crossspectrum-based semi-analytic likelihood approximation scheme, called momento. In Sec. 2.2.2 we present a polarisation-only pixelbased likelihood.

Low multipole likelihood-approximation scheme
We present a cross-spectrum-based likelihood-approximation scheme (momento) that is based on the 'General Likelihood Approximate Solution Scheme' ( ; Gratton 2017). allows to approximate the sampling distribution by computing moments and cumulants of the data. This is relevant in situations where the exact 5 One can easily analytically quantify the effect of a tight Gaussian prior combined with a large uncertainty from the data. Since the data is underconstraining the problem, the tight prior will lead to spectral index maps centred around the input prior value with only very low scatter. Adopting priors that are too tight can lead to results which are highly misleading. likelihood is unknown or computationally too expensive to compute. momento was extensively tested on simulations and data in the context of low-ℓ CMB analysis in B21. In practice, the moments can be either computed analytically or through forward simulations, to compute the gradient (denoted by , ) of the log-likelihood: This procedure is generally applicable to a multitude of statistics , ∈ [1, . . . , ], and denotes a vector of powers of varying degree of these statistics. The ideally should be chosen to be as close to summary statistics as possible, e.g. elements of the power spectrum in the present case. It is also important to be able to obtain the required moments of the as a function of the model either analytically or with forward simulations.
For momento we choose quadratic cross-spectra (QCS), compressing the CMB data down to a set of multipoles, here ℓ = 2 − 29. A fiducial power spectrum fid ℓ is assumed for the purpose of computing the QCS. Subsequently, for each evaluation of the likelihood the following actions are taken by momento: (i) take input set of theory ℓ 's; (ii) compute the residual between fiducial and theory spectra Δ ℓ = ℓ − fid ℓ ; (iii) perform a path integral for the change of the log-likelihood in parameter space from fid ℓ to ℓ . We choose the path fid ℓ + Δ ℓ in QCS space, where ∈ [0, 1]. This requires computing the gradients of the log-likelihood, given in Eq. (8), along this path.
To apply momento, pixel-pixel covariance matrices are required to compute moments of the QCS of the input maps. In B21, we used NCMs of the frequency maps to compute . In the present paper, momento uses the CMB 'block' of the inverse of the curvature covariance matrix which incorporates the uncertainties from the foreground cleaning scheme.
In Sec. 4.3 we will show results for using momento on templatecleaned and CMB maps recovered from the presented component separation scheme. In the present paper, we compute only using the QCS spectra.

Polarisation-only pixel-based likelihood
Under the assumptions that signal and noise are Gaussian and the noise covariance matrix is known accurately, it is possible to develop a pixel-based likelihood for low resolution CMB maps, see e.g. Page et al. (2007). This likelihood is formally exact, even for incomplete sky coverage, and is conceptually straightforward. A drawback of this likelihood is that the NCM has to be accurately known. Further, when using a pixel-based likelihood, it is difficult to assess whether or not a given result on is being driven by outliers at one or two multipoles.
The pixel-based polarisation-only likelihood is given by the probability of the data given the theoretical noise and signal model: where d is the data vector and d its smoothed mean, both containing the polarisation and maps. This procedure can be trivially extended to account for temperature as well. However, in the present context we only focus on polarisation data. C is the covariance matrix which consists of a signal S and noise N component C = S + N. The signal matrix can be constructed following the procedure outlined in e.g. Tegmark & de Oliveira-Costa (2001). with sky = 0.85 is used as a 'processing' mask to perform the degrading of the high resolution maps to low resolution (as discussed in the text). The mask with sky = 0.52 is binary and used for foreground cleaning and power spectrum computations.

DATA
We use various map products based on the time-

Map compression
For the degradation process, we first apply an apodised mask with sky = 0.85 (as plotted in Fig. 1) to the high resolution Planck and maps to suppress the Galactic plane region and avoid smearing high amplitude foreground emission to high Galactic latitudes. The maps are then smoothed using the harmonic-space smoothing operator: with ℓ 1 = lr side and ℓ 2 = 3 lr side . To be compatible with certain covariance matrix products, discussed below, we further smooth the maps with a HEALPix pixel window function. Finally, we repixelise the maps at lr side = 16 (with pix = 3072 pixels).

Pixel-pixel noise covariance matrices
The fidelity of our component separation procedure depends on how well the noise covariance matrices (NCMs) model the noise and systematics in the data. We closely follow the procedure used in B21, extending the treatment from 100 and 143 GHz to all the LFI and HFI frequencies used in this paper. We summarise only the key steps here.
For the NCMs, we first fit a parametric model to the covariance matrix computed from end-to-end simulations for LFI and HFI data 7 . The empirical noise and systematics only covariance matrix 6 Recently, two new sets of Planck maps have been introduced for LFI and HFI jointly (NPIPE, Planck Collaboration 2020g) and for LFI only (BeyondPlanck, BeyondPlanck Collaboration 2020). With a computation of corresponding noise covariance matrices, an analysis such as that presented here could be repeated with these maps. 7 We omit 100 end-to-end simulations, to test the robustness of our NCMs.  Table 1. Summary of mean and standard deviations of the spectral indices for SRoll1 and SRoll2 maps using a mask with sky 0.52. We show the three different cases that differ from each other by their degree of inter-pixel correlations in the noise covariance and prior matrices.
is constructed from simulations viaN = −1 (n −n) (n −n) where n are the simulated sky maps andn a smoothed template for and Stokes parameters. The smoothed template corresponds to a smoothed mean of the maps (i.e. a maximum likelihood solution of the map-making equation, see e.g. Tegmark (1997)).
We approach the problem of fitting a model toN as a maximum likelihood inference problem. We assume a Gaussian probability distribution for each noise realisation where M is the model for the noise covariance matrix. We assume that M consists of two terms where the parameter scales the low resolution (lr) map-making covariance matrix N 0 and Y Y models additional large-scale effects. N 0 is constructed from FFP8 end-to-end simulations (Planck Collaboration 2016b,c) to capture the scanning strategy, detector white noise and '1/f'-type noise. The term Y Y models largescale modes with parameters as described in more detail in B21.
To obtain a NCM, we minimise the negative log-likelihood of Eq.
(11) with respect to and M.

RESULTS
In this section, we present results of our component separation method on low-resolution polarised Planck 2018 LFI and HFI data. In Sec. 4.1 we discuss the recovered spectral index maps for electron synchrotron and thermal dust. We compare in Sec. 4.3 the effect of the presented foreground scheme to a spatially-uniform template cleaning method (TC) on cosmological parameter estimation. Therefore, we measure the optical depth to reionization, , with a crossspectrum based likelihood-approximation scheme (momento) and a pixel-based likelihood (pixLike), using our recovered CMB products.

Synchrotron and dust spectral index maps
In Table 1 we present the mean and standard deviation ( ) of the recovered spectral index maps. These are obtained for synchrotron and dust from our component separation method when using either SRoll1 or SRoll2 maps for HFI. Since the same Planck 2018 LFI data is used in all cases, we expect differences in dust and only minor changes in sync . We present our results for three cases with varying degrees of inter pixel correlations: (i) diagonal NCM and diagonal prior matrix (diagNCM + diagPrior); (ii) inter pixel correlations in the NCM and a diagonal prior matrix (corrNCM + diagPrior); (iii) inter pixel correlations in both the NCM and prior matrix (corrNCM + corrPrior). In our analyses, we impose wide priors of sync = −3.0 with ( sync ) = 2.0 and dust = 1.5 with ( dust ) = 1.0. All the recovered spectral indices given in Table 1 are in good agreement and return average values of sync ≈ −2.93 and dust ≈ 1.42 with large average spreads of 0.69 and 0.32 respectively. We adopt the results of the corrNCM + corrPrior case applied to full-mission HFI SRoll2 maps as our fiducial analysis. For this case, we find the following means and standard deviations for the spectral indices: We show the histograms of the recovered synchrotron and dust spectral index maps in Fig. 2 for the three cases when using HFI SRoll2 data. (The results for the same analyses using SRoll1 for HFI are shown in appendix B.) As expected, assuming independent pixels (diagNCM + diagPrior) yields the largest dispersion in values. Introducing inter-pixel correlations, we find statistically insignificant shifts in the mean spectral indices, though the dispersions are reduced. The effects of additionally including inter pixel-correlations in the prior matrix are minor compared to the case where these are only present in the NCM.
Our mean recovered spectral indices given in Eqs. (13a)-(13b) are in agreement with, yet lower by ∼ 0.8 than, dust ≈ 1.60 obtained for LFI (Svalheim et al. 2020) and dust ≈ 1.55 for HFI data (Planck Collaboration 2016d, 2020c. Similarly, the SRoll2 collaboration finds a narrow distribution for dust (Delouis et al. 2021). However, the spread of our spectral index maps is wider by a factor of ∼ 10 − 15 for synchrotron and ∼ 4 for dust. These differences stem mostly from our large prior widths that allow the 's to explore a wider parameter space. Our priors are significantly broader than those adopted in previous studies, see e.g. Dunkley et al. (2009b); Armitage-Caplan et al. (2011);Fuskeland et al. (2014);Svalheim et al. (2020), which are in alignment with the Planck priors (Planck Collaboration 2016e, 2020c) 8 . These analyses use tight Gaussian priors of sync = −3.10±0.10 and dust = 1.56±0.10. We test our component separation method for these priors and recover similar priordriven results for SRoll2 full-mission data¯s ync = −3.000 ± 0.008 and¯d ust = 1.494 ± 0.028 using the pixel-by-pixel method. The recovered spectral index maps are almost identical when including inter pixel correlations. Such priors are so tight that they do no allow any significant exploration of the parameter space, which results in almost automatically rejecting the hypothesis of spectral index variations across the sky. Therefore we emphasise that previous results are strongly prior dominated. We next investigate angular structure in our recovered spectral index maps.

Spectral index variations
In Fig. 3 we show the recovered synchrotron and dust spectral index maps. One immediately notices the apparent feature in the upper right quadrant with particularly noticeable outliers in the spectral indices in the independent pixel case 9 . These pixels are located in the tails of their respective histograms. To investigate further, we look at the uncertainties our method yields for all of the fitted parameters with their covariance being given by the inverse of the curvature matrix evaluated at the solution. In Fig. 4 we plot the square root of the diagonals of the spectral index covariance matrices. We note the increased uncertainties in the recovered spectral indices in the upper right quadrant of the sky map, decreasing their apparent significance. In fact, in the diagNCM + diagPrior case one can see that the recovered uncertainties are approaching their maximal values (that of the prior) over much of the sky for the synchrotron spectral index; this indicates that the data is not constraining sync at all in the upper right quadrant.
More physically, the synchrotron spectral index maps in Fig. 3 display ring-like structures with steeper sync . These features correlate with the positions of radio loops, shown in Fig. 5, which are nearby structures within a few hundred parsecs (for a review, see e.g. Dickinson (2018)). Electron synchrotron radiation in these loops comes from young, fast cosmic rays which are remnants of supernovae. The uncertainties from the curvature matrix, shown in Fig. 4, are smaller within the loops, effectively increasing the statistical significance of spectral index variations within the loops.
To quantify these features further, we quote the mean and standard deviations of the spectral index maps within the radio loops in Table  2. For synchrotron, the means of loops I, III, IIIs, VIIb, and IV deviate from the mean across the unmasked sky by just under one standard deviation for a single pixel. With each loop containing ∼ 30 − 50 pixels and if each pixel's spectral index is assumed to be independent, one would expect the mean for each loop to agree with that across the whole sky up to a standard deviation / √ . Thus, we assess the evidence for spatial variations to be at the twoto four-sigma level. As expected, the pixel-by-pixel procedure shows stronger spatial variations of sync than the correlated methods. This comes from the inter pixel correlations effectively smoothing out features and penalising neighbouring pixels with (strongly) dissimilar -values. In contrast to synchrotron, visual inspection of the dust 9 This sky region is where differences in the polarization maps between different processing pipelines are particularly noticeable at 30 GHz; see figs. 42 and 50 of Planck Collaboration (2020g). Indeed, the time-dependent gain calibration at 30 GHz is difficult given the limited number of radiometers. This may lead to temperature-to-polarisation leakage in regions of the sky where the temperature fluctuations are largest, as set by the CMB dipole. Additionally, the upper right quadrant of the sky map coincides partially with effects stemming from analogue-to-digital converter non-linearities (ADCNL; PCP18) in the HFI data. Thus it would be interesting to see how our recovered spectral index maps might change when using alternative data products as discussed in Sec. 3.
corrNCM + corrPrior corrNCM + diagPrior diagNCM + diagPrior shows the fully correlated 'corrNCM + corrPrior' case, the middle column shows the 'corrNCM + diagPrior' case, and the right one shows the pixel-by-pixel 'diagNCM + diagPrior' case. Masked pixels are shown in grey. The corresponding histograms of the spectral index maps are shown in Fig. 2 and the associated means and standard deviations are given in Table 1.   Table 2. Summary of spectral index distributions for synchrotron and dust in six relevant radio loops, shown in Fig. 5, for full-mission SRoll1 and SRoll2 data. We compare three different cases: (i) corrNCM + corrPrior, (ii) corrNCM + diagPrior, and (iii) diagNCM + diagPrior. The spectral index maps are shown for SRoll2 in Fig. 3 and for SRoll1 in Fig. B1. spectral index maps, shown in Fig. 3, and the quantitative test in the loop regions, quoted in Table 2, reveal no clear evidence for dust spectral index variations across the sky.
Spectral index variations have been measured in e.g. Fuskeland et al. (2014) who with tight priors find evidence for weak spatial variations in polarised synchrotron emission using the nine-year Wilkinson Microwave Anisotropy Probe data (WMAP, Bennett et al. 2003). Their analysis, similar to the Gibbs sampling technique in Eriksen et al. (2008) and Dunkley et al. (2009a), finds an increasingly steep spectral index for higher Galactic latitudes, ranging from plane = −2.98 ± 0.01 to high−lat = −3.12 ± 0.04. Recent results on LFI Planck data by Svalheim et al. (2020) report similar findings. In appendix D we test for latitudinal spectral index variations but do not see a significant effect for either sync or dust .

Effect on cosmological parameter estimation
An important aim of the present investigation is to quantify the effect of spatially varying Galactic polarised foregrounds on cosmological parameter investigation. Therefore, we measure the optical depth to reionization parameter, , comparing our findings here using our Newton-Raphson component separation procedure (henceforth denoted by NR) to B21 which used simple template cleaning (henceforth denoted by TC). The TC method cleaned the SRoll2 signal channels at 100 and 143 GHz with uniformly scaled versions of 30 GHz (for synchrotron) and 353 GHz (for dust) half-mission maps using the coefficients of table 1 in B21. The NR analysis uses sets of either full-mission or half-mission maps. To explore the cosmological parameter space, we perform one-dimensional parameter scans in . Therefore, we generate theoretical power spectra with the base ΛCDM cosmology of (PCP18): 0 = 67.04, Ω b ℎ 2 = 0.0221, Ω c ℎ 2 = 0.12, Ω ℎ 2 = 0.00064, * = 1.0411, s = 0.96, and keeping 10 9 −2 = 1.870 fixed. A fiducial = 0.060 is assumed for forming the quadratic cross-spectra (QCS) and for the scans ranges from 0.01 to 0.2 with Δ = 0.001.
In Fig. 6 we compare the posteriors for from NR-cleaned SRoll2 data using the two likelihoods introduced in Sec. 2.2. We produce the two independent CMB maps required for the likelihoodapproximation scheme momento, which uses cross-spectra, by running the NR procedure twice, once on half-mission 1 (HM1) input maps and once on half-mission 2 (HM2) input maps. We label the obtained posterior HM1 × HM2. The pixel-based likelihood pixLike uses and full-mission maps to constrain , giving a posterior we label 'full, pixlike'. We use the appropriate block of the inverse of the recovered curvature matrix as the effective noise covariance matrix for each map in the likelihoods, propagating foreground-cleaning uncertainties into the estimates. We also show two posteriors obtained using TC cleaning of SRoll2 from B21. The 100 × 143 GHz full-mission posterior was the baseline re-  Figure 6. Posteriors for the optical depth to reionization using SRoll2 data for HFI. We compare the CMB maps obtained from the presented Bayesian component separation method (labelled NR) and template-cleaned (labelled TC) maps. We perform a one-dimensional parameter scan over using a pixel-based likelihood (pixLike) and a cross-spectrum-based likelihood (momento).
sult, and the 100 HM1×143 HM2 posterior relates to the NR-cleaned half-mission posterior (HM1 × HM2) in its use of half-mission data.
The means and standard deviations of the four posteriors are: It is reassuring to see that the two likelihoods give consistent results when using NR cleaning, and that both of these analyses are consistent with the TC-based results presented in B21. We stress again our use of wide priors in the NR method. Indeed, applying tight priors leads to results on that are higher than the TC-cleaned ones by several . Such results are expected as the freedom to reduce foreground residuals in the recovered CMB maps is constrained by a tight prior.

SUMMARY AND CONCLUSION
Foregrounds pose a major challenge to cosmological parameter estimation from large angular scale (2 ≤ ℓ ≤ 30) polarised CMB data. We have presented a Bayesian parametric foreground fitting framework to disentangle the primary CMB and the Galactic polarised foregrounds from one another. We modelled electron synchrotron radiation as a power law and thermal dust emission as a modified blackbody in frequency space. Our method allows for inter-pixel correlations in the noise covariance matrices of the input maps and for a physically motivated spatial correlation length in the spectral index prior matrices. We used six frequency channels from polarised Planck 2018 LFI and HFI data and found evidence for a spatially varying spectral index for synchrotron but not for dust. The means and dispersions of both spectral indices are given in Eqs. (13a) and (13b). From an analysis of our synchrotron spectral index maps, we find significantly steeper indices in the locations of radio loops. We have investigated whether the optical depth to reionization measured from the Planck 2018 LFI and SRoll2 HFI polarisation maps is sensitive to the method of polarized foreground removal. We have applied the Bayesian parametric foreground component separation method developed in this paper, using both a cross-spectrum-based likelihood momento and a pixel-based likelihood pixLike, to derive the posteriors given in Eqs. (14c)-(14b). These results are in excellent agreement with the simple uniform template cleaning method adopted in B21.
We do not anticipate any further methodological advances that would lead to substantive improvements in Planck polarization maps at large angular scales. Next-generation space-borne CMB experiments, however, will measure the sky with higher signal-to-noise and their analysis will be sensitive to the spatial variations of polarised foreground properties. Foreground cleaning and inference methods such as those presented here will thus be necessary for accurate cosmological inference from upcoming polarised CMB surveys with their ambitious science targets such as the measurement of primordial -modes. This will be especially acute at lower CMB frequencies in which synchrotron is a significant contaminant which has, as we have seen, a complex angular structure.

APPENDIX A: SPECTRAL INDEX MODELS
If radiation at a frequency has brightness , one defines a brightness temperature from this as where is the Boltzmann constant, and the speed of light. For a blackbody with = ℎ / , where ℎ is Planck's constant. Thus one can relate thermodynamic and brightness temperature variations through We model the synchrotron radiation as a power law and so take Δ sync B ∝ s . We characterise the dust spectral energy distribution (SED) by a modified blackbody where we assume d = 22.7 K  for the dust temperature. Therefore, the brightness temperature is where we have defined ( ) ≡ 1/( − 1) with = ℎ / d . The numerical values are given in Table A1, where we adopt a fiducial frequency for dust of d 0 = 100 GHz. The full foreground model in thermodynamic temperature units is given in Eq. (3).

APPENDIX B: SPECTRAL INDEX VARIATIONS IN SRoll1 DATA
In Fig. B1, we show the spectral index maps obtained in the three cases introduced in Sec. 2.1.2 using SRoll1 data in place of SRoll2 data for HFI. The corresponding histograms are shwon in Fig. B2.

APPENDIX C: MASK DEPENDENCE OF THE SRoll1 RESULTS
This appendix discusses the effect of increasing the fraction of unmasked sky pixel sky in the presented component separation method.  Figure B1. Distribution of synchrotron (upper panel) and dust (lower panel) spectral index maps using SRoll1 for HFI for the three inter-pixel correlation cases considered. Masked pixels are shown in grey. The means and standard deviations of the maps are given in Table 1. cause high amplitude pixels closer to the Galactic plane to affect areas away from the Galactic plane, potentially leading to convergence issues. Ignoring these correlations in the noise and the prior matrices, we test using a mask with sky fraction sky 0.66 (and a simple noise variance pattern suitable for this mask).
In Fig. C1 we show the recovered spectral index maps for sync and dust with the corresponding histograms shown in Fig. C2. The recovered mean and standard deviation for the diagNCM+diagPrior case with sky = 0.66 is sync = −3.055 ± 0.916 , and are consistent with the corresponding SRoll1 diag-NCM+diagPrior result presented in Table 1. We note that the re-βsync: diagNCM+diagPrior -5 -1 βdust: diagNCM+diagPrior 0.5 2.5 Figure C1. Distribution of synchrotron (left) and thermal dust (right) spectral indices across the sky using SRoll2 data for HFI. The masked pixels are shown in grey with a fraction of sky = 0.66 unmasked pixels. covered scatter in the spectral index maps is larger but attribute this to the simplified noise model used for this test.
Corresponding to the results discussed in Sec. 4.1, we find evidence for synchrotron spectral index variations in the radio loops I, III, IIIb, VIIb and IV, given in Table C1.  Figure C2. Histograms of spectral index maps shown in Fig. C1 using a mask with sky = 0.66. The input priors are shown as dashed vertical black lines for both spectral indices. Figure D1. Four stripes in Galactic latitude, labelled 1 to 4, used to investigate potential spectral index variation. Each stripe has a width of 13 • in latitude and they are centred at ±38 • (stripes 2 and 3) and ±70 • (stripes 1 and 4  Table D1. Summary of spectral index distributions for synchrotron and dust in four stripes in Galactic latitude, shown in Fig. D1, for SRoll2 data.

APPENDIX D: TESTING FOR LATITUDINAL SPECTRAL INDEX VARIATIONS
In this appendix, we test for spectral index variations as a function of Galactic latitude using SRoll2 data for HFI. We investigate four disjoint regions, shown in Fig. D1. Two 'stripes' are located in the northern and two stripes in the southern hemisphere, avoiding the Galactic plane and the Galactic poles. For the three pixel-correlation cases we consider, we obtain the results quoted in Table D1. We do not see any significant spectral index variations in either dust or synchrotron between these stripes. This paper has been typeset from a T E X/L A T E X file prepared by the author.