QUIJOTE scientific results - VIII. Diffuse polarized
foregrounds from component separation with
QUIJOTE-MFI

E. de la Hoz, R. B. Barreiro, P. Vielva, E. Martínez-González, J. A.

Rubiño-Martín, B. Casaponsa, F. Guidi, M. Ashdown, R. T. Génova-Santos,
E. Artal, et al.

► To cite this version:

E. de la Hoz, R. B. Barreiro, P. Vielva, E. Martínez-González, J. A. Rubiño-Martín, et al. QUIJOTE scientific results - VIII. Diffuse polarized foregrounds from component separation with QUIJOTE-MFI. Monthly Notices of the Royal Astronomical Society, 2023, 519, pp.3504-3525. 10.1093/mn-ras/stac3020. insu-03947303

HAL Id: insu-03947303

https://insu.hal.science/insu-03947303

Submitted on 10 Jul 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# QUIJOTE scientific results – VIII. Diffuse polarized foregrounds from component separation with QUIJOTE-MFI

E. de la Hoz<sup>1,2*</sup> R. B. Barreiro,<sup>1</sup> P. Vielva,<sup>1</sup> E. Martínez-González,<sup>1</sup> J. A. Rubiño-Martín<sup>1,3,4</sup>

B. Casaponsa,<sup>1</sup> F. Guidi<sup>1,3,4,5</sup> M. Ashdown,<sup>6,7</sup> R. T. Génova-Santos,<sup>3,4</sup> E. Artal,<sup>8</sup> F. J. Casas,<sup>1</sup>

R. Fernández-Cobos,<sup>9</sup> M. Fernández-Torreiro,<sup>3,4</sup> D. Herranz<sup>1</sup>,<sup>1</sup> R. J. Hoyland,<sup>3,4</sup> A. N. Lasenby,<sup>6,7</sup>

M. López-Caniego,<sup>10,11</sup> C. H. López-Caraballo,<sup>3,4</sup> M. W. Peel<sup>1,3,4</sup> L. Piccirillo,<sup>12</sup> F. Poidevin<sup>1,3,4</sup>

R. Rebolo,<sup>3,4,13</sup> B. Ruiz-Granados,<sup>3,4,14</sup> D. Tramonte,<sup>3,4,15,16</sup> F. Vansyngel<sup>1,3,4</sup> and R. A. Watson.<sup>12</sup>

*Affiliations are listed at the end of the paper*

Accepted 2022 October 14. Received 2022 October 14; in original form 2022 July 27

## ABSTRACT

We derive linearly polarized astrophysical component maps in the Northern Sky from the QUIJOTE-MFI data at 11 and 13 GHz in combination with the Wilkinson Microwave Anisotropy Probe K and Ka bands (23 and 33 GHz) and all Planck polarized channels (30–353 GHz), using the parametric component separation method B-SECRET. The addition of QUIJOTE-MFI data significantly improves the parameter estimation of the low-frequency foregrounds, especially the estimation of the synchrotron spectral index, $\beta_s$. We present the first detailed $\beta_s$ map of the Northern Celestial Hemisphere at a smoothing scale of 2°. We find statistically significant spatial variability across the sky. We obtain an average value of –3.08 and a dispersion of 0.13, considering only pixels with reliable goodness of fit. The power-law model of the synchrotron emission provides a good fit to the data outside the Galactic plane but fails to track the complexity within this region. Moreover, when we assume a synchrotron model with uniform curvature, $c_s$, we find a value of $c_s = –0.0797 \pm 0.0012$. However, there is insufficient statistical significance to determine which model is favoured, either the power law or the power law with uniform curvature. Furthermore, we estimate the thermal dust spectral parameters in polarization. Our cosmic microwave background, synchrotron, and thermal dust maps are highly correlated with the corresponding products of the PR4 Planck release, although some large-scale differences are observed in the synchrotron emission. Finally, we find that the $\beta_s$ estimation in the high signal-to-noise synchrotron emission areas is prior-independent, while, outside these regions, the prior governs the $\beta_s$ estimation.

**Key words:** cosmology: observations – methods: data analysis – polarization – cosmic microwave background.

## 1 INTRODUCTION

Currently, most of the efforts of the cosmic microwave background (CMB) community are devoted to the search for primordial B modes. These predicted B modes at large scales can only be produced by tensor modes, and their detection would constitute compelling evidence of an inflationary phase. The intensity of this primordial signal is determined by the tensor-to-scalar ratio $r$, the relative amplitude between the tensor and scalar modes at a given pivot scale. The current best upper bound on the tensor-to-scalar ratio is $r < 0.032$ at 95 per cent CL, set by the combination of Planck, BICEP2/KeckArray, and baryon-acoustic-oscillation data (Tristram et al. 2022).

The weakness of the primordial B modes makes its detection a tremendous experimental challenge, requiring high-sensitivity experiments as well as an exquisite control of systematics. Indeed, a large effort is currently on-going with the aim to detect, or at least to

constrain, $r$ with a sensitivity $\sigma_r(r = 0) \le 10^{-3}$. This includes many planned ground-based experiments, e.g. GroundBIRD (Lee et al. 2020), LSPE-Strip (Lamagna et al. 2020), CMB-S4 (Abazajian et al. 2016), Simons Observatory (Ade et al. 2019), and BICEP array (Hui et al. 2018), as well as satellite missions [e.g. LiteBIRD (LiteBIRD Collaboration et al. 2022) and PICO (Hanany et al. 2019)].

The detectability of the primordial B modes could be improved by removing the secondary B-mode component induced by weak gravitational lensing. Several delensing procedures have been proposed in the literature (Planck Collaboration 2016b; Millea, Anderes & Wandelt 2019) and have been applied to data from current CMB experiments (Planck Collaboration 2016b; Carron, Lewis & Challinor 2017; BICEP/Keck Collaboration 2021), and in forecasts of future CMB experiments (Diego-Palazuelos et al. 2020; Namikawa et al. 2022).

It is necessary to disentangle the CMB polarization signal from those coming from other microwave emissions, such as Galactic synchrotron, thermal dust, and extragalactic point sources. Thus, the problem of component separation is a crucial step in order to detect the primordial B mode of CMB polarization. This process

<footnote>*E-mail: delahoz@ifca.unican.es</footnote>
benefits from the characterization of foreground emissions, using complementary frequency ranges that provide unique information about the contaminants.

The main diffuse polarized contaminants are the synchrotron emission (at low frequencies) and the thermal dust emission (at high frequencies). The best characterization of these diffuse foregrounds has been done by Planck (Planck Collaboration IV 2020d), using a data set covering frequencies from 30 to 353 GHz. This frequency range limited strongly the estimation of the synchrotron spectral parameters. In Planck Collaboration IV (2020d), it is shown that, with Planck data only, one cannot test the spatial variability of the synchrotron spectral index due to limited sensitivity and frequency coverage. The data only allows a measurement of a global spectral index of $\beta_s = -3.1 \pm 0.1$. The synchrotron spectral index has also been estimated using other data sets (e.g. Fuskeland et al. 2014; Krachmalnicoff et al. 2018; Fuskeland et al. 2021).

The Q-U-I JOint Tenerife Experiment (QUIJOTE; Rubiño-Martín et al. 2010) is a polarimetric ground-based CMB experiment whose main scientific goal is the characterization of the polarization of the CMB and other Galactic and extragalactic physical processes in the frequency range of 10–40 GHz and at large angular scales ($\gtrsim 1^\circ$). The experiment is located at the Teide Observatory (at $\sim 2400$ m above sea level) in Tenerife. It is composed of two telescopes equipped with three instruments: the Multi-Frequency Instrument (MFI), the Thirty-GHz Instrument (TGI), and the Forty-GHz Instrument (FGI), operating at 10–20, 26–36, and 39–49 GHz, respectively.

The MFI instrument has been operating from 2012 November to 2018 October. It conducted two different surveys: (i) a shallow Galactic survey (called ‘wide survey’) covering all the visible sky from Tenerife at elevations larger than 30°, and (ii) a deep cosmological survey covering approximately 3000 deg² in three separated sky patches in the northern sky. In this work we use the QUIJOTE-MFI wide survey maps. This survey provides an average sensitivity in polarization of $\sim 35–40$ µK per 1-deg beam in four bands centred around 11, 13, 17, and 19 GHz (Rubiño-Martín et al. 2022). Those frequencies are crucial to achieving a better characterization of the low-frequency foregrounds. In intensity, this additional information helps breaking degeneracies between the synchrotron, free-free, and anomalous microwave emissions (AMEs) while, in polarization, the QUIJOTE-MFI channels are key to characterize the synchrotron spectral dependence.

In this work, we perform a component separation analysis to obtain more information about the polarized sky using the QUIJOTE-MFI data<fnref n="1" /> (Rubiño-Martín et al. 2022) in combination with the publicly available Planck (Planck Collaboration I 2020a; Planck Collaboration LVII 2020f) and Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2013) data. To perform component separation analysis, we use B-SeCRET (Bayesian-Separation of Components and Residual Estimation Tool), a parametric maximum-likelihood method described in de la Hoz et al. (2020).

The paper is organized as follows. In Section 2, we provide details of the main components in the polarized microwave sky and the corresponding parametric models used to characterize them. Section 3 describes briefly the B-SeCRET method. The data used in the analysis are presented in Section 4. Then, the main component separation results obtained are shown in Section 5. Finally, the main conclusions from the analysis are given in Section 6. In Appendix A, we provide maps of the synchrotron spectral index obtained from independent fits in linear Stokes parameters Q and U. Appendix B

compares the variations on the synchrotron spectral index due to rotations of the polarized angle with Faraday rotation.

## 2 THE MICROWAVE SKY MODEL

The polarized microwave sky is composed primarily of photons from the CMB, synchrotron, and thermal dust. As stated before, the synchrotron emission dominates at low-frequencies while the thermal dust is the principal component at higher frequencies. The contribution from other components, discussed in Section 2.5, is expected to be insignificant and not taken into account. Apart from these astronomical signals, the measured sky signal maps have another contribution from the instrumental noise. The characteristics of this noise depend on the specifications of the experiment. Furthermore, contaminants such as the atmosphere and artificial signals from satellites also contribute to the microwave sky, see Rubiño-Martín et al. (2022) for more details. Thus, the measured polarized sky signal for a given v channel can be expressed as the following sum:

$$ \left(\begin{array}{c} Q \\ U \end{array}\right)_{\nu} = \left(\begin{array}{c} Q_{\text{cmb}} \\ U_{\text{cmb}} \end{array}\right)_{\nu} + \left(\begin{array}{c} Q_{s} \\ U_{s} \end{array}\right)_{\nu} + \left(\begin{array}{c} Q_{d} \\ U_{d} \end{array}\right)_{\nu} + \left(\begin{array}{c} Q_{n} \\ U_{n} \end{array}\right)_{\nu}, \quad (1) $$

where $X_{cmb}$, $X_s$, and $X_d$ are the CMB, synchrotron, and thermal dust signals, respectively, and $X_n$ is the instrumental noise ($X \in \{Q, U\}$). In the subsequent subsections, we describe the main physical components that encompass the sky signal as well as some effects that alter this signal. Moreover, we present the parametric models that we use in the component separation analysis for each polarized astronomical component.

### 2.1 Synchrotron

The synchrotron emission arises from relativistic particles (cosmic rays) passing through the Galactic magnetic field. Its emissivity depends both on the magnetic field strength and energy distribution of the relativistic particles (generally electrons). These quantities are not uniform in the Galactic disc. For instance, the free electrons are more predominant in compact regions as supernovae remnants. On the other hand, the magnetic field is amplified in some compact regions and can have different strength and direction across the sky.

The synchrotron spectral energy distribution (SED) is generally described as a power law (Rybicki & Lightman 2008):

$$ \left(\begin{array}{c} Q_s \\ U_s \end{array}\right)_{\nu} = \left(\begin{array}{c} A_s^Q \\ A_s^U \end{array}\right)_{\nu} \left(\frac{\nu}{\nu_s}\right)^{\beta_s}, \quad (2) $$

where $A_s$ is the amplitude in brightness temperature at the pivot frequency $\nu_s = 30$ GHz and $\beta_s$ is the spectral index which is assumed to be equal for both Q and U Stokes parameters.

Previous works dedicated to the estimation of the spectral index, found values around $\beta_s \simeq -3.1$ (Planck Collaboration IV 2020d). However, the spectral index is expected to vary spatially due to its dependence on the energy distribution of the cosmic rays $N(E)$. Studies such as Fuskeland et al. (2014), Vidal et al. (2015), Krachmalnicoff et al. (2018), Martire, Barreiro & Martínez-González (2022), and Weiland et al. (2022) indicate that different polarized regions present different spectral indices. Here, we conduct a more detailed analysis of the $\beta_s$ spatial variations in the Northern Hemisphere by performing a pixel-by-pixel component separation analysis using the QUIJOTE MFI polarized maps.

The S-PASS survey (Carretti et al. 2019) has provided the most sensitive reconstruction of the $\beta_s$ variations of the South Celestial Hemisphere (Krachmalnicoff et al. 2018). They found large variability over the sky, and a mean value of $-3.22 \pm 0.08$. Those results

<footnote n="1">¹This is one of the papers which are part of the MFI wide survey data release.</footnote>
were further confirmed in the analysis of Fuskeland et al. (2021) that estimated the spectral index taking into account the Faraday Rotation effect. They also studied the Galactic plane and found compatible results to those where only WMAP data were used, finding a flatter index in the Galactic plane than at high Galactic latitudes.

We have also considered an extension of equation (2), where we include a possible curvature in the synchrotron’s SED:

$$ \left( \frac{Q_s}{U_s} \right)_v = \left( \frac{A_s^Q}{A_s^U} \right)_v \left( \frac{\nu}{\nu_s} \right)^{\beta_s + c_s \log\left(\frac{\nu}{\nu_s}\right)}, \quad (3) $$

where $c_s$ is the parameter that represents the curvature. This extension is worth studying since a curved spectrum can account for steepening or flattening of the SED due to different effects, e.g. cosmic ray aging effect or multiple synchrotron components along the line of sight. This model could also account for the presence of polarized AME.

## 2.2 Thermal dust

The thermal dust radiation comes from dust grains present in the interstellar medium. Those grains absorb ultraviolet light and re-emit as a grey body. In general, these dust grains are not perfectly spherical and typically have their minor axis aligned with the direction of the local magnetic field. This effect yields polarized thermal dust emission. The SED of this radiation is often described as a modified blackbody with emissivity index $\beta_d$ and dust temperature $T_d$:

$$ \left( \frac{Q_d}{U_d} \right)_v = \left( \frac{A_d^Q}{A_d^U} \right) \left( \frac{\nu}{\nu_d} \right)^{\beta_d+1} \frac{e^{\gamma \nu_d} - 1}{e^{\gamma \nu} - 1}, \quad (4) $$

where $A_d$ is the amplitude of the dust in brightness temperature evaluated at the pivot frequency $\nu_d = 143$ GHz and $\gamma = \frac{h}{k_B T_d}$. The amplitude is well characterized by the higher frequency channels, where the other components are clearly subdominant. The current temperature map of the dust grains ($T_d$) is obtained from temperature analysis and has values mostly between 14 and 26 K. The polarized dust emissivity evaluated with Planck data is $\beta_d = 1.55 \pm 0.05$ (Planck Collaboration IV 2020d).

Several works support the idea that a single component dust model is too simplistic and more components might be required to fully characterize this emission (e.g. McBride, Bull & Hensley 2022; Ritacco et al. 2022). None the less, since this paper is focused on the low frequency foregrounds, we keep the model used in Planck Collaboration VI (2020d) which seems to provide a good description at the Planck polarized frequencies ($30 \text{ GHz} < \nu < 353 \text{ GHz}$).

## 2.3 CMB

The CMB radiation has a thermal blackbody spectrum with a temperature of $T_o = 2.7255 \pm 0.0006$ K (Fixsen 2009). CMB photons are linearly polarized due to the Thomson scattering experienced with the hot electron gas at the last scattering surface. Unlike in intensity, where the CMB can be the dominant contribution at intermediate frequencies (70–150 GHz) and high Galactic latitudes, in polarization, the foreground contribution cannot be overlooked. Thus, in order to detect the primordial B-mode, experiments with very high sensitivity, exquisite control of systematics and a careful removal of foregrounds are mandatory.

The CMB signal at each pixel is given by its amplitude $A_{cmb}$, which is the only free parameter for this component. Since the rest

of the components are given in brightness temperature, we convert the thermodynamic temperature of the CMB to the same units:

$$ \left( \frac{Q_{cmb}}{U_{cmb}} \right)_v = \left( \frac{A_{cmb}^Q}{A_{cmb}^U} \right) \frac{x^2 e^\nu}{(e^\nu - 1)^2}, \quad (5) $$

where $x = \frac{h\nu}{k_B T_o}$.

## 2.4 Faraday rotation

Another issue intrinsic to the polarization signal is the Faraday rotation effect, i.e. the rotation of the plane of polarization that occurs when light passes through the interstellar medium in the presence of a magnetic field. The magnitude of this effect scales with the square of the wavelength; hence, its repercussions are more significant at low frequencies. To properly account for this effect we require a broad knowledge of the Galactic magnetic field as well as the interstellar medium, in order to recognize the regions where the effect is more significant. Moreover, since the instrumental beam has a finite size, the measured signal is an average of the emission from various directions within the beam with slightly different rotation angles. This results in a ‘beam depolarization’ of the signal.

Hutschenreuter et al. (2022) show that the possible Faraday Rotation effects at the QUIJOTE-MFI frequencies (10–20 GHz) are very small in most of the sky, and particularly at high Galactic latitudes. Thus, in this work we have not considered any Faraday Rotation effect. Nevertheless, in Appendix B we study variations on the synchrotron spectral index due to rotations of the polarized angle and compare it to Faraday Rotation models such as the one proposed in Hutschenreuter et al. (2022).

## 2.5 Other contributions

It is well known that there are other foreground components whose emissions are important for intensity analyses. In particular, at low frequencies, one needs to consider two additional Galactic emission components: the bremsstrahlung radiation generated from electron ion scattering in interstellar plasma (free-free), and AME, whose physical origin still is not fully clear. At high frequencies, in addition to thermal dust, we find an isotropic extragalactic emission called the cosmic infrared background (CIB), coming from different sources, e.g. dusty star-forming galaxies, quasars, intergalactic stars, and intercluster dust in the Local group. We also have other contributions such as CO line emission or Sunyaev-Zeldovich effect (SZ) from clusters of galaxies (Sunyaev & Zeldovich 1972) that should be taken into account in intensity analyses (Planck Collaboration 2016a). In addition, emission from extragalactic point sources, both at radio and infrared frequencies is an important contaminant at small scales. In polarization the problem is simplified since several of these emissions (free-free, CIB, SZ) are not expected to be polarized (at least significantly); therefore, we do not consider them.

The polarization of the AME is still under study because its nature is still uncertain (Dickinson et al. 2018). Several models have been proposed such as spinning dust particles (Ali-Haïmoud 2013), magnetic dipole emission (Draine & Lazarian 1999), or more recently the proposal of spinning nano-diamonds (Greaves et al. 2018). The predicted polarization fraction of the AME emission for most of these models is below 5 per cent. From the data analysis point of view, no evidence of polarization has been found in compact region studies [the most stringent constraints on the polarization fraction, $\Pi$, have been provided by Génova-Santos et al. (2017), $\Pi < 0.22$ per cent at 41 GHz]. Due to this lack of evidence, we do not take into account the AME component in this work.

<footnote>$^2h$ and $k_B$ are Planck and Boltzmann constants, respectively.</footnote>
On the other hand, point sources present some degree of polarization, which is in general small (a few percent). However, at the resolutions considered in this work, they are subdominant with respect to Galactic foregrounds. Thus, we do not include them in our analysis. We note however that in the data, a few polarized point sources are present that are not taken into account in the component separation analysis (see Herranz et al. 2022).

# 3 COMPONENT SEPARATION METHODOLOGY

In this work, we apply the parametric component separation method B-SeCRET to extract the polarized astrophysical signals. Parametric methods are very powerful since they provide physical information of each sky component. However, they require a profound theoretical understanding of the nature of the foregrounds and accurate knowledge of the experiment's characteristics to avoid biases in the analysis.

Below, in Section 3.1, we outline the component separation technique applied in this work. Then, in Section 3.2, we describe the prior information that is used in the Bayesian analyses.

## 3.1 Bayesian analyses

The B-SeCRET methodology is a parametric pixel-based maximum-likelihood method, which relies on an Affine-Invariant Markov Chain Monte Carlo Ensemble sampler to draw samples from a posterior distribution (Foreman-Mackey et al. 2013). This methodology has already been applied in previous studies (e.g. de la Hoz et al. 2020, 2022).

B-SeCRET applies Bayesian inference to determine the best-fitting model parameters given some prior information. In Bayesian statistics, the probability of the set of model parameters $\theta_p$ given the signal data $d_p$ at the pixel $p$ is proportional to the probability of the $d_p$ given $\theta_p$ times the probability of $\theta_p$, i.e.

$$ P(\theta_p | d_p) \propto P(d_p | \theta_p) P(\theta_p). \quad (6) $$

$P(\theta_p)$ is commonly known as the prior information, whereas $P(d_p|\theta_p)$ is usually referred to as the likelihood. Assuming Gaussian noise, the likelihood of the data can be expressed as

$$ P(d_p|\theta_p) = \frac{\exp\left(-\frac{1}{2}(d_p - S_p)^T \mathbf{C}^{-1} (d_p - S_p)\right)}{\sqrt{(2\pi)^N \det(\mathbf{C})}}, \quad (7) $$

where $\mathbf{C}$ is the noise covariance matrix, $N$ is the number of elements in the $d_p$ array, and $S_p$ is the parametric model considered, which has been described in detail in Section 2.

To draw samples from the posterior probability, we use the PYTHON implementation EMCEE (Foreman-Mackey et al. 2013) of an affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC; Goodman & Weare 2010). In each pixel, the best-fitting parameters and their uncertainties are obtained as the median and the standard deviation of their respective marginalized posterior probability.

## 3.2 Priors

In this work, we benefit from prior information about astrophysical foregrounds to help with convergence and computational time reduction. For example, the synchrotron spectral index is known to be around -3.1, although experiments such as S-PASS found a more negative value. Here, we use the estimated value obtained with Planck polarization data by the SMICA method, $\beta_s = -3.1 \pm 0.06$ (Planck

<figure><img src="image_2.png" /><figcaption>**Figure 1.** QUIJOTE observed sky after removing the geostationary satellite band and the region around the north celestial pole, which is affected by high atmospheric air-mass [$f_{sky} = 51$ per cent, Galactic coordinates centred on (0,0)].</figcaption></figure>

Collaboration IV 2020d) and use a broad Gaussian distribution $N(-3.1, 0.3)^3$ as a prior on $\beta_s$. When we include a curvature in the synchrotron model we apply a Gaussian prior $N(0, 0.1)$ on $c_s$. Moreover, we apply Gaussian priors $N(1.55, 0.1)$ and $N(21, 3)$ on both $\beta_d$ and $T_d$, respectively. Finally, flat priors are used in the characterization of the amplitude parameters.

# 4 DATA

The aim of this work is to obtain a better characterization of the low-frequency foregrounds by including the newly released QUIJOTE-MFI maps in component separation analyses. In this section, we summarize the basic details of these maps as well as those from the other experiments used in the analysis, i.e. the $K$ and $Ka$ bands from WMAP and Planck's third and fourth public releases (PR3 and PR4, respectively). We also discuss some technical issues related to the instruments such as the estimated noise, RFI, and the colour corrections.

## 4.1 Data sets

In this analysis, we have used the data from the following experiments:

(i) QUIJOTE: We have used the low frequency QUIJOTE MFI 11 and 13 GHz channels (MFI) (Rubiño-Martín et al. 2022) due to their better signal-to-noise ratio. Although QUIJOTE has observed 70 per cent of the sky there are regions with poorer sensitivity due to the presence of artificial satellites and high atmospheric masses in some directions. Thus, in this analysis we have considered the mask shown in Fig. 1, as the observable sky. This mask (satband + NCP) is described in Rubiño-Martín et al. (2022).

(ii) WMAP: We have used the low-frequency Nine-Year WMAP $K$ (22.8 GHz) and $Ka$ (33.1 GHz) bands (Bennett et al. 2013).$^4$

(iii) Planck: We have used the full set of Planck polarization maps, i.e. the low-frequency instrument (LFI) 30, 44, and 70 GHz frequency maps and the high-frequency instrument (HFI) 100, 143, 217, and 353 GHz maps. We have obtained results from both PR3$^5$ (Planck

<footnote>$^3$N($x, \sigma$) represents a normal distribution with mean $x$ and variance $\sigma^2$.</footnote>

<footnote>$^4$The other bands were not included since they have a much lower synchrotron signal-to-noise ratio and do not contribute to the determination of the synchrotron characteristics.</footnote>

<footnote>$^5$We used the Planck maps corrected from bandpass leakage.</footnote>
Collaboration II 2020b; Planck Collaboration III 2020c) and PR4 (Planck Collaboration LVII 2020f) data releases.

Before component separation analyses, the frequency maps are all convolved (taking appropriately into account the beam window function of each particular frequency map) with a common beam, a Gaussian beam of full width at half-maximum (FWHM) = 2°, and downgraded to the same resolution through spherical harmonics, given by the HEALPix parameter $N_{side} = 64$. The procedure followed is described below:

(i) We calculate the spherical harmonics coefficients ($t_{\ell m}$, $e_{\ell m}$, $b_{\ell m}$) using the healpy routine map2alm.

(ii) To convolve all channels with the same beam we multiply the ($t_{\ell m}$, $e_{\ell m}$, $b_{\ell m}$) by $b_{\ell}(2^\ell) p_\ell(64)/(b_{i,\ell} p_\ell(N_{side}))$, where $b_{\ell}(\alpha)$ is a Gaussian beam window function whose FWHM is $\alpha$, $b_{i,\ell}$ is the $i$th channel beam window function, and $p_\ell(N_{side})$ is the pixel window function at the resolution $N_{side}$.

(iii) We obtain the downgraded maps at $N_{side} = 64$, applying the healpy routine alm2map to the new ($t_{\ell m}$, $e_{\ell m}$, $b_{\ell m}$).

Several combinations of the previous data sets have been tested. Each configuration's name is given by the 'sum' of the sets of maps included in the analysis. For example, the configuration composed of PR4 channels in combination with WMAP's K and Ka bands is referred as K/Ka+PR4, or MFI-QUIJOTE low-frequency channels in combination with PR4 and WMAP channels is specified as MFI+K/Ka+PR4.

## 4.2 Instrumental effects

Real data present different instrumental effects that need to be accounted for. For example, an important contribution to the observed signal is the noise produced by the detectors of each experiment. A proper characterization of the noise levels is key for component separation analyses. In this work, we have calculated the covariance matrix among the frequency channels per pixel, required by the parametric component separation method, using realistic noise simulations specific to each instrument. Each experiment's noise simulations are obtained as follows:

(i) QUIJOTE: We have used the correlated noise simulations described in Rubiño-Martín et al. (2022). They account for the 1/f noise present in the maps, and the correlated noise component between 11 and 13 GHz.

(ii) WMAP: We have generated a set of white noise simulations using the RMS noise per pixel provided by the WMAP collaboration (Hinshaw et al. 2003). The RMS noise $\sigma$ is calculated as $\sigma = \sigma_0/\sqrt{N_{obs}}$.<fnref n="6" />

(iii) Planck: For PR3, we have used the FFP10 simulations generated by the Planck Collaboration (Planck Collaboration I 2020a). In the case of the PR4, we have employed the noise simulations described in Planck Collaboration LVII (2020f).<fnref n="7" />

While the frequency channels of different experiments are uncorrelated, there might be correlations between channels of a given instrument. This is the case for the 11 and 13 GHz low-frequency MFI channels. On the other hand, we have assumed no correlations between frequency channels for WMAP and Planck. Thus, for a

given pixel $p$, the Planck and WMAP frequency covariance matrices are diagonal while QUIJOTE's has non-zero off-diagonal terms. For a given configuration, the covariance matrix is obtained as a block matrix, where each block corresponds to the frequency covariance matrix of each instrument included in that configuration.

To obtain the experiments' frequency covariance matrices, first we pre-process the noise simulations in the same manner as the data maps. Then, for Planck and WMAP, the diagonal terms are calculated as the variance of the noise simulations at the corresponding pixel for each frequency. Each pixel covariance matrix between QUIJOTE 11 and 13 GHz is calculated as the sample covariance matrix using the values of the 11 and 13 GHz noise simulations at that specific pixel.

One test to verify that our covariance matrices are well estimated is the following. We obtained a distribution of $\chi^2_{n,i}$ values as

$$ \chi^2_{n,i} = n_i^T C_i^{-1} n_i, \qquad (8) $$

where $n_i$ is a noise simulated map<fnref n="8" /> at the frequency $i$ and $C_i$ is the noise covariance matrix described above. The $\chi^2_{n,i}$ distributions should have the expected form with $N_{pix}$ degrees of freedom (dof). This is consistent with the values obtained for Planck and WMAP. In the case of QUIJOTE, the distribution deviates slightly from the expected $N_{pix}$ dof $\chi^2$-distribution, since they are not end-to-end noise simulations and hence not as accurate (see Rubiño-Martín et al. 2022 for details). However, as subsequent analyses will show, we find that when the astrophysical emission is included, the obtained $\chi^2$ is correct as expected, i.e. in the regions where the model properly explains the data (outside the Galactic plane). Thus, QUIJOTE's noise simulations are accurate enough to perform scientific analyses.

We explored the possibility of including correlations among neighbouring pixels within a 1 deg radius.<fnref n="9" /> The smoothing process of the maps induces noise correlations among different pixels and, although this does not affect our pixel-by-pixel analyses, it can affect analyses where we assume a uniform parameter value within one region. Therefore, for each pixel, we calculated the covariance matrix among its neighbouring pixels from noise simulations. Then we generated a sparse covariance matrix where the only non-zero values in each row were the diagonal element and the correlation with the neighbouring pixels. In this case, the distribution does not follow a $N_{pix}$ degrees-of-freedom $\chi^2$ distribution as one would expect. The recovered values were smaller than expected, more notably for Planck maps. This is a consequence of not having enough noise simulations, which prevents us from obtaining a good characterization of the noise correlations. Therefore, we use the covariance matrices that do not take into account possible noise correlations among neighbouring pixels in the following.

As explained in Rubiño-Martín et al. (2022), in order to correct residual RFI signals emerging after co-adding all data in the mapping process of the QUIJOTE-MFI data, the polarization maps are corrected using a function of the declination (FDEC). This correction is equivalent to applying a filter to QUIJOTE data, which removes the zero mode in lines of constant declination. In Appendix C, we studied whether this correction affects the recovery of foregrounds spectral parameters such as $\beta_s$. We found that if only QUIJOTE maps are filtered with FDEC the recovered $\beta_s$ map is biased in regions such as the North Polar Spur. When all data maps are

<footnote n="6">⁶$\sigma_0$ and $N_{obs}$ are given in https://lambda.gsfc.nasa.gov/product/wmap/dr5/skymap_info.html.</footnote>

<footnote n="7">⁷Simulations available at NERSC under /global/cfs/cdirs/cmb/data/planck2020.</footnote>

<footnote n="8">⁸The noise simulations used in this test are different from the noise simulations used to calculate the noise covariance matrices.</footnote>

<footnote n="9">⁹The pixels contained within this radius are the ones with the largest correlations induced by the smoothing process.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. Synchrotron spectral index (top row) and uncertainty maps (middle row) obtained after component separation with four different data sets. The synchrotron emission is modelled with a power law. Bottom row: Reduced $\chi^2$ map obtained for each data set.</figcaption></figure>

filtered in the same way this bias disappears. Thus, for this analysis we have filtered all signal maps with their corresponding FDEC function.

Another important instrumental effect arises from detectors having a finite bandwidth. This issue has to be taken into account when dealing with foreground components whose amplitude varies within that frequency band. This effect can be corrected by adding a multiplicative factor, called colour correction, to the signal that depends on the spectral behaviour. We have used the fastcc PYTHON code (Peel et al. in preparation; Genova-Santos et al. in preparation) to obtain the colour corrections of each experiment considered here. Therefore, our model for the sky signal presented in Section 2 is corrected as follows:

$$X_{\nu} = X_{\nu, \text{cmb}} + \frac{X_{\nu, \text{s}}}{C_{\text{s}}(\alpha, \nu)} + \frac{X_{\nu, \text{d}}}{C_{\text{d}}(\beta_{\text{d}}, T_{\text{d}}, \nu)}, \quad (9)$$

where X is either Q or U, $C_s(\alpha, \nu)$ is synchrotron colour correction whose spectral behaviour is modelled as a power law with $\alpha = \beta_s + 2$. The spectral behaviour of dust colour correction $C_d(\beta_d, T_d, \nu)$ is assumed to be a modified blackbody and it is determined by its $\beta_d$ and $T_d$ parameters. The colour correction values are updated in each MCMC iteration.

# 5 RESULTS

In this section, we present the component separation products obtained using the recently released QUIJOTE low-MFI data along with the already available Planck and WMAP data. We have focused primarily on the synchrotron spectral parameters since those are the parameters where a greater improvement is found (see Sections 5.1 and 5.2). Moreover, we show the recovered amplitudes of the CMB, synchrotron and thermal dust and, compare them with those obtained by Commander using PR4 data in Section 5.3. In Section 5.4, we present the spectral parameters of the thermal dust. Finally, we evaluate the robustness of these results in Section 5.5.

## 5.1 Synchrotron spectral index

The major improvement obtained from including the low-frequency QUIJOTE-MFI channels is having the sufficient sensitivity to study the synchrotron spectral index with great accuracy. Here, we have conducted a deep study on several aspects with regard to this parameter. First, we have compared the recovered $\beta_s$ maps using different combinations of the available data sets (Section 5.1.1). Section 5.1.2 studies the spatial variability of $\beta_s$. Finally, we compare our results to the available $\beta_s$ models that are often exploited in simulations used in CMB science forecasts in Section 5.1.3.

### 5.1.1 Data sets

We have obtained different $\beta_s$ maps from component separation analyses using the four following data sets: WMAP K and Ka bands with PR4 (K/Ka+PR4); QUIJOTE-MFI 11 and 13 GHz channels with PR4 (MFI+PR4); QUIJOTE-MFI 11 and 13 GHz channels, WMAP K and Ka bands and PR4 (MFI+K/Ka+PR4) and QUIJOTE-MFI 11 and 13 GHz channels, WMAP K and Ka bands and PR3 (MFI+K/Ka+PR3). The results are shown in Fig. 2. It is clear from the comparison of the synchrotron spectral index uncertainty maps obtained in the K/Ka+PR4 case (first column) with respect to the MFI+K/Ka+PR4 case (third column), that the inclusion of QUIJOTE channels significantly improves the estimation of $\beta_s$. Moreover, we observe that, outside the Galactic plane, the estimation of $\beta_s$ is very close to the mean value of the prior set on this parameter, in this case -3.1. In other words, the information contained in that fraction of the data, i.e. the likelihood, is very poor and the estimation is driven by the prior.

This improvement does not come from the inclusion of more channels, but from channels where the synchrotron contribution is larger. This is evident from the comparison of the results from K/Ka+PR4 with respect to MFI + PR4, where the number of frequency channels is the same but the results are significantly better for the latter.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Reduced $\chi^2$, $\chi^2_{red}$, obtained using the MFI+$K/Ka$+PR4 data set versus the $\chi^2_{red}$ obtained using $K/Ka$+PR4 (left), MFI+PR4 (centre), and MFI+$K/Ka$+PR3 (right). The colour scale is related to the density of points; redder (bluer) corresponds to denser (sparser) regions. The orange rectangle shows the $\chi^2_{red}$ within a 95 per cent confidence region. The slope calculated with the points within this 95 per cent confidence region is $m = 0.686 \pm 0.004$ (left column), $m = 0.732 \pm 0.003$ (centre column), and $m = 0.731 \pm 0.003$ (right column), shown with a green dashed line. The orange dashed line shows the one-to-one line. The synchrotron emission is modelled with a power law.</figcaption></figure>

Finally, we have compared also the results obtained with MFI+$K/Ka$+PR3 and MFI+$K/Ka$+PR4 (fourth and third columns, respectively). In this case, the recovered uncertainty maps are virtually the same but there are some distinct differences between the $\beta_s$ maps that should be ascribed to changes in Planck maps.

One of the advantages of using a parametric component separation method is that we can evaluate the goodness of the fit with certain estimators. In this work, we use the reduced $\chi^2$ estimator, whose value at a given pixel $p$ is calculated as

$$ \chi^2_{\text{red},p} = \frac{1}{N_{\text{dof}}} \sum_{i \in \{Q, U\}} (\mathbf{d}_{p,i} - \mathbf{S}_{p,i}) \mathbf{C}_{p,i}^{-1} (\mathbf{d}_{p,i} - \mathbf{S}_{p,i}), \quad (10) $$

where the sum is over all Q and U frequency channels, and $N_{\text{dof}}$ is the number of dof. The bottom row of Fig. 2 shows the $\chi^2_{red}$ maps obtained for each data set combination. These maps show that our default model, i.e. a power law and a modified blackbody to model the synchrotron and thermal dust emission, respectively, provides a good fit (low values of $\chi^2_{red}$) outside the Galactic plane. Within the Galactic plane, this model is not able to capture all the physical complexity and the $\chi^2_{red}$ values are quite large. However, we note that in this analysis we have considered statistical uncertainties but not calibration errors, which in QUIJOTE are of 5 per cent. Apart from the higher complexity of the Galactic plane emission, the higher $\chi^2_{red}$ in this region could also be due, in part, to having neglected calibration errors.

We have also used the $\chi^2_{red}$ estimator to select the data set that is used as the default for further tests between the MFI+$K/Ka$+PR3 and the MFI+$K/Ka$+PR4 data sets, i.e. the only combinations that include all the channels considered. In Fig. 3, the $\chi^2_{red}$ obtained using the MFI+$K/Ka$+PR4 data set is plotted against the $\chi^2_{red}$ obtained with $K/Ka$+PR4, MFI+PR4, and MFI+$K/Ka$+PR3. The 95 per cent confidence regions are delimited by orange lines. These lines indicate the $\chi^2$ values, from the reduced $\chi^2$-distribution with $N_{\text{dof}}$ dof,<fnref n="10" /> that satisfy that the normalized area covered to their left is equal to 0.95. We have also fitted the points within this confidence regions to a straight line to determine which data set has more pixels with smaller $\chi^2_{red}$. If the slope is larger than unity, the data set on the horizontal axis has more pixels with smaller $\chi^2$. If the slope is smaller than

unity, the data set on the vertical axis is the one which satisfies that condition.

Although it is not clear from the left plot of Fig. 3 that data set is better, the slope $m = 0.686 \pm 0.004$ indicates that the MFI+$K/Ka$+PR4 data set provides a better fit. Moreover, the $K/Ka$+PR4 data set has larger uncertainties which can mask model inconsistencies. On the other hand, from the middle plot of Fig. 3, we observe that the inclusion of the K and Ka WMAP bands to the MFI+PR4 data set improves the goodness of the fit. Finally, comparing MFI+$K/Ka$+PR4 with MFI+$K/Ka$+PR3, we see that PR3 provides a better fit in the Galactic plane, while PR4 fits better outside the Galactic plane (Fig. 2). Since the fit in the Galactic plane is bad in both cases we have chosen the MFI+$K/Ka$+PR4 as our default data set, as it retrieves better fits within the 95 per cent confidence regions (pixels outside the Galactic plane, Fig. 4).

### 5.1.2 Spatial variability

We have also studied the spatial variability of the synchrotron spectral index in several high signal-to-noise regions of the sky (see Fig. 5). These connected regions satisfy the condition that $\beta_s$ is estimated with a signal-to-noise ratio larger than 15. In particular, R1 is associated with the North Polar Spur (NPS), and R2 encompasses the Galactic plane. R3, R4, and R5 are other sky regions where the polarized synchrotron intensity has a large signal-to-noise ratio.

Fig. 6 shows the estimated synchrotron spectral index against the uncertainty on the estimation of all the pixels within a given region. We have limited this study to those pixels with a $\chi^2_{red}$ within the 95 per cent confidence region. The area delimited by the dotted lines contains the values that are consistent within $3\sigma$ with the weighted mean in each region. The top left-hand panel indicates that $\beta_s$ has a large spatial variability across the whole available QUIJOTE-MFI sky (QS). Therefore, a constant value of $\beta_s$ is not a good model of the synchrotron emission. On the contrary, the R1, R3, R4, and R5 pixels values are well within those lines, i.e. a uniform $\beta_s$ value could be a good model for all pixels within each region. Finally, R2 (the top right-hand panel) shows a significant spatial variability which is consistent with the large heterogeneity observed in the $\beta_s$ map.

<footnote n="10">¹⁰The $\chi^2$-distribution with $N_{\text{dof}}$ divided by $N_{\text{dof}}$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Synchrotron spectral index (top), its uncertainty (middle), and reduced $\chi^2$ (bottom) maps obtained after component separation with the default data set MFI+$K$/Ka+PR4. The synchrotron emission is modelled with a power law.</figcaption></figure>

The study of regions with uniform $\beta_s$ values helps with improving the detectability of primordial B modes. Allowing spatial variations of the spectral parameters at the pixel level results in a very robust parametrization of the signal sky. However, this robustness comes at the expense of an increase in the statistical uncertainty of the parameters as less information is provided in the fit (Errard & Stompor 2019). Thus, several approaches have been proposed in the literature to define sky regions with uniform spectral parameters. For example, in Errard & Stompor (2019), these regions are chosen

<figure><img src="image_6.png" /><figcaption>**Figure 5.** R1, R3, R4, and R5 are sky regions where $\beta_s$ is assumed uniform in Section 5.1.2 and R2, which encompasses the Galactic plane seen by QUIJOTE, is a very heterogeneous region. These regions satisfy that $\beta_s$ is recovered with a signal-to-noise ratio larger than 15.</figcaption></figure>

as super-pixels at a lower HEALPix maps resolution, whereas in Grumitt, Jew & Dickinson (2020), the regions are obtained using clustering algorithms such as the mean-shift clustering algorithm. Recently, Puglisi et al. (2022) has presented a new methodology based on spectral clustering to define geometrical affine regions with similar spectral parameters. It is worth noting that if the assumption of uniform spectral parameters within those regions does not hold, the modelling errors introduced might bias cosmological parameters measurements obtained from the output CMB map after component separation, as well as foreground model parameters.

We have calculated the value of $\beta_s$ in some of these regions assuming a constant value within each region. We have performed the fit in the following manner:

(i) First, we fix $\beta_s$ to a given value and fit the rest of the model parameters in each pixel of the region.

(ii) Then, the rest of the parameters are fixed to the estimation from the previous fit, and we fit $\beta_s$ assuming a unique value in the whole region under study.

(iii) $\beta_s$ is fixed to the new obtained value and the process is repeated until it reaches convergence.

We have chosen the median of the $\beta_s$ values (obtained pixel-wise) within that region as the initial guess of $\beta_s$. The results are shown in Table 1. Notice that the uncertainty on the recovered $\beta_s$ has dramatically decreased. This is simply a result of having $N_{pix}^R$ (the number of pixels contained within the region $R$) times more information to fit the parameter. The $\beta_s$ values recovered in each region (R1, R3, R4, and R5) are not consistent among them. These results further showcase the spatial variability of the synchrotron's spectral parameter.

### 5.1.3 Comparison with current $\beta_s$ models

In this section, we compare our $\beta_s$ map with the currently most used $\beta_s$ template,<fnref n="11" /> the 'Model 4' Miville-Deschenes et al. template, which was constructed with Haslam and WMAP observations in temperature (Miville-Deschênes et al. 2008). Fig. 7 shows the distribution of the spectral index value for this model (blue) and for our analysis (orange), considering only those QUIJOTE-MFI pixels

<footnote n="11">¹¹Used for example in the Planck Sky Model (Ashdown et al. 2012), or in the Python Sky Model (PYSM) a PYTHON library to simulate foregrounds (Thorne et al. 2017).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Synchrotron spectral index estimate against its uncertainty within different sky regions: QUIJOTE-MFI sky (QS) (Fig. 1); R1, R2, R3, R4, and R5 are shown in Fig. 5. The solid, dashed, and dotted lines enclose the values of $\beta_s$ within $1\sigma$, $2\sigma$, and $3\sigma$ of the weighted mean, respectively. The study is limited to those pixels whose $\chi^2_{\text{red}}$ lies within the 95 per cent confidence region.</figcaption></figure>

<table><caption>**Table 1.** Synchrotron spectral index estimation $\beta_s^R$ and its uncertainty $\sigma(\beta_s^R)$ obtained, assuming uniform value across the regions R1, R3, R4, and R5 shown in Fig. 5.</caption><thead><tr><th>Region</th><th>$f_{\text{sky}}$ (per cent)</th><th>$\beta_s^R$</th><th>$\sigma(\beta_s^R)$</th></tr></thead><tbody><tr><td>R1</td><td>4.84</td><td>-3.028</td><td>0.002</td></tr><tr><td>R3</td><td>0.96</td><td>-2.945</td><td>0.008</td></tr><tr><td>R4</td><td>0.56</td><td>-3.319</td><td>0.011</td></tr><tr><td>R5</td><td>0.21</td><td>-3.228</td><td>0.019</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Distribution of the synchrotron spectral index from ‘Model 4’ of Miville-Deschênes et al. (2008) and from our estimation using the MFI+K/Ka+PR4 data set. Vertical dashed lines indicate the mean value for each distribution.</figcaption></figure>

that lie within the 95 per cent confidence region of the $\chi^2$ (QS$_{95}$). In the QS$_{95}$ region, the mean and the standard deviation from the ‘Model 4’ of Miville-Deschênes et al. (2008) template are $-3.00 \pm 0.05$ while those from our estimate are $-3.08 \pm 0.13$. It is interesting to note that the variability observed in our analysis is significantly larger. A direct comparison of the dispersion of both maps (using the same mask) indicates an increment of the spatial variability in our study around a factor of 2.6, i.e. $\sigma(\beta_s^{\text{MFI+K/Ka+PR4}})/\sigma(\beta_s^{\text{Model4}}) \sim 2.6$.

One may wonder if this result can be affected by the considered prior, since the estimated spectral indices for low signal-to-noise pixels are significantly constrained by it (see Section 5.5.2). In order to test this point, we have repeated the previous analysis considering only those pixels satisfying that the recovered $\beta_s$ values have a signal-to-noise larger than 15 (i.e. where the synchrotron signal-to-noise ratio is high and thus the results are not driven by the prior) and lie within the 95 per cent confidence region of the $\chi^2$ (HS2N$_{95}$). In this case, we find that the mean value and dispersion of the distribution of $\beta_s$ are $-3.12 \pm 0.15$ for our analysis (see green histogram in Fig. 7) versus $-3.00 \pm 0.05$ for ‘Model 4’ in the same region, confirming our finding. Although our estimations can be affected by the presence of noise, the results show that the variability of the synchrotron spectral index assumed in current templates is underestimated. A similar increment in the variability was also noted by analysing the S-PASS data in the Southern hemisphere (Krachmalnicoff et al. 2018).

Recently, Weiland et al. (2022) published a composite map of $\beta_s$, using publicly available data covering approximately 44 per cent of the sky. In the region covered in our study, they obtained $\beta_s$ estimates in the Galactic plane and the North Polar Spur using information from WMAP K and Ka band, and estimates at latitudes larger than 40° using K, Ka, and DRAO 1.41 GHz map (Wolleben et al. 2006). From a visual inspection, our results are compatible within the North Polar Spur. We find that our derived spectral indices are steeper at the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top row: Synchrotron curvature estimate (left) and uncertainty (right) maps obtained after component separation using the default data set (MFI+*K*/Ka+PR4). The synchrotron emission is modelled using a power law with spatially varying curvature (pixel-wise). Bottom row: Reduced $\chi^2$ map (left) and $c_s$ signal-to-noise map (right).</figcaption></figure>

Galactic plane. Weiland et al. (2022) found discrepancies between the $\beta_s$ values obtained in the Fan Region when they performed the analysis using WMAP K and Ka band versus WMAP K band and Planck LFI 30 GHz channel. In the latter case, the recovered $\beta_s$ were significantly steeper. We repeated our analysis excluding the PR4 30 GHz channel and did not observe a discrepancy concerning the $\beta_s$ recovered from the default analysis in Fan Region. This results from the fact that the $\beta_s$ recovery is mainly driven by QUIJOTE-MFI data. At high latitudes, we cannot make a reasonable comparison since our $\beta_s$ estimates are driven by the prior. They also show that DRAO data have some unexplained systematics and can be affected by Faraday rotation depolarization.

Other studies, such as those presented in Vidal et al. (2015), Fuskeland et al. (2014, 2021), Martire et al. (2022), also find variability of the spectral index analysing different regions of the sky. However, it is difficult to compare the same regions in our map, since they compute a global spectral index for large areas, while we work pixel by pixel. For example, near the centre of the Galactic plane we see a fair amount of structure that cannot be accounted for in the T–T scatter plots analyses carried out in some of the cited papers, that use several pixels to obtain a single $\beta_s$ value. In that sense, the methodology followed here is more complete given that we perform a full component separation in each pixel, retrieving information at smaller scales for a large fraction of the sky.

Rubiño-Martín et al. (2022) obtain an estimate of the synchrotron spectral index map directly from the comparison of the QUIJOTE-MFI 11 GHz map with the WMAP K-band map. The results obtained there are fully consistent with the ones from this work.

## 5.2 Synchrotron curvature

We have also explored a synchrotron model with curvature, i.e. the model presented in equation (3), using the MFI+$K$/Ka+PR4 data set. Fig. 8 shows the estimation and uncertainty maps of the curvature parameter as well as the $\chi^2_{red}$ map and the $c_s$ signal-to-noise map.

We observe from the signal-to-noise map that curvature is detected at more than $3\sigma$ in the Galactic plane, in regions where the fit is not good as it can be seen from the $\chi^2_{red}$ map. Even though the inclusion of a curvature parameter is not able to explain the complexity of this region, this parameter can account for some effects along the Galactic plane, e.g. Faraday rotation.

Outside the Galactic plane the estimated $c_s$ values are close to zero and their uncertainties are around 0.1, which are the expected value and the spread of the prior set on $c_s$. Moreover, the recovered $\beta_s$ map in this case is very similar to the one obtained when the synchrotron is model with a power law. This means that we do not have enough sensitivity to detect a spatially varying curvature. Hopefully, joint analyses with future releases of the Northern celestial hemisphere data like the new MFI2 instrument and C-BASS at 5 GHz (Jones et al. 2018) might elucidate more details on changes of the power-law spectrum.

In Fig. 9, we compare the goodness of fit using a power law versus a power law with curvature as the synchrotron model. We see that there are more points located below the bisector. Besides, the slope 0.9227 ± 0.0005 calculated at the 95 per cent confidence region, shows that, given the current data, the power-law model is slightly preferred over the power law plus curvature model.

Furthermore, we have considered modelling the synchrotron emission with a power law with uniform curvature. We have assumed a constant $c_s$ in four regions: RC1, RC2, and the Haze and North bubble (Fig. 10). The recovered curvature values are shown in Table 2. RC1 encompasses all the pixels whose $\chi^2_{red}$ is within 95 per cent confidence region. RC2 is composed of the RC1 pixels that also satisfy that the synchrotron polarized intensity signal-to-noise ratio at 30 GHz is larger than 5. We detect curvature in all regions. The detection is more evident in RC1 and RC2, mostly due to the higher sensitivity (lower $\sigma_C$) in these regions. However, it is important to highlight that there is no physical reasoning behind the definition of RC1 and RC2, and the assumption of uniform curvature in all synchrotron
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Reduced $\chi^2$ calculated using a power law as a model of the synchrotron emission ($\chi^2_{red,pl}$) versus $\chi^2_{red}$ when the model is a power law with spatially varying curvature ($\chi^2_{red,plc}$). The colour scale is related to the density of points, redder (bluer) corresponds to denser (sparser) regions. The red rectangle shows the $\chi^2_{red}$ within a 95 per cent confidence region. The slope at the 95 per cent confidence region is $m = 0.9227 \pm 0.0005$, shown with a green dashed line. The orange dashed line shows the one-to-one line.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Regions where $c_s$ has been assumed to be uniform.</figcaption></figure>

<table><caption>**Table 2.** Estimated values of the curvature and its uncertainty obtained assuming the curvature is uniform within the region.</caption><thead><tr><th>Region</th><th>$f_{sky}$ (%)</th><th>$c_s^R$</th><th>$\sigma_{c_s}^R$</th><th>$|c_s^R| / \sigma_{c_s}^R$</th></tr></thead><tbody><tr><td>RC1</td><td>45.48</td><td>-0.0797</td><td>0.0012</td><td>63.75</td></tr><tr><td>RC2</td><td>5.93</td><td>-0.2768</td><td>0.0017</td><td>161.57</td></tr><tr><td>Haze</td><td>0.94</td><td>0.041</td><td>0.010</td><td>4.23</td></tr><tr><td>North bubble</td><td>0.63</td><td>-0.083</td><td>0.007</td><td>11.43</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison between the pixel $\beta_s$ values obtained when fitting the synchrotron emission with a spatially varying curvature model (y-axis) versus with a model with uniform curvature (x-axis) in the regions RC1 and RC2 using the MFI+K/Ka+PR4 data set.</figcaption></figure>

high signal-to-noise regions is arbitrary.<fnref n="12" /> In the Haze and North bubble, we find a curvature value different from zero at more than $3\sigma$. These regions are studied in greater detail in Guidi et al. (2022).

We have studied how $\beta_s$ changes when we impose the constraint of having a uniform $c_s$ value within each region. The results are displayed in Fig. 11. For RC2, we observe that $\beta_s$ steepens considerably. The weighted mean value of $\beta_s$ in RC2 is $\langle\beta_s\rangle = -3.022 \pm 0.011$ in the pixel-wise analysis and, $\langle\beta_s\rangle = -3.375 \pm 0.002$ when $c_s$ is imposed to be uniform in RC2. For RC1, this effect is not as considerable. The weighted mean values are $\langle\beta_s\rangle = -3.079 \pm 0.002$ and $\langle\beta_s\rangle = -3.1651 \pm 0.0014$ when $c_s$ varies pixel-wise and is uniform, respectively. The steepening of $\beta_s$ leads to values of the exponent $\beta_s + c_s \log(v/v_s)$ within [-3.04, -3.10] at 11 GHz which are compatible with the average value of $\beta_s$ when we fit to a power-law model. From these results, we infer that the $\beta_s$ and $c_s$ parameters are not independent. More sensitive data at the QUIJOTE frequencies and at lower and/or higher frequencies are required to break the degeneracy.

In order to test which model provides a better goodness of fit, we calculate the reduced $\chi^2$ of a given region R as follows:

$$ \chi^2_{\text{red,R}} = \frac{1}{N_{\text{dof}}} \sum_{p=1}^{N_{\text{pix}}} \sum_{I \in \{Q,U\}} (\mathbf{d}_{p,i} - \mathbf{S}_{p,i}) \mathbf{C}_{p,i}^{-1} (\mathbf{d}_{p,i} - \mathbf{S}_{p,i}), \quad (11) $$

where we sum over all pixels $N_{pix}^R$ within R. The dof are given as $N_{dof} = N_{pix}^R(2N - N_\theta)$ when all model parameters are allowed to vary pixel-wise, and $N_{dof} = N_{pix}^R(2N - (N_\theta - 1)) - 1$ when $c_s$ is assumed uniform in the analysis, where $N_\theta$ is the number of model parameters. We calculated the value of this estimator in three cases: (i) when the model parameters are allowed to vary spatially using a power-law model for the synchrotron component, (ii) when the model parameters vary from pixel-to-pixel using a power law with curvature model, (iii) when we fit the data assuming uniform curvature using a power law with curvature model. The results are given in Table 3.

<footnote n="12">¹²Any curvature will be more easily detected in high signal-to-noise regions than in low signal-to-noise regions.</footnote>
<table><caption>**Table 3.** Reduced $\chi^2$ obtained using either a power law or a power law with curvature model in different regions, R. We have considered two curvature models: one where $c_s$ varies spatially (spatial) and other where $c_s$ is assumed constant in R.</caption><thead><tr><th>Model</th><th>Curvature</th><th>Region</th><th>χ<sup>2</sup><sub>red,R</sub></th></tr></thead><tbody><tr><td>Power law</td><td>–</td><td>RC1</td><td>0.892</td></tr><tr><td>Power law + curvature</td><td>Spatial</td><td>RC1</td><td>0.965</td></tr><tr><td>Power law + curvature</td><td>Uniform</td><td>RC1</td><td>0.936</td></tr><tr><td>Power law</td><td>–</td><td>RC2</td><td>1.010</td></tr><tr><td>Power law + curvature</td><td>Spatial</td><td>RC2</td><td>1.088</td></tr><tr><td>Power law + curvature</td><td>Uniform</td><td>RC2</td><td>1.081</td></tr><tr><td>Power law</td><td>–</td><td>Haze</td><td>0.845</td></tr><tr><td>Power law + curvature</td><td>Spatial</td><td>Haze</td><td>0.936</td></tr><tr><td>Power law + curvature</td><td>Uniform</td><td>Haze</td><td>0.885</td></tr><tr><td>Power law</td><td>–</td><td>North bubble</td><td>0.961</td></tr><tr><td>Power law + curvature</td><td>Spatial</td><td>North bubble</td><td>1.041</td></tr><tr><td>Power law + curvature</td><td>Uniform</td><td>North bubble</td><td>0.986</td></tr></tbody></table>

The $\chi^2_{red}$ results show that the models we used, i.e. power law and power law with curvature, are compatible with the data. However, there is not enough statistical significance to discern which model suits better the data. Especially, considering that we have not been able to take into account possible correlations between pixels and that the power law with curvature model is degenerate.<fnref n="13" />

## 5.3 Recovered amplitudes and comparison with Planck results

We have compared our baseline results, i.e. using the MFI+$K/Ka$+PR4 data set and a power law as the synchrotron model, to those obtained from the Commander pipeline (Eriksen et al. 2008) applied to PR4 data.<fnref n="14" /> We have only considered this pipeline among those used by Planck, since it is the reference method with regard to the recovery of foreground components. In Figs 12–14, we show a comparison of the CMB, the synchrotron emission at 30 GHz, and the thermal dust emission at 353 GHz between Commander and our results. In order to perform a direct comparison we have filtered Commander results with FDEC. The left column shows the Q and U Commander amplitudes, the centre column our amplitudes and the right column the corresponding uncertainties. A visual inspection shows that both estimates are very similar, especially the synchrotron and thermal dust emissions that are the dominant contributions in polarization.

### 5.3.1 CMB

Regarding CMB, the left column of Fig. 15 shows the pixel-to-pixel comparison for the recovered CMB map from our analysis and from Commander both in Q and U. We have applied a combination of the QUIJOTE observed sky and the common polarization confidence mask provided by the Planck Collaboration<fnref n="15" /> (Planck Collaboration IV 2020d).

We observe from the maps that there is a discrepancy. We found that the application of the FDEC filter, before the component

separation process, leads to a decrease of the amplitude in the power spectra of our recovered CMB map. This power reduction appears only when Planck and WMAP are filtered with FDEC, since the CMB information is extracted mainly from those channels. Instead of applying the FDEC filter, one could apply a filter that suppresses the large scales. This would be equivalent to applying a linear function to the CMB and there would not be a reduction of power. However, since we want to study all scales we decided to apply the FDEC filter. Since the aim of this work is the study of the foregrounds, we keep the results obtained with all the data filtered with FDEC to recover the $\beta_s$ map without any bias. One can in principle recover the unbiased CMB following one of the approaches described below:

(i) Perform the component separation analysis without filtering the data with FDEC and including the FDEC correction in QUIJOTE-MFI data as part of the model; or

(ii) Given the unbiased $\beta_s$ map<fnref n="16" /> and Planck data, one can construct a template with the modes that QUIJOTE-MFI data are missing after being filtered with FDEC. Then perform the analysis with the reconstructed QUIJOTE-MFI maps.

Since the estimation of the CMB is out of the scope of this paper, we leave this analysis for future works.

### 5.3.2 Synchrotron

Fig. 16 shows the difference between the synchrotron amplitude maps obtained using the MFI+$K/Ka$+PR4 and the Commander reconstruction using the PR4 data. The largest differences observed are located in the Galactic plane where the model fails to reproduce the sky signal. We also observe large-scale structures in the difference map. These structures can originate from the fact that we have obtained a more accurate estimation of the scaling law, as our fit is performed using additional frequencies. However, overall, the correlation between both methods is very good.

This can also be seen in the centre column of Fig. 15, where a pixel-to-pixel comparison is given, showing that both methods present a synchrotron amplitude at 30 GHz highly correlated for Q and U except in some pixels where the synchrotron emission is very large. Those pixels are located primarily in the Galactic plane. These discrepancies are likely to arise from differences in the amplitude of the polarized intensity instead of from differences in the polarization angles. In Fig. 15, we observe that both the slopes, in the Q and U plots, are higher than unity. If the discrepancies were originated from differences in the polarization angle, one slope would be higher than unity and the other lower.

### 5.3.3 Dust

Regarding thermal dust emission, this foreground strongly dominates the 353 GHz Planck frequency map and, therefore, the recovered amplitude is very much determined by this channel. This was also the case in the Commander analysis done by the Planck Collaboration. Thus, our recovered Q and U components of the thermal dust are strongly correlated with those obtained using Commander, see the right column of Fig. 15.

<footnote n="13">¹³We considered applying other statistics such as the Bayesian evidence to do model selection. However, since the QUIJOTE-MFI noise simulations are not end-to-end and the Bayesian evidence is very computationally expensive we did not perform any model selection analysis. This is left for future work.</footnote>

<footnote n="14">¹⁴Data available at NERSC under /cmb/daa/planck20.</footnote>

<footnote n="15">¹⁵Available at https://pla.esac.esa.int/#maps.</footnote>

<footnote n="16">¹⁶Obtained in the component separation analysis using the data filtered with FDEC.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left column: Commander *Q* (top) and *U* (bottom) CMB maps at *N*<sub>side</sub> = 64, smoothed with a Gaussian beam to a final resolution of FWHM = 2°. Centre column: CMB *Q* and *U* maps using the MFI+K/Ka+PR4 data set. Right column: Uncertainty of the CMB maps. Maps are in thermodynamic temperature (μK). We apply the common polarization confidence mask provided by *Planck*.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Left column: Commander *Q* (top) and *U* (bottom) synchrotron amplitude maps at 30 GHz at *N*<sub>side</sub> = 64, smoothed with a Gaussian beam to a final resolution of FWHM = 2°. Centre column: Our estimate of the synchrotron amplitude at 30 GHz, using the MFI+K/Ka+PR4 data set. Right column: Uncertainty of the estimated synchrotron amplitude. Maps are in antenna temperature (μK).</figcaption></figure>

## 5.4 Dust spectral parameters

Although the frequencies of QUIJOTE-MFI do not overlap with the spectral range where the thermal dust is more dominant, we have studied whether the inclusion of this data set in the analysis can help with the thermal dust characterization due to an improvement on the determination of the rest of the polarized foreground parameters. Fig. 17 shows the thermal dust spectral index $\beta_d$ recovered with the

default data set, modelling the synchrotron emission as a power law, in two cases:

(i) $T_d$ is included as a model parameter.

(ii) $T_d$ is fixed to Commander's estimation of the thermal dust temperature from the component separation analysis in intensity (Planck Collaboration 2016a) like Commander did in their polarization analysis. Fixing $T_d$ helps breaking its degeneracy with $\beta_d$ in
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Left column: Commander *Q* (top) and *U* (bottom) thermal dust amplitude maps at 353 GHz at *N*<sub>side</sub> = 64, smoothed with a Gaussian beam to a final resolution of FWHM = 2°. Centre column: Our estimate of the thermal dust amplitude at 353 GHz, using the MFI+*K*/K*a*+PR4 data set. Right column: Uncertainty of the estimated thermal dust amplitude. Maps are in antenna temperature (μK).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Comparison of CMB (left), synchrotron at 30 GHz (centre), and thermal dust at 353 GHz (right) amplitudes recovered using the MFI+*K*/K*a*+PR4 data set and the ones obtained by Commander using PR4 data. The correlation factors are ρ<sup>Q</sup> = 0.543 and ρ<sup>U</sup> = 0.817 (CMB), ρ<sup>Q</sup> = 0.992 and ρ<sup>U</sup> = 0.973 (synchrotron), and ρ<sup>Q</sup> = 1.000 and ρ<sup>U</sup> = 0.997 (thermal dust).</figcaption></figure>

the Rayleigh–Jeans part of the thermal dust spectrum, which is the
one observed with *Planck* in polarization.

In both maps, we find that the recovered β<sub>d</sub> values are close to
the expected value of the prior, i.e. 1.55, except close to the Galactic

plane where the thermal dust signal is larger.<fnref n="17" /> The results differ

<footnote n="17">¹⁷Notice that the uncertainty does not improve in the regions where the β<sub>d</sub> values are close to the mean value of the prior when we fix one parameter. The uncertainty in those pixels is the spread of the prior.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Difference between the synchrotron amplitude $a_s^Q$ ($a_s^U$) obtained with the MFI+K/Ka+PR4 and the Commander estimate, top row (bottom row). Maps are in antenna temperature (µK).</figcaption></figure>

significantly along the Galactic plane (see Fig. 18). This difference originates since our recovered $T_d$ map does not resemble the used $T_d$ template, as shown in Fig. 19. We remark that although in the first case $T_d$ is estimated from the polarization analysis, the $T_d$ recovered values lie close to the expected value of the prior (22 K) except along the Galactic plane where the fit is not good. Moreover, it is

<figure><img src="image_19.png" /><figcaption>**Figure 18.** $\beta_d$ (top row) and $\beta_s$ (bottom row) relative difference map between the maps obtained when we include $T_d$ as a model parameter and when we fix it.</figcaption></figure>

very difficult to fit $T_d$ from polarization data only, as the highest frequency is 353 GHz, and thus we are not able to trace the thermal dust peak.

In Fig. 18, we show the relative difference between spectral index map of the thermal dust and synchrotron obtained when $T_d$ is included

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Left column: Estimate (top) and uncertainty (bottom) of thermal dust spectral index obtained when $T_d$ is included as a model parameter. Right column: Estimate (top) and uncertainty (bottom) of thermal dust spectral index obtained when the $T_d$ template obtained by Commander in the intensity analysis is used to fix $T_d$ in the component separation process.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 19. Top row: Thermal dust temperature map recovered in the default case. Bottom row: Difference map between the top row map and the $T_d$ template used in the analysis. Maps are in Kelvin.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 20. $\beta_d$ relative difference map between the map obtained using the MFI+$K/Ka$+PR4 and the one obtained with $K/Ka$+PR4 data sets when we include $T_d$ as a model parameter (top row) and when we fix it (bottom row).</figcaption></figure>

as a model parameter and when it is fixed. The relative difference is calculated as follows:

$$ \widehat{\Delta\beta}_{p,1,2} = \frac{\beta_{p,1} - \beta_{p,2}}{\sqrt{\sigma_{p,\beta_1}^2 + \sigma_{p,\beta_2}^2 - 2\sigma_{p,\beta_1,\beta_2}}}, \quad (12) $$

where $\sigma_{\beta_1}^2$ ($\sigma_{\beta_2}^2$) is the variance of the $\beta_1$ ($\beta_2$) map, and $\sigma_{\beta_1,\beta_2}$ is the covariance between the $\beta_1$ and $\beta_2$ maps that are being compared. As expected from Fig. 17 the differences close to the Galactic plane are significantly large in the case of $\beta_d$. On the other hand, we find that, the $\beta_s$ maps recovered in both cases are compatible and the differences resemble Gaussian noise except along the Galactic plane where the model fails.

We also studied the relative difference between the $\beta_d$ map obtained with the MFI+$K/Ka$+PR4 and $K/Ka$+PR4 data sets in Fig. 20. The top panel shows the relative difference when $T_d$ is included as a model parameter and the bottom panel when $T_d$ is fixed. We observed that both maps are compatible except in regions where the fit is not good. Moreover, when we compare the uncertainty maps we find that there is not a significant improvement when we include QUIJOTE-MFI channels. Thus, we conclude that the improvement in the characterization of low-frequency foregrounds does not help necessarily with the estimation of thermal dust spectral parameters.

## 5.5 Goodness of fit

In this section, we study in depth the quality of the results obtained using the default data set. In Section 5.5.1, we analyse the $\chi^2$ distribution of the results as well as the Q and U residuals of each channel. Section 5.5.2 investigates the robustness of our results regarding the estimation of the synchrotron spectral index with respect to the prior applied to this parameter.

<figure><img src="image_23.png" /><figcaption>Figure 21. $\chi_p^2$ distribution obtained using the default data set. The orange curve shows the theoretical $\chi^2$ probability density function with $N_{dof} = 13$. The area to the left of the grey dashed line shows values within the 95 per cent confidence region.</figcaption></figure>

### 5.5.1 $\chi^2$ distribution and residuals

We have studied the pixel $\chi^2$ distribution obtained from the fit using MFI+$K/Ka$+PR4 (see Fig. 21):

$$ \chi_p^2 = N_{dof} \cdot \chi_{red}^2. \quad (13) $$

Moreover, we have also calculated the residuals per channel involved in the analysis:

$$ r_{p,v} = \frac{(d_{p,v} - S_{p,v})}{\sigma_{p,v}}. \quad (14) $$
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Q and U residual maps for each frequency channel at $N_{\text{side}} = 64$. Maps are displayed in thermodynamic temperature (µK).</figcaption></figure>

In the perfect scenario, residuals maps are consistent with instrumental noise alone. Therefore, they are a valuable tool to look for either systematic effects or mismatches in the foreground modelling.

First of all, we recall that the number of dof for this analysis is 13 [11 channels × 2 (Q and U) minus 9 free parameters]. We find $\langle\chi_p^2\rangle = 14.3$ and $\sigma = 8.9$ slightly larger than what is expected for the theoretical number of dof. Fig. 21 shows that the $\chi_p^2$ values follow a $\chi^2$-like distribution, whose peak lies close to $N_{\text{dof}} = 13$. However, there is an excess of pixels at large values of $\chi^2$ with respect to the $\chi^2_{N_{\text{dof}}}$ distribution. That excess appears since there are pixels where the model is not able to track the true sky emission, mainly in the Galactic plane. Thus, those pixels are highly inconsistent with this $\chi^2$ distribution.

Fig. 22 shows the Q and U residuals maps of every frequency channel from the MFI+K/Ka+PR4 data set. We find that Planck and WMAP residuals maps are reasonably consistent with the expected noise except along the Galactic plane. The residuals in this region are a consequence of an incorrect modelling of the sky as we saw in the $\chi^2_{\text{red}}$ maps. For the MFI channels we observe that the largest residuals are located in compact regions along the Galactic plane. We observe in the 11 GHz U channel a redder region in the NPS's closest part to the Galactic centre. This region overlaps with the area where we obtain a better goodness of fit if Faraday rotation effects are taken into account (see Fig. B2 in Appendix B). Furthermore, artefacts that resemble the FDEC morphology are present in MFI 13 GHz.

In light of these tests, we are confident of the results obtained in those pixels that are properly modelled by our assumed parametric model. The pixels outside the confidence region are located mainly in the Galactic plane, probably because our model fails to account for the complexity of this region. It would be convenient to study these regions in more detail with more complex models. However, the aim of this work is to study the diffuse components and the study of specific regions has been conducted in other works (Watson et al. in preparation; Ruiz-Granados et al. in preparation; Guidi et al. 2022).

### 5.5.2 Robustness with respect to the prior

As previously stated, the use of prior information is essential in Bayesian analysis, helping with convergence and computational time reduction. Besides, when the data do not have enough sensitivity, i.e. there is not enough information to obtain a reliable estimation of the spectral index, the prior tends to provide a value close to the mean value of the distribution. In other words, a conservative value is assigned to the spectral index in those pixels. Thus, in order to detect which pixels are prior dependent, we have also performed component separation using two additional Gaussian priors on $\beta_s$, $\mathcal{N}(-3.1, 0.6)$ and $\mathcal{N}(-3.0, 0.3)$. The $\beta_s$ estimation and uncertainty maps with these new priors are shown in Fig. 23 together with those obtained with the prior used in the default analysis (left columns).

Comparing the results using the default prior, i.e. $\mathcal{N}(-3.1, 0.3)$, versus a less restrictive prior, i.e. $\mathcal{N}(-3.1, 0.6)$, we observe that the uncertainty on the recovered $\beta_s$ increases at the prior-dominated pixels. On the other hand, in those regions where the synchrotron emission is very intense the uncertainty remains the same. Likewise, the estimated $\beta_s$ in the latter pixels are very similar whereas the other pixels are visually different. The $\beta_s$ distribution of the pixels outside the low-uncertainty regions are compatible with the prior distribution. This is the reason why the estimated values are different and the spread is larger when the prior is relaxed.

When we use a prior with a different expected value, i.e. $\mathcal{N}(-3.0, 0.3)$, but equal standard deviation we obtain a similar uncertainty map. The estimated $\beta_s$ is almost the same in the low-uncertainty regions, i.e. the high-intensity synchrotron regions. However, a flatter spectrum (closer to -3.0 instead of -3.1) is recovered outside those areas. This is more evident from the bottom panel of Fig. 24 where the difference between the $\beta_s$ map estimated with the default prior and the $\mathcal{N}(-3.0, 0.3)$ prior is shown. Outside the regions where the synchrotron emission is the largest, the difference is close to -0.1 which is the difference between the expected value of the priors. In other words, when there is not enough information from the data the recovered $\beta_s$ is close to the expected value of the prior. This is an advantage of using prior information, since it assigns a
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Synchrotron spectral index estimate (top row) and uncertainty (bottom row) obtained using different Gaussian prior distributions and the default data set (MFI+*K*/Ka+PR4). The synchrotron emission is modelled as a power law.</figcaption></figure>

conservative value to the spectral index instead of unphysical values or simply failing to perform the fit.

## 6 CONCLUSIONS

In this work, we have presented the component separation products in polarization obtained from combining the QUIJOTE-MFI data at 11 and 13 GHz, with the WMAP *K* and *Ka* bands and all *Planck* polarized channels. We have seen that the inclusion of the QUIJOTE-MFI data is crucial to improve the parameter estimation of the low-frequency foregrounds, in particular for the estimation of the synchrotron spectral index.

We have obtained the first detailed $\beta_s$ map of the Northern Celestial hemisphere at a scale of 2°, assuming the synchrotron emission is modelled as a power law. This model represents well the data except in the Galactic plane where the physics might be more complex. We find, using the pixels whose $\chi^2_{red}$ lies within the 95 per cent confidence region, an average value of -3.08 and a dispersion of 0.13. The latter is broader than the dispersion of commonly used $\beta_s$ templates. Moreover, we have found that the spectral index is not compatible with a uniform value, i.e. there are statistical significant differences of $\beta_s$ across the observable sky.

We have also modelled the synchrotron emission as a power law with curvature. The pixel-based analysis of the curvature shows that $c_s$ is only detected in some regions in the Galactic plane where the fit is bad. When we assume a model with uniform curvature in RC1 (the region that includes all pixels whose $\chi^2_{red}$ is within the 95 per cent confidence region for the power law with curvature model) we found a $c_s = -0.0797 \pm 0.0012$. We found that both models, i.e. power law and power law with uniform curvature, provide a good fit given the available data. However, there is not enough statistical significance to distinguish which model is better. A more thorough study is left for future work.

We found that our recovered synchrotron and thermal dust maps are highly correlated with the maps presented by the *Planck* collaboration using *Commander*, even though we found some large-scale difference between the synchrotron emission maps that arise from better estimation of the SED due to the addition of more frequency

<figure><img src="image_26.png" /><figcaption>**Figure 24.** Difference map between the estimated $\beta_s$ using the default prior, i.e. $\mathcal{N}(-3.1, 0.3)$, and the one obtained using an alternative prior, see Fig. 23. Top: $\mathcal{N}(-3.1, 0.6)$. Bottom: $\mathcal{N}(-3.0, 0.3)$.</figcaption></figure>

channels. On the other hand, we recovered a CMB with less power when we use the filtered *K*, *Ka*, and PR4 with FDEC. Since our analysis focuses on the characterization of foregrounds, we keep the results obtained with the filtered maps. However, as commented in Section 5.3.1 an unbiased CMB map can be recovered following other approaches.
We have also performed different analyses to test the validity of our results. First, we found that our results are compatible with a $\chi^2$ distribution in those pixels where the power-law model fits well the data. Furthermore, we have calculated the normalized residuals of the pixels with an acceptable goodness of fit of all frequency channels and they are all consistent within the $3\sigma$ level. Finally, we have evaluated the robustness of the estimated $\beta_s$ varying the prior imposed in this parameter. We found that the estimations in the high signal-to-noise synchrotron areas are prior independent, while outside these regions the prior governs the $\beta_s$ estimation.

## ACKNOWLEDGEMENTS

We thank the staff of the Teide Observatory for invaluable assistance in the commissioning and operation of QUIJOTE. The QUIJOTE experiment is being developed by the Instituto de Astrofísica de Canarias (IAC), the Instituto de Física de Cantabria (IFCA), and the Universities of Cantabria, Manchester, and Cambridge. Partial financial support was provided by the Spanish Ministry of Science and Innovation under the projects AYA2007-68058-C03-01, AYA2007-68058-C03-02, AYA2010-21766-C03-01, AYA2010-21766-C03-02, AYA2014-60438-P, ESP2015-70646-C2-1-R, AYA2017-84185-P, ESP2017-83921-C2-1-R, AYA2017-90675-REDC (co-funded with EU FEDER funds), PGC2018-101814-B-I00, PID2019-110610RB-C21, PID2020-120514GB-I00, IACA13-3E-2336, IACA15-BE-3707, EQC2018-004918-P, the Severo Ochoa Programs SEV-2015-0548 and CEX2019-000920-S, the María de Maeztu Program MDM-2017-0765, and by the Consolidador-Ingenio project CSD2010-00064 (EPI: Exploring the Physics of Inflation). We acknowledge support from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias, and the European Regional Development Fund (ERDF) under grant with reference ProID2020010108. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement number 687312 (RADIOFOREGROUNDS). EdH acknowledges financial support from the Concepción Arenal Programme of the Universidad de Cantabria. DT acknowledges the support from the Chinese Academy of Sciences (CAS) President's International Fellowship Initiative (PIFI) with grant no. 2020PM0042. FP acknowledges support from the Spanish State Research Agency (AEI) under grant number PID2019-105552RB-C43. The authors acknowledge the computer resources, technical expertise, and assistance provided by the Spanish Supercomputing Network (RES) node at Universidad de Cantabria. Some of the presented results are based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA) and the Planck Legacy Archive (PLA). Support for LAMBDA is provided by the NASA Office of Space Science. Some of the results in this paper have been derived using the HEALPIX package (Górski et al. 2005), and the HEALPY (Zonca et al. 2019), NUMPY (Harris et al. 2020), EMCEE (Foreman-Mackey et al. 2013), and MATPLOTLIB (Hunter 2007) PYTHON packages.

## DATA AVAILABILITY

The parameter maps obtained from the component separation analysis in the default case, i.e. with the MFI+$K/K_a$+PR4 data set using a power law to model the synchrotron emission, are included in

the released data products associated with the QUIJOTE-MFI wide survey.

These data products as well as the maps can be freely downloaded from the QUIJOTE web page,<fnref n="18" /> as well as from the RADIO-FOREGROUNDS platform.<fnref n="19" /> They include also an Explanatory Supplement describing the data formats. Any other derived data products described in this paper are available upon request to the QUIJOTE collaboration.

## REFERENCES

Abazajian K. N. et al., 2016, preprint (arXiv:1610.02743)

Ade P. et al., 2019, *J. Cosmology Astropart. Phys.*, 2019, 056

Ali-Haïmoud Y., 2013, *Adv. Astron.*, 2013, 462697

Ashdown M. et al., 2012, Astrophysics Source Code Library, record ascl:1208.005

Bennett C. L. et al., 2013, ApJS, 208, 20

BICEP/Keck Collaboration, 2021, *Phys. Rev. D*, 103, 022004

Carretti E. et al., 2019, MNRAS, 489, 2330

Carron J., Lewis A., Challinor A., 2017, *J. Cosmol. Astropart. Phys.*, 2017, 035

de la Hoz E., Vielva P., Barreiro R. B., Martínez-González E., 2020, *J. Cosmol. Astropart. Phys.*, 2020, 006

de la Hoz E., Diego-Palazuelos P., Martínez-González E., Vielva P., Barreiro R. B., Bilbao-Ahedo J. D., 2022, *J. Cosmol. Astropart. Phys.*, 2022, 032

Dickinson C. et al., 2018, *New Astron. Rev.*, 80, 1

Diego-Palazuelos P., Vielva P., Martínez-González E., Barreiro R. B., 2020, *J. Cosmol. Astropart. Phys.*, 2020, 058

Draine B. T., Lazarian A., 1999, ApJ, 512, 740

Eriksen H. K., Jewell J. B., Dickinson C., Banday A. J., Górski K. M., Lawrence C. R., 2008, ApJ, 676, 10

Errard J., Stompor R., 2019, *Phys. Rev. D*, 99, 043529

Fixsen D. J., 2009, ApJ, 707, 916

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Fuskeland U., Wehus I. K., Eriksen H. K., Næss S. K., 2014, ApJ, 790, 104

Fuskeland U. et al., 2021, A&A, 646, A69

Génova-Santos R. et al., 2017, MNRAS, 464, 4107

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

Greaves J. S., Scaife A. M. M., Frayer D. T., Green D. A., Mason B. S., Smith A. M. S., 2018, Nature Astron., 2, 662

Grumitt R. D. P., Jew L. R. P., Dickinson C., 2020, MNRAS, 496, 4383

Guidi F. et al., 2022, MNRAS, 519, 3460

Hanany S. et al., 2019, BAAS, 51, 194

Harris C. R. et al., 2020, Nature, 585, 357

Herranz D. et al., 2022, MNRAS, 519, 3526

Hinshaw G. et al., 2003, ApJS, 148, 63

Hui H. et al., 2018, in *Zmuidzinas J., Gao J.-R., eds, SPIE Conf. Ser. Vol. 10708, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX*, SPIE, Bellingham, p. 1070807

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Hutschenreuter S. et al., 2022, A&A, 657, A43

Jones M. E. et al., 2018, MNRAS, 480, 3224

Krachmalnicoff N. et al., 2018, A&A, 618, A166

Lamagna L. et al., 2020, J. Low Temp. Phys., 200, 374

Lee K. et al., 2020, J. Low Temp. Phys., 200, 384

Lewis A., Challinor A., 2011, Astrophysics Source Code Library, record ascl:1102.026

LiteBIRD Collaboration, 2022, preprint (arXiv:2202.02773)

Martire F. A., Barreiro R. B., Martínez-González E., 2022, *J. Cosmol. Astropart. Phys.*, 2022, 003

<footnote n="18">¹⁸http://research.iac.es/proyecto/cmb/qui-jote.</footnote>

<footnote n="19">¹⁹http://www.radioworlds.eu/.</footnote>
McBride L., Bull P., Hensley B. S., 2022, preprint (arXiv:2207.14213)

Millea M., Anderes E., Wandelt B. D., 2019, Phys. Rev. D, 100, 023509

Miville-Deschênes M. A., Ysard N., Lavabre A., Ponthieu N., Macías-Pérez
J. F., Aumont J., Bernard J. P., 2008, A&A, 490, 1093

Namikawa T. et al., 2022, Phys. Rev. D, 105, 023511

Planck Collaboration X, 2016a, A&A, 594, A10

Planck Collaboration XV, 2016b, A&A, 594, A15

Planck Collaboration I, 2020a, A&A, 641, A1

Planck Collaboration II, 2020b, A&A, 641, A2

Planck Collaboration III, 2020c, A&A, 641, A3

Planck Collaboration IV, 2020d, A&A, 641, A4

Planck Collaboration VI, 2020e, A&A, 641, A6

Planck Collaboration LVII, 2020f, A&A, 643, A42

Puglisi G., Mihaylov G., Panopoulou G. V., Poletti D., Errard J., Puglisi P.
A., Vianello G., 2022, MNRAS, 511, 2052

Ritacco A., Boulanger F., Guillet V., Delouis J.-M., Puget J.-L., Aumont J.,
Vacher L., 2022, preprint (arXiv:2206.07671)

Rubiño-Martín J. A. et al., 2010, in Highlights of Spanish Astrophysics V.
Springer, Berlin/Heidelberg, p. 127

Rubiño-Martín J. A. et al., 2022, MNRAS, 519, 3383

Rybicki G. B., Lightman A. P., 1985, Radiative Processes in Astrophysics.
John Wiley & Sons, New York, NY

Sunyaev R. A., Zeldovich Y. B., 1972, Comm. Astrophys. Space Phys., 4,
173

Thorne B., Dunkley J., Alonso D., Næss S., 2017, MNRAS, 469, 2821

Tristram M. et al., 2022, Phys. Rev. D, 105, 083524

Vidal M., Dickinson C., Davies R. D., Leahy J. P., 2015, MNRAS, 452,
656

Weiland J. L., Addison G. E., Bennett C. L., Halpern M., Hinshaw G., 2022,
ApJ, 936, 24

Wolleben M., Landecker T. L., Reich W., Wielebinski R., 2006, A&A, 448,
411

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K.,
2019, J. Open Source Softw., 4, 1298

APPENDIX A: INDEPENDENT Q AND U
SYNCHROTRON SPECTRAL INDEX

In order to test the assumption of having the same $\beta_s$ in both $Q$
and $U$, we fit $Q$ and $U$ signals independently. Fig. A1 shows the
spectral index, the uncertainty of the spectral index as well as the
reduced $\chi^2$ maps obtained from the three independent fits using the
MFI+$K$/Ka+PR4 data set. We infer from the $\chi^2_{red}$ maps that the fit
outside the Galactic plane is better when $Q$ and $U$ are fitted together.
When we fit just $U$ we observe that the goodness of fit improves
significantly in the Galactic plane. However, this effect is due to the
low signal-to-noise ratio in that area, not due to a better modelling
of the signal.

The $\beta_s^Q$ and $\beta_s^U$ maps are distinctly different. The $\beta_s^{QU}$ map resembles more the $\beta_s^Q$ map. This is expected, since $Q$ has more signal than $U$ in Galactic coordinates. That is also the reason why the uncertainty on the recovered $\beta_s$ is smaller when we fit just $Q$ compared to $U$. However, in those regions where $\sigma_{\beta_s^U}$ is smaller than $\sigma_{\beta_s^Q}$, i.e. regions where $U$ has more signal than $Q$, the $\beta_s^{QU}$ values obtained are closer to those of $\beta_s^U$. This is clearly seen in Fig. A2 where the relative difference between $\beta_s^{QU}$ with respect to $\beta_s^Q$ (top row) and $\beta_s^U$ (bottom row) is shown. The largest differences shown in the top (bottom) panel are located in regions where the signal-to-noise ratio is larger in $U$ ($Q$). On the other hand, the relative difference decreases significantly in the regions where the uncertainty on $\beta_s^Q$ (top) or $\beta_s^U$ (bottom) is smaller.

<figure><img src="image_27.png" /></figure>

**Figure A1.** Synchrotron spectral index estimate (top row) and uncertainty maps (second row) obtained after component separation using the MFI+K/Ka+PR4 data set. The left column shows the $\beta_s$ recovered when we assume that $Q$ and $U$ share the same spectral index, while the centre and right columns depict the $Q$ and $U$ $\beta_s$ when they are assumed to be independent. Bottom row: Reduced $\chi^2$ map for each case study considered. The synchrotron emission is modelled as a power law.
<figure><img src="image_28.png" /><figcaption>**Figure A2.** Relative difference map between the $\beta_s$ map obtained when we assume the same $\beta_s$ in both $Q$ and $U$, and $\beta_s$ recovered from the fit with just $Q$ (top) and just $U$ (bottom). The synchrotron emission is modelled as a power law.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure B1.** Relative difference map between the $\beta_s$ map from the independent $Q$ and $U$ fit using the MFI+$K/Ka$+PR4 data set (top), and using the MFI(FR)+$K/Ka$+PR4 data set (bottom). The synchrotron emission is modelled as a power law.</figcaption></figure>

## APPENDIX B: FARADAY ROTATION

We have also studied the significance of the difference between the $\beta_s^Q$ and $\beta_s^U$ maps, see the top row of Fig. B1. The discrepancies larger than $3\sigma$ are concentrated in the Galactic plane, close to the Galactic centre. This could be a tracer of Faraday rotation. If Faraday rotation is non-negligible at QUIJOTE frequencies, there will be a difference between the polarization angles at QUIJOTE frequencies and those at WMAP/Planck frequencies. This yields a $\beta_s^Q$ map different from $\beta_s^U$ due to the bias introduced by the change in angle. That bias is reasonably cancelled out when combining both $Q$ and $U$ to obtain a single index.

We have studied the possibility of correcting the Faraday rotation effect in the QUIJOTE MFI maps using the model from Hutschenreuter et al. (2022). The rotation of the polarization plane experienced due to the Faraday Rotation effect can be described by

$$ \Delta\phi = RM\lambda^2, \qquad (B1) $$

where $\lambda$ is the wavelength and RM is the rotation measure. We use the RM map estimated by Hutschenreuter et al. (2022) to calculate the rotation angle maps at 11 and 13 GHz QUIJOTE frequencies. Then, QUIJOTE $Q$ and $U$ maps at a given frequency $v$ are de-rotated as follows:

$$ \left( \begin{matrix} Q_{FR} \\ U_{FR} \end{matrix} \right)_v = \left( \begin{matrix} \cos(2\Delta\phi_v) - \sin(2\Delta\phi_v) \\ \sin(2\Delta\phi_v) \cos(2\Delta\phi_v) \end{matrix} \right)_v \left( \begin{matrix} Q \\ U \end{matrix} \right)_v, \qquad (B2) $$

The variance of the de-rotated $Q_{FR}$ and $U_{FR}$ is

$$ \sigma_{Q_{FR}}^2 = \cos^2(2\Delta\phi_v)\sigma_Q^2 + \sin^2(2\Delta\phi_v)\sigma_U^2 + 4[\sin(2\Delta\phi_v)Q + \cos(2\Delta\phi_v)U]^2 \sigma_{\Delta\phi}^2, \qquad (B3) $$

<figure><img src="image_30.png" /><figcaption>**Figure B2.** Difference map between the $\chi_{red}^2$ obtained with the MFI+$K/Ka$+PR4 data set with respect to the $\chi_{red,FR}^2$ obtained with MFI(FR)+$K/Ka$+PR4 data set. In both fits, we have assumed that $Q$ and $U$ share the same spectral indices. The synchrotron emission is modelled as a power law.</figcaption></figure>

$$ \sigma_{U_{FR}}^2 = \sin^2(2\Delta\phi_v)\sigma_Q^2 + \cos^2(2\Delta\phi_v)\sigma_U^2 + 4[\cos(2\Delta\phi_v)Q - \sin(2\Delta\phi_v)U]^2 \sigma_{\Delta\phi}^2. \qquad (B4) $$

Therefore, we have repeated the same analysis but using the MFI(FR)+$K/Ka$+PR4 data set, where MFI(FR) indicates that the QUIJOTE 11 and 13 GHz maps have been de-rotated using the angle obtained from the Hutschenreuter et al. (2022) model, to correct any possible mismatch due to the Faraday Rotation effect (see the bottom row of Fig. B1).
<figure><img src="image_31.png" /><figcaption>**Figure C1.** Relative difference map between the $\beta_s$ template used in the simulation and the $\beta_s$ map from the fit using the simulated data with an FDEC filter applied to all maps (top), and an FDEC filter applied only to QUIJOTE-MFI frequencies (bottom).</figcaption></figure>

We compare these maps (Fig. B1) with the different map between the reduced $\chi^2$ map ($\chi^2_{\text{red}}$) obtained with the MFI+$K$/Ka+PR4 data set with respect to the reduced $\chi^2$ ($\chi^2_{\text{red,FR}}$) obtained with MFI(FR)+$K$/Ka+PR4 data set shown in Fig. B2. The sky regions where the absolute value of the relative difference $\Delta\beta_{Q(\text{FR}),U(\text{FR})}$ is smaller than $\Delta\beta_{Q,U}$ are correlated to those regions where the $\chi^2_{\text{red,FR}}$ is smaller than $\chi^2_{\text{red}}$ (reddish regions) and vice versa (bluish regions). This result suggests that Faraday rotation might be playing a role in some of the significant differences areas observed between $\beta_s^Q$ and $\beta_s^U$.

## APPENDIX C: FUNCTION-OF-DECLINATION CORRECTION SIMULATIONS

We studied using simulations if the application of a function-of-declination (FDEC) filter to QUIJOTE-MFI maps biases the $\beta_s$ map obtained from component separation. We generated sky simulation maps with the following components at the QUIJOTE-MFI 11 and 13 GHz, K and Ka, and PR4 frequencies:

(i) CMB: Generated as Gaussian random samples using the power spectra obtained from CAMB (Lewis & Challinor 2011) with the latest Planck cosmological parameters (Planck Collaboration 2020e).

(ii) Synchrotron: Generated using the s1 model of the Python Sky Model (PYSM; Thorne et al. 2017).

(iii) Thermal dust: Generated using the d1 model of the PySM.

(iv) Realistic noise simulations: For each experiment, we use the ones described in Section 4.

All components are either generated or downgraded to $N_{\text{side}} = 512$. Then the components maps are added and we apply the corresponding FDEC filter to each signal map. Finally, all maps are downgraded to $N_{\text{side}} = 64$ and smoothed with a Gaussian beam of FWHM = 2 deg following the procedure described in Section 4.

We perform the component separation analysis in two scenarios: (i) when only the QUIJOTE-MFI frequency signal maps are filtered, and (ii) when all maps are filtered. Fig. C1 shows the relative difference (equation 12) between the $\beta_s$ map recovered from the component separation analysis and the $\beta_s$ template (equation 12 taking into account that the uncertainty of the template map is set to zero, $\sigma_{\beta_s} = 0$). We find that when only QUIJOTE-MFI channels are filtered (bottom panel) the relative differences are larger in regions such as the North Polar Spur or the R3 region than when all maps are filtered. Moreover, in those regions the $\beta_s$ relative differences are larger than $3\sigma$ with respect to the template. In the case when all maps are filtered (top panel), these biases are reduced significantly.

<footnote n="1">¹Instituto de Física de Cantabria (IFCA), CSIC-Univ. de Cantabria, Avda. los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="2">²Dpto. de Física Moderna, Universidad de Cantabria, Avda. de los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="3">³Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="4">⁴Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="5">⁵Institut d'Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="6">⁶Astrophysics Group, Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="7">⁷Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="8">⁸Dpto. de Ingenieria de COMunicaciones (DICOM), Edif. Ingenieria de Telecomunicacion, Pl. de la Ciencia s/n, E-39005 Santander, Spain</footnote>

<footnote n="9">⁹Dpto. de Matemáticas, estadística y computación, Univ. de Cantabria, Avda. de los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="10">¹⁰Aurora Technology for the European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="11">¹¹Universidad Europea de Madrid, E-28670, Madrid, Spain</footnote>

<footnote n="12">¹²Jodrell Bank Centre for Astrophysics, Alan Turing Building, Department of Physics and Astronomy, School of Natural Sciences, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="13">¹³Consejo Superior de Investigaciones Científicas, E-28006 Madrid, Spain</footnote>

<footnote n="14">¹⁴Dpto. de Física, Facultad de Ciencias, Univ. de Córdoba, Campus de Rabanales, Edif. C2. Planta Baja, E-14071 Córdoba, Spain</footnote>

<footnote n="15">¹⁵Purple Mountain Observatory, CAS, No.10 Yuanhua Road, Qixia District, Nanjing 210034, China</footnote>

<footnote n="16">¹⁶NAOC-UKZN Computational Astrophysics Center (NUCAC), University of Kwazulu-Natal, Durban 4000, South Africa</footnote>

<footnote>This paper has been typeset from a TeX/LA TeX file prepared by the author.</footnote>