© The Authors 2026

1. Introduction

The increasing sensitivity of Cosmic Microwave Background (CMB) observations has revolutionised our understanding of the Universe in the past three decades. With the tremendous success of three space-based CMB missions, Cosmic Background Explorer (Mather et al. 1990), Wilkinson Microwave Anisotropy Probe (WMAP) (Bennett et al. 2013), and Planck (Planck Collaboration I 2020), and plenty of ground-based and sub-orbital CMB experiments (Carlstrom et al. 2011; Keck Collaboration 2024; Madhavacheril et al. 2024; Ade et al. 2025), we have reached an era where the primary focus of the cosmological observations is to understand the polarisation of the CMB in the upcoming decades. In particular, a hunt for the inflationary gravitational waves through the potential detection of the primordial B-mode polarisation of CMB is one of the key targets of future CMB polarisation missions (Kamionkowski & Kovetz 2016). Future experiments, including LiteBIRD (LiteBIRD Collaboration 2022), Simons Observatory (Ade et al. 2019), CMB-Bharat (Adak et al. 2022), and PICO (Hanany et al. 2019) will provide the most sensitive polarisation maps that will significantly tighten the constraints on the physics of inflation.

The weak signal of B-mode CMB polarisation as compared to the polarised foregrounds makes the experiments very challenging. At frequencies below ∼70 GHz, the dominant foreground contamination is diffuse synchrotron emission generated from gyrating cosmic ray electrons around Galactic magnetic fields (Rybicki & Lightman 1979; Planck Collaboration IV 2020). At frequencies > 100 GHz, polarised thermal emission from dust grains is the key obstacle for B-mode observations (Planck Collaboration IV 2020). Therefore, understanding the spectral and spatial behaviour of these diffuse emissions and subtracting them from the data is an obvious step to reliably map the fluctuation of CMB polarisation. Another potential foreground contamination for B-mode studies is Anomalous Microwave Emission (AME) seen at the frequency range ∼10–60 GHz (Remazeilles et al. 2016; Dickinson et al. 2018; Cepeda-Arroita et al. 2025). The favoured model of AME emission predicts that the radiation is due to electric dipole emission from spinning dust grains of size in the nanometer range (Draine & Hensley 2013; Ali-Haïmoud et al. 2009) populated at different environments of the Interstellar medium (Hoang & Lazarian 2016; Hensley et al. 2022). However, the true origin of AME emission is still ambiguous. Although AME is a well-motivated foreground component for intensity data, the polarisation property of AME is not well understood from available data. The current constraints on AME polarisation are set to below 1% from various observations (Génova-Santos et al. 2017; Herman et al. 2023; Fernández-Torreiro et al. 2024; González-González et al. 2025). Therefore, it is generally believed that AME is not a major obstacle for CMB B-modes. However, it cannot be totally disregarded because: 1) there is still uncertainty on the theoretical models, and 2) the constraints on the polarisation fraction are obtained on specific regions and could be affected by observational effects such as beam depolarisation. The recent studies of Herman et al. (2023) has found inference of the AME polarisation fraction depends on the choice of the prior on the synchrotron spectral index, concluding that the current data sets are not strong enough to simultaneously and robustly constrain both polarised synchrotron emission and AME. In this work, therefore, we investigate using simulations whether AME polarisation biases the synchrotron spectral index measurement using our methodology.

The CMB and various astrophysical foregrounds arise from distinct physical processes and therefore exhibit different spectral behaviours. These spectral differences allow us reliable modelling of the sky emission and enables the extraction of CMB polarisation maps from data obtained by multiple experiments (Eriksen et al. 2008; Stompor et al. 2008; de la Hoz et al. 2023; Brilenkov et al. 2023). The method requires fitting at least six parameters for the case of the simplest polarised sky model. With the sensitivity of current data, the best-fit results at a significant fraction of the sky are degenerate, prior-dependent, and the signal-to-noise ratio (S/N) is low for various foreground parameters (Planck Collaboration XXV 2016; Jew et al. 2019; de la Hoz et al. 2023). For instance, to reduce the inter-dependency among the foreground spectral parameters, in the main analysis of Planck official data using Commander, the synchrotron spectral index has been fixed to the best-fit value derived from the Planck 2015 temperature data, corresponding to β_s = −3.1 (Planck Collaboration IV 2020). Many attempts have been made to exploit the current available data that yield the synchrotron spectral indices ranging between −5 to −1 (see de Belsunce et al. 2022 and references therein), which is inconsistent with predictions from the energy distribution of cosmic ray electrons (Rybicki & Lightman 1979; Yang & Aharonian 2017; Neronov et al. 2017). Therefore, exploiting all pieces of information from independent foreground spectrum analysis is always helpful for parametric component separation to set the priors. For example, synchrotron spectral index estimation from experiments such as S-PASS, C-BASS, QUIJOTE, and CLASS (Krachmalnicoff et al. 2018; Dickinson et al. 2019; de la Hoz et al. 2023; Watts et al. 2024; Eimer et al. 2024) is useful since their channels are synchrotron-dominated. Furthermore, for the recent development of semi-blind component separation methods (Remazeilles et al. 2021; Adak 2021; Carones et al. 2023; Carones & Remazeilles 2024), prior knowledge of the spatial distribution of the spectral index is a piece of very useful information to deproject a few moments (Chluba et al. 2017; Vacher et al. 2023) of the synchrotron to reduce residual leakage to recovered CMB maps.

In this work, we establish a new semi-blind approach to map the synchrotron spectral index distribution using low-frequency data from WMAP K and Ka bands (Bennett et al. 2013), Planck LFI 30 GHz data (Planck Collaboration II 2020), and QUIJOTE MFI (Rubiño-Martín et al. 2023), where the synchrotron is the dominant diffuse emission. For this endeavour, we use a generalised method of constrained Internal Linear Combination (cILC, Remazeilles et al. 2011) for spin-2 fields (Fernández-Cobos et al. 2016; Adak 2021). The cILC is designed to estimate the map of the desired component from a combination of multi-frequency data, deprojecting other emissions. In this work, we aim to reconstruct the synchrotron moment maps by applying this technique and estimate the average spectral index from the correlation of the first two moments. As illustrated in Figure 1, the first two synchrotron moments exhibit distinct spectral energy distributions. Consequently, frequency channels below 40 GHz, where thermal dust emission remains subdominant compared to the synchrotron moments, are expected to aid in disentangling the two moments with this method. In literature, the cILC method has been extensively employed on spin-0 fields for reconstruction of various cosmological signals (Remazeilles et al. 2011; Remazeilles & Chluba 2020; Remazeilles et al. 2021; Rotti & Chluba 2020; Carones et al. 2023; Carones & Remazeilles 2024). The straightforward generalisation of this method for spin-2 fields is to decompose the Q, U maps to E and B-mode polarisation and apply cILC on E and B-mode maps separately. However, the decomposition of the Q, U maps to E, B-mode polarisation maps on incomplete sky data leads to leakage of E to B-mode. This motivates us to apply the constrained polarisation ILC (cPILC hereafter, Adak 2021) directly on Q, U maps instead. The cPILC method was developed in Adak (2021), which is a generalisation of the polarisation ILC (PILC) method developed in Fernández-Cobos et al. (2016) for the case of more than one constraint (Remazeilles et al. 2011). As we are working with Stokes Q and U parameters in Galactic coordinates, one major caveat of the cPILC method is that since it simultaneously minimises Q and U variance, the method is sub-optimal for U mode if the variance is largely driven by Q in general. However, in the case of synchrotron moment recovery, since U has a significant contribution in the variance, the reconstruction is expected to be reasonably good. Furthermore, cPILC is designed to preserve the coherence between the two spinorial components by estimating weights to minimise the variance of a spin-0 field composed of a combination of Q and U (Fernández-Cobos et al. 2016). On the contrary, recovery of synchrotron moments independently in Q and U requires multiplying the Stokes parameters by different weights. Since Q and U depend on the local coordinate frame, this implies that weights are determined by changing arbitrarily the polarisation angle and modulus, spoiling the physical meaning of the spin-2 fields. Therefore, in this work, we employ the cPILC method for our purpose.

Fig. 1. Polarisation intensity, $P=\sqrt{Q^2+U^2}$, of the synchrotron zeroth moment (dash-dotted blue line) and the first moment (green) for two pivot values of βs = −3.0 (dash-dotted green line) and −2.97 (dashed green lines) at a pivot frequency, νs = 22.8 GHz. Emissions from dust, CMB, and AME (0.5% polarisation fraction) are presented by black, cyan, and magenta dash-dotted lines, respectively. The grey areas represent the frequency bands of QUIJOTE MFI bands, WMAP K and Ka-band, and Planck LFI 30 GHz.

This paper is organised as follows. We explain the methodology to estimate the synchrotron moments and inference of spectral index using their correlation in Section 2. In Section 3, we describe the sky models of the simulation used in the validation of the proposed method. In Section 4, we describe the data sets used in this work. We validate the method on simulated sky maps and present the inferred results in Section 5. In Section 6, we provide the inferred spectral index map at N_side = 32 at low and intermediate latitudes where the S/N of synchrotron in data is high and compare our results with the previous analyses in literature. We repeat a similar analysis over the patches defined in Fuskeland et al. (2014) that include regions of low S/N located at high latitudes and discuss the obtained results in Section 7. We conclude the results in Section 8.

2. Methodology

2.1. Estimation of synchrotron moment maps

The Galactic foreground emissions are considered to be the superposition of emissions from many emission blocks along the line-of-sight (LOS). For the synchrotron and thermal dust, the emission follows a power-law and modified blackbody spectrum (MBB), respectively, with variable spectral parameters along the LOS. Therefore, the integrated total emission along the LOS and within the beam size alters the spectral properties of foregrounds, e.g., the sum of multiple power laws is not a power law, and likewise, the sum of multiple MBB spectra is not an MBB. Considering there is an infinite number of emission blocks along LOS, foreground emissions are statistically expressed as a summation of infinite moments following the moment expansion method of Chluba et al. (2017). The synchrotron Stokes parameters $Q_{\nu }^{\mathrm{sync}}(p),U_{\nu }^{\mathrm{sync}}(p)$ at frequency ν can be modelled with an amplitude Q_νs(p),U_νs(p) at reference frequency ν_s and a statistical average of SEDs along LOS and inside the beam, f_sync(ν, β_s(p)) = ∫(ν/ν_s)^β_s(p(r))𝒫(β_s(p(r)))dr,

$$\begin{array}{c}\hfill \left[\begin{array}{c}{Q}_{\nu }^{\mathrm{sync}}(p)\\ U_{\nu }^{\mathrm{sync}}(p)\end{array}\right]=\left[\begin{array}{c}{Q}_{{\nu }_s}(p)\\ U_{{\nu }_s}(p)\end{array}\right]f_{\mathrm{sync}}(\nu ,{\beta }_s(p)),\end{array}$$(1)

where 𝒫 denotes the probability density of having the spectral index β_s(p(r)) at pixel p along the LOS distance r. Following the moment expansion of synchrotron emission in Chluba et al. (2017), we can express Equation (1) as a summation of moments around some pivot value of the spectral index, ${\overline{\beta }}_s$, as

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{Q}_{\nu }^{\mathrm{sync}}(p)\\ U_{\nu }^{\mathrm{sync}}(p)\end{array}\right]& =\left[\begin{array}{c}{Q}_{{\nu }_s}(p)\\ U_{{\nu }_s}(p)\end{array}\right]f_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)\hfill \\ \hfill & +\left[\begin{array}{c}{Q}_{{\nu }_s}(p)\\ U_{{\nu }_s}(p)\end{array}\right]\mathrm{\Delta }{\beta }_s(p){\partial }_{{\beta }_s}{f}_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)+\mathcal{O}({\beta }_s^2),\hfill \end{array}$$(2)

where $\mathrm{\Delta }{\beta }_s(p)={\beta }_s(p)-{\overline{\beta }}_s$, and

$$\begin{array}{cc}\hfill f_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)& ={\left(\frac{\nu }{{\nu }_s}\right)}^{{\overline{\beta }}_s},\hfill \end{array}$$(3a)

$$\begin{array}{cc}\hfill {\partial }_{{\overline{\beta }}_s}{f}_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)& =ln\left(\frac{\nu }{{\nu }_s}\right)f_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s).\hfill \end{array}$$(3b)

Here we do not consider the higher order moments assuming their contribution can be significantly reduced with a proper choice of pivot value of ${\overline{\beta }}_s$ (Rotti & Chluba 2020; Vacher et al. 2023). In principle, since spatial distribution of the synchrotron spectral index is non-Gaussian, achieving high accuracy in synchrotron modelling using moment expansion requires adding more moments, especially where the spectral index deviates further from its pivot value. However, this also requires deprojection of higher-order moments to reduce bias in the moments we aim to recover. Given the sensitivities of the data used here, this would introduce a noise penalty. Therefore, we restrict our moment expansion up to the first order.

Similarly, we can model the thermal dust emission with moments expressed with respect to the pivot values of dust temperature, ${\overline{T}}_d$, and spectral index, ${\overline{\beta }}_d$, as

$$\begin{array}{cc}\hfill [Q_{\nu }^{\mathrm{dust}}(p),U_{\nu }^{\mathrm{dust}}(p)]& =[Q_{{\nu }_d}(p),U_{{\nu }_d}(p)]f_{\mathrm{dust}}(\nu ,T_d(p),{\beta }_d(p))\hfill \\ \hfill & =[Q_{{\nu }_d}(p),U_{{\nu }_d}(p)]f_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d)\hfill \\ \hfill & +[Q_{{\nu }_d}(p),U_{{\nu }_d}(p)]\mathrm{\Delta }{\beta }_d(p){\partial }_{{\beta }_d}{f}_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d)\hfill \\ \hfill & +[Q_{{\nu }_d}(p),U_{{\nu }_d}(p)]\mathrm{\Delta }{T}_d(p){\partial }_{T_d}{f}_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d)\hfill \\ \hfill & +\dots ,\hfill \end{array}$$(4)

where $\mathrm{\Delta }{\beta }_d(p)={\beta }_d(p)-{\overline{\beta }}_d$, $\mathrm{\Delta }{T}_d(p)=T_d(p)-{\overline{T}}_d$, and

$$\begin{array}{cc}\hfill f_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d)& ={\left(\frac{\nu }{{\nu }_d}\right)}^{{\overline{\beta }}_d+1}\frac{e^{{\overline{x}}_d}-1}{e^{\overline{x}}-1},\hfill \\ \hfill {\partial }_{{\beta }_d}{f}_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d)& =ln\left(\frac{\nu }{{\nu }_d}\right)f_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d),\hfill \\ \hfill {\partial }_{T_d}{f}_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d)& =\frac{1}{{\overline{T}}_d}[\frac{\overline{x}{e}^{\overline{x}}}{e^{\overline{x}}-1}-\frac{{\overline{x}}_d{e}^{{\overline{x}}_d}}{e^{{\overline{x}}_d}-1}]f_{\mathrm{dust}}(\nu ,{\overline{\beta }}_d,{\overline{T}}_d),\hfill \\ \hfill \end{array}$$(5)

are the spectra of different moments. Here, $x_d=\frac{h{\nu }_d}{K_B{T}_d}$ and $\overline{x}=\frac{h\nu }{K_B{\overline{T}}_d}$.

The spin-2 fields, P^±(p) = Q(p)±iU(p) of the microwave polarisation data at frequency ν can be described as

$$\begin{array}{cc}\hfill P_{\nu }^{\pm }(p)& =P_{{\nu }_s}^{\pm }(p)f_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)+P_{{\nu }_s}^{\pm }(p)\mathrm{\Delta }{\beta }_s(p){\partial }_{{\overline{\beta }}_s}{f}_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)+n_{\nu }^{\pm }(p)\hfill \\ \hfill & ={\displaystyle \sum _{c=1}^{N_c}}{A}_{\nu c}{P}_c^{\pm }(p)+n_{\nu }^{\pm }(p),\hfill \end{array}$$(6)

where A_νc is the mixing matrix associated with the components $P_c^{\pm }(p)$, N_c is the number of components and $n_{\nu }^{\pm }(p)$ comprises the total nuisance emission from CMB, dust moments, AME and instrument noise. In order to take into account the effect of instrument bandpass in real data, we perform a bandpass integration over the respective spectrum of the components while computing the mixing matrix in real data analysis. This prescription is complementary of multiplying colour correction factors to the data. In Equation (6), we can consider the spin-2 fields of synchrotron amplitude or zeroth moment, P_νs^±(p), and first moment, P_νs^±(p)Δβ_s(p) to be two independent components following the spectrum, $f_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)$ and ${\partial }_{{\overline{\beta }}_s}{f}_{\mathrm{sync}}(\nu ,{\overline{\beta }}_s)$ respectively. In Figure 1, we display polarisation intensity $P=\sqrt{Q^2+U^2}$ of synchrotron zeroth moment (dash-dotted blue line) and first moments for two choice of pivot spectral parameter, ${\overline{\beta }}_s=-3.0$ (dash-dotted green line) and −2.97 (dashed green line) and pivot frequency ν_s = 22.8 GHz. We also present other nuisance emission from CMB (cyan), thermal dust (black) and AME (magenta) with 0.5% polarisation fraction.

In this work, our first step is to estimate the zeroth moment or synchrotron amplitude, P_νs^±(p) and the first moment which is P_νs^±(p)Δβ_s(p). Since the first two moments of synchrotron emission are highly correlated but follows two different spectrum as shown in Figure 1, in order to deproject the possible contamination from one component while estimating another as well as other nuisance foregrounds, we apply the cPILC method as motivated in Section 1. Here we briefly summarise the method.

Assuming the mixing matrix to be constant over some domain 𝒟(p)¹ of the sky for a reasonable choice of the pivot values of the synchrotron and dust spectral parameters, Equation (6) for the data vector ${\mathbf{P}}^{\pm }(p)=\{P_{\nu }^{\pm }\}(p)$ at all channels can be generalised as

$$\begin{array}{c}\hfill {\mathbf{P}}^{\pm }(p)=\mathbf{A}{\mathit{P}}_c^{\pm }(p)+{\mathbf{n}}^{\pm }(p).\end{array}$$(7)

The estimation of the desired component P_c^± deprojecting others in the cPILC method consists of a weighted sum of available channel maps,

$$\begin{array}{c}\hfill {\widehat{P}}_c^{\pm }(p)={\mathit{w}}^T{\mathbf{P}}^{\pm }(p).\end{array}$$(8)

The weights w = {w_ν} are determined by minimising the covariance matrix, C_νν′ = ⟨Q_ν(p)Q_ν′(p)+U_ν(p)U_ν′(p)⟩ estimated using the all pixels within domain 𝒟(p) while the weights have unit response to the spectrum of the component of interest and zero response to the spectrum of the components desired to deproject. In particular, for estimation of first synchrotron moment in Equation (6), the conditions applied to the weights are:

$$\begin{array}{cc}\hfill {\mathit{w}}^T{\partial }_{{\overline{\beta }}_s}{\mathit{f}}_{\mathrm{sync}}& =1,\hfill \end{array}$$(9a)

$$\begin{array}{cc}\hfill {\mathit{w}}^T{\mathit{f}}_{\mathrm{sync}}& =0,\hfill \end{array}$$(9b)

$$\begin{array}{cc}\hfill {\mathit{w}}^T{\mathit{f}}_{\mathrm{dust}}& =0.\hfill \end{array}$$(9c)

In order to minimise possible contamination from thermal dust, we deproject its zeroth moment by applying the constraint given in Equation (9c). Introducing additional constraints can further reduce contamination from higher-order moments, but this comes at the cost of increased noise in the recovered map (Remazeilles et al. 2021). Hence, the number of imposed constraints must be chosen carefully to balance the trade-off between residual foreground contamination and residual noise amplification. In our specific case, the higher-order moments of the synchrotron component are several orders of magnitude smaller than the first two moments, and thermal dust contributes negligibly within the frequency range considered. Therefore, we do not impose additional constraints to deproject higher-order moments of either component. The condition of minimum covariance of C_νν′ with a set of constraints in Equations (9a), (9b), and (9c) using the Lagrange’s undetermined multipliers method yields the optimised weights,

$$\begin{array}{c}\hfill {\mathit{w}}^T={\mathbf{e}}^T{({\mathbf{F}}^T{\mathbf{C}}^{-1}\mathbf{F})}^{-1}{\mathbf{F}}^T{\mathbf{C}}^{-1},\end{array}$$(10)

where $\mathbf{F}=[{\partial }_{{\overline{\beta }}_s}{\mathit{f}}_{\mathrm{sync}},{\mathit{f}}_{\mathrm{sync}},{\mathit{f}}_{\mathrm{dust}}]$ and e = [1, 0, 0]^T. The minimum variance estimated first synchrotron moment is then given by

$$\begin{array}{c}\hfill {\widehat{s}}_1^{\pm }(p)=P_{{\nu }_s}^{\pm }(p)\mathrm{\Delta }{\beta }_s(p)+{\mathit{w}}^T{\mathbf{n}}^{\pm }(p),\end{array}$$(11)

which is an unbiased estimation of the first synchrotron moment with zero contamination of synchrotron amplitude or zeroth moment and some residual contamination from remaining nuisance components in Equation (6). Similarly, if we interchange the constraints in Equations (9a) and (9b), we can obtain the unbiased estimation of the zeroth synchrotron moment,

$$\begin{array}{c}\hfill {\widehat{s}}_0^{\pm }(p)=P_{{\nu }_s}^{\pm }(p)+{\overline{\mathit{w}}}^T{\mathbf{n}}^{\pm }(p).\end{array}$$(12)

From Equations (11) and (12), we finally estimate the map of spectral index modulated synchrotron amplitude, P_νs^±β_s(p) following

$$\begin{array}{c}\hfill {\widehat{s}}^{\pm }(p)={\widehat{s}}_1^{\pm }(p)+{\overline{\beta }}_s{\widehat{s}}_0^{\pm }(p)=P_{{\nu }_s}^{\pm }{\beta }_s(p)+{\widehat{\mathbf{w}}}^T{\mathbf{n}}^{\pm },\end{array}$$(13)

where $\widehat{\mathbf{w}}=\mathit{w}+{\overline{\beta }}_s\overline{\mathit{w}}$. We derive the Stokes parameter maps, Q(p) and U(p), from their corresponding recovered spin-2 fields, ${\widehat{s}}_0^{\pm }(p)$ and ${\widehat{s}}^{\pm }(p)$, obtained from Equations (12) and (13), respectively. For convenience, we denote these maps as m₀(p) = {Q(p),U(p)} and m₁(p) = {Q(p)β_s(p),U(p)β_s(p)}, and we will use these notations in the following sections.

2.2. Estimation of spectral index map

The derived Stokes parameters corresponding to two moment maps in Equations (12) and (13) are related as

$$\begin{array}{c}\hfill m_1(p)+n_1(p)={\widehat{\beta }}_s{m}_0(p)+n_0(p),\end{array}$$(14)

where the corresponding noise vector are denoted by n₀(p) and n₁(p) respectively. In Equation (14), we assume both Q and U follow the same average spectral index, ${\widehat{\beta }}_s$ over some region. We perform the linear regression between two moments using the ‘T-T plot’ technique (Turtle et al. 1962) over the sub-domain 𝒮(p) ∈ 𝒟(p) to estimate the average spectral index within 𝒮(p)². It is worth noticing in Equation (14) (also in Equations (12) and (13)) that some residual leakage from nuisance components and instrument noise contaminates both estimated moments. Therefore, to propagate the uncertainties in both dimensions in the spectral index estimation, we adopt the ‘effective variance method’ (Orear 1982; Petrolini 2014) that supports the standard deviation in both variables, minimising the error function,

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{p\in \mathcal{S}(p)}}\frac{{(m_1(p)-{\widehat{\beta }}_s{m}_0(p))}^2}{{\sigma }_1{(p)}^2+{(\partial m_1(p)/\partial m_0(p))}^2{\sigma }_0{(p)}^2},\end{array}$$(15)

where σ₀(p) and σ₁(p) are the vector of noise standard deviation corresponding to n₀(p) and n₁(p) respectively. The maximum contribution of thermal dust and CMB is at WMAP Ka-band. Their typical level is 1–2% of that of synchrotron. At channels having frequencies lower than WMAP Ka-band, this level is several times smaller. In the recovered maps, the dominant source of total residual leakage is instrumental noise. Consequently, when estimating the standard deviation vectors σ₀(p) and σ₁(p) in Equation (15), we account only for the contribution from noise. Specifically, we compute these vectors from multiple noise residual maps obtained by applying the same optimised weights to 100 independent noise realisations. The contribution from residual leakage of other nuisance components is therefore neglected.

In order to minimise the systematic uncertainties related to the orientation angle used to pixelise the Stokes parameters with respect to the Galactic coordinate system, following Fuskeland et al. (2021), we marginalise over polarisation angles. That is, we perform the T-T plot for Stokes parameters at a new coordinate system by rotating the plane of polarisation by an angle α that ranges from 0° to 85° in bins of 5°. The final inferred spectral index is the inverse-variance-weighted mean of spectral index for all α, i.e., ${\widehat{\beta }}_s={\displaystyle \sum _{\alpha }}\frac{{\widehat{\beta }}_{s,\alpha }/{\sigma }_{{\beta }_{s,\alpha }}^2}{1/{\sigma }_{{\beta }_{s,\alpha }}^2}$ and the error bar is the quadratic sum of systematic uncertainty, $[max{\widehat{\beta }}_{s,\alpha }-min{\widehat{\beta }}_{s,\alpha }]/2$ and statistical uncertainty (Watts et al. 2024). The minimum uncertainty found for different α rotations is considered to be the statistical uncertainty at respective super-pixels.

3. Sky simulation

The polarised microwave sky is primarily composed of CMB, synchrotron, and thermal dust emission. The polarisation fraction of AME is limited below 1% as suggested in various observational limits that are discussed in Section 1. We simulate these components using PySM (Thorne et al. 2017; Zonca et al. 2021) at N_side = 64 and expressed at Rayleigh–Jeans (μK_RJ) units. In simulations, we use the bandpass to be a delta function for all channels. In addition to these sky components, the data also contains the instrument noise. Therefore, we add random realisation of noise maps to the sky simulations at respective channels. The simulation of noise realisation maps for different experiments is discussed in Section 4. In the following sections, we discuss the model of different sky components used in our simulation.

3.1. CMB

We use a Gaussian random realisation of the lensed CMB map simulated using a set of C_ℓs. The C_ℓs are generated from CAMB (Lewis et al. 2000) for best-fit cosmological parameters from Planck 2018 results (Planck Collaboration VI 2020). We ignore the primordial B-mode polarisation in CMB map by setting the tensor-to-scalar ratio parameter r = 0.

3.2. Synchrotron

The spectral energy distribution (SED) of synchrotron is commonly described in μK_RJ units by a power-law model (Rybicki & Lightman 1991),

$$\begin{array}{c}\hfill \left(\begin{array}{c}{Q}_{\nu }^{\mathrm{sync}}\\ U_{\nu }^{\mathrm{sync}}\end{array}\right)=\left(\begin{array}{c}{Q}_{{\nu }_s}\\ U_{{\nu }_s}\end{array}\right){\left(\frac{\nu }{{\nu }_s}\right)}^{{\beta }_s},\end{array}$$(16)

where Q_νs, U_νs are the synchrotron polarisation templates at a reference frequency of ν_s and β_s is the synchrotron spectral index. We use s1 model of synchrotron of PySM where the synchrotron Q and U maps obtained from WMAP 9-year 23 GHz (Bennett et al. 2013) is used as templates at reference frequency ν_s = 23 GHz. The synchrotron spectral index β_s used in simulation varies across the sky. This β_s map is obtained in Miville-Deschênes et al. (2008) fitting the Haslam 408 MHz and WMAP 23 GHz map. The template maps are smoothed to 5° and small-scale information is added to the final template maps following the method discussed in Thorne et al. (2017).

3.3. Thermal dust

In the microwave frequency range, the thermal dust emission is well described by a MBB in μK_RJ units,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{Q}_{\nu }^{\mathrm{dust}}\\ U_{\nu }^{\mathrm{dust}}\end{array}\right)=\left(\begin{array}{c}{Q}_{{\nu }_d}\\ U_{{\nu }_d}\end{array}\right){\left(\frac{\nu }{{\nu }_d}\right)}^{{\beta }_d+1}\frac{e^{\gamma {\nu }_d}-1}{e^{\gamma \nu }-1},\end{array}$$(17)

where Q_νd, U_νd are the dust polarisation templates at reference frequency ν_d, β_d is the dust spectral index and $\gamma =\frac{h}{K_B{T}_D}$ and T_d is the dust temperature. We use one component dust model, d1 of PySM. The model uses the dust polarisation maps obtained from Commander as a template at reference frequency ν_d = 353 GHz (Planck Collaboration X 2016). The template maps are smoothed to FWHM = 2.6° and small-scale information is added following the method discussed in Thorne et al. (2017). The spatially varying dust temperature T_d and spectral index β_d used in this model are obtained from Commander.

3.4. AME

We simulate the AME maps from a2 AME model of PySM,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{Q}_{\nu }^{\mathrm{AME}}\\ U_{\nu }^{\mathrm{AME}}\end{array}\right)=p^{\mathrm{AME}}{I}_{\nu }^{\mathrm{AME}}\left(\begin{array}{c}cos(2{\gamma }_d)\\ sin(2{\gamma }_d)\end{array}\right),\end{array}$$(18)

where p^AME is the polarisation fraction. We modify the PySM code to set p^AME to 0.1% and 0.5% as motivated in Section 1. The AME intensity $I_{\nu }^{\mathrm{AME}}$ has two components which follow two different emissivities that are obtained from the SpDust2 code (Ali-Haïmoud et al. 2009). The templates of spatially varying two components are obtained from Commander fitting (Planck Collaboration X 2016). The AME intensity maps are translated to polarisation maps using the same dust polarisation angle γ_d obtained from Q,  U dust maps of Commander fitting.

4. Data

We use frequency maps below 40 GHz which can be treated as synchrotron tracers. In this paper, the data products we use are from the following three surveys.

(i) QUIJOTE: The QUIJOTE (Q-U-I JOint TEnerife) MFI survey³ (Rubiño-Martín et al. 2023) has observed approximately 70% of the sky across four frequency bands: 11, 13, 17, and 19 GHz. Among the four channels, the 11 GHz and 13 GHz channels provide the highest S/N to synchrotron emission. However, the noise at 13 GHz is observed to be ∼30% correlated with that at 11 GHz (Rubiño-Martín et al. 2023) for polarisation. Therefore, to maximise the effective S/N, we obtain our main results in our analysis using a weighted combined map of the 11 GHz and 13 GHz channel data. We assign weights of 55% and 45% respectively for 11 and 13 GHz channels based on their sensitivity. In addition to this baseline data combination, we also repeat the full analysis using the 11 GHz and 13 GHz maps individually to assess the impact of each channel on the inferred results. Certain regions of the survey suffer from reduced sensitivity due to factors such as the presence of geostationary satellites and increased atmospheric contamination in specific directions. The S/N is particularly high near the Galactic plane and the North Polar Spur (NPS) regions. To ensure the robustness of our analysis, we restrict ourselves to the high S/N regions that encompass 30% of the sky of the QUIJOTE MFI survey area. This selected area includes prominent features such as the Fan region and the NPS extending north Galactic plane and near the Galactic centre. Additionally, following Watts et al. (2024), we apply the P06 polarisation mask provided by WMAP⁴ to exclude the brightest emission regions. The final mask used in our analysis is shown in Figure 2.

Fig. 2. Mask at Nside = 64 used in this analysis that retains 30% high S/N regions of the QUIJOTE MFI survey. The bright source positions are masked using P06 polarisation mask provided by WMAP.

(ii) WMAP: We consider WMAP 9-year data⁵ at K and Ka bands (Bennett et al. 2013). We do not consider other bands of WMAP data due to the lower S/N of the synchrotron signal.

(iii) Planck: We use Planck low-frequency instrument (LFI) data at 30 GHz from PR3 data release⁶ (Planck Collaboration I 2020; Planck Collaboration II 2020). We verified in Appendix A whether the recovered spectral index remains essentially unchanged when using PR4 (Planck Collaboration LVII 2020) instead.

An important consideration in the use of QUIJOTE MFI data is the possible impact of Faraday rotation on synchrotron emission (Beck et al. 2013). Faraday rotation, which alters the polarisation angle of incoming radiation, becomes increasingly significant at lower frequencies due to its dependence on the square of the wavelength (Beck et al. 2013). This effect can lead to a flattening of the synchrotron spectral index, particularly near the Galactic plane, where the electron density and magnetic fields are stronger. To assess the potential influence of Faraday rotation in our analysis, we examined the spatial correlation between the Stokes Q and U maps from WMAP K-band and QUIJOTE MFI 11 GHz data. We found a high degree of correlation between the two frequency channels, which suggests the presence of a coherent magnetic field largely aligned with the Galactic plane. This also indicates that the level of Faraday depolarisation near the Galactic plane is negligible, which is our region of interest. Furthermore, since Rubiño-Martín et al. (2023) reported in their analysis that there are no changes in the final inference after correcting for Faraday rotation using the rotation measure from Hutschenreuter & Enßlin (2020). Therefore, we do not apply any correction for Faraday rotation in this analysis.

Before applying the cPLIC method, all data are processed as follows:

• (i)

  The residuals of Radio Frequency Interference (RFI) mainly due to emissions from geostationary satellites projected in the map-making process appear at fixed azimuth locations in QUIJOTE MFI frequency maps. This residual RFI is subtracted using a function of the declination (FDEC⁷). This function applies a filter to the maps to remove the median value of all pixels at constant declination. de la Hoz et al. (2023) reported that if we apply the same FDEC correction to the other data sets, the potential bias in β_s gets reduced significantly. Therefore, and in order to maintain the consistency of zero levels to the other data sets, the same FDEC filter is applied to WMAP and Planck data.

• (ii)

  Our analysis requires all the maps to have the same resolution and to be pixelised on the same HEALPix grid (Górski et al. 2005; Zonca et al. 2019). Therefore, we apply the FDEC correction to all frequency maps at the common one degree resolution and N_side = 512, and after this, we further convolve the maps to a final Gaussian beam of FWHM of 2° to increase the signal to noise. Finally, all maps are downgraded to N_side = 64 and converted to μK_RJ units.

We also use the noise realisation maps at these frequencies in order to estimate the variance of residual leakage to recovered synchrotron moment maps using the cPILC method. We consider QUIJOTE MFI simulated noise maps that include the 1/f noise (Rubiño-Martín et al. 2023). We simulate the WMAP noise using Gaussian random realisations for standard deviation maps, $\sigma (p)={\sigma }_0/\sqrt{N_{\mathrm{obs}}(p)}$, where σ₀ is RMS noise at respective frequencies (Hinshaw et al. 2003) and N_obs is the hit map of WMAP mission⁸. For Planck channels, we use noise maps of end-to-end noise simulation (Planck Collaboration XII 2016). All the noise realisation maps are also smoothed to a common beam resolution of FWHM = 2°, downgraded to N_side = 64 and converted to μK_RJ units.

5. Validation on simulated maps

We first apply our method to two different sets of simulations for validation. We start with the simulation without a polarised AME component. We discuss the obtained results analysing this set of simulations in Section 5.1 and Section 5.2. To investigate the impact of AME polarisation, next, we incorporate a polarised AME component into the simulation. We discuss the corresponding results in Section 5.3.

5.1. Reconstruction of synchrotron moment maps

We apply the cPILC method to simulated data to recover the synchrotron moment maps while cleaning the nuisance components described in Equation (6). Specifically, we obtain maps of the synchrotron Stokes parameters, Q_νs(p) and U_νs(p), as well as the spectral index modulated Stokes parameters, Q_νs(p)β_s(p) and U_νs(p)β_s(p), at a reference frequency of ν_s = 22.8 GHz, following Equations (11) and (13), respectively. These two sets of maps have zero contamination from one to another. To minimise the contribution from synchrotron higher-order moments to residual leakage, we carefully select the pivot value of the synchrotron spectral index. This is achieved by iteratively varying the pivot value within the range of −3.2 to −2.6. Figure 3 shows the mean absolute difference between the inferred and input β_s, normalised to the corresponding mean uncertainty σ_βs, as a function of the chosen pivot value. We find that the deviation of the inferred spectral index is smallest when the pivot value is close to the central value of the inferred β_s distribution. The deviation remains nearly unchanged for pivot values between −3.0 and −2.9. Therefore, we adopt an iterative scheme that updates the pivot until the inferred distribution stabilises. Without deprojecting thermal dust emission, we find that the recovered synchrotron moments; particularly the first-order moment maps contain noticeable residuals from thermal dust. To mitigate this, we impose an additional constraint as described in Equation (9c) to deproject the zeroth moment of thermal dust. We adopt the Planck best-fit values of thermal dust spectral parameters at ${\overline{T}}_d=19.8$K and ${\overline{\beta }}_d=1.55$ (Planck Collaboration IV 2020) to be the pivot values of the MBB parameters while deprojecting the dust moment. In Figure 4, we compare the input (upper panel) and cPILC-recovered (lower panel) maps of Q_νs(p),U_νs(p) and Q_νs(p)β_s(p),U_νs(p)β_s(p) at the reference frequency of 22.8 GHz. All the maps are displayed at a resolution of 2°. The method enables us to recover both the moments with high fidelity, as evidenced by the strong spatial correlation between the input and recovered maps. The recovered Q maps appear cleaner than the U maps, reflecting the fact that the variance minimisation in the cPILC method is dominated by the stronger Q signal. Finally, we examine the correlation between the two recovered moment maps. Figure 5 shows the T–T correlation between them, yielding correlation coefficients of 0.99 for Stokes Q and 0.95 for Stokes U. The black dashed line indicates a linear fit with β_s = −2.96, corresponding to the mean value of the synchrotron spectral index template used in the simulation.

Fig. 3. Changes of the mean absolute deviation of the inferred spectral index with respect to the input synchrotron spectral index template, normalised with respect to the mean error of estimated βs with a different choice of the pivot value of the synchrotron spectral index.

Fig. 4. Comparison between input (upper panel) and recovered (lower panel) maps of Qνs(p), Uνs(p) and Qνs(p)βs(p), Uνs(p)βs(p) at the reference frequency of 23 GHz. All maps are shown at resolution of beam FWHM = 2°. The correlation between input and recovered maps is found to be 99% for Q and 95% for U. A similar correlation is found for Qνs(p)βs(p) and Uνs(p)βs(p) between input and recovered maps.

Fig. 5. Correlation between recovered Stokes Q(p) (green) and U(p) (magenta) maps corresponding to zeroth and first order moments. The black dashed line is the linear fit with βs = −2.96 which corresponds to mean value of the input spectral index map.

5.2. Inference of synchrotron spectral index

The two synchrotron moment maps recovered in the previous section are found to be anti-correlated. We exploit this anti-correlation to infer the synchrotron spectral index, β_s, by performing a linear regression (“T–T plot”) between the two derived synchrotron moment maps, as described in Section 2.2. The regression is carried out over sub-pixels at N_side = 64 contained within each super-pixel at N_side = 32. We estimate the β_s for 18 rotation angles within 0–85° and then compute the inverse-variance-weighted mean of the spectral index from measurements of all rotation angles. The resulting map of spectral index at N_side = 32 is presented in the middle panel of Figure 6. We present the input reference β_s map in the left panel of this figure for comparison. The right panel presents the corresponding uncertainty map, ${\sigma }_{{\beta }_s}^{\mathrm{stat}+\mathrm{sys}}$, which combines both statistical and systematic errors (stat + sys) in quadrature, as discussed in Section 2.2. Visually, the recovered spectral index map closely resembles the input β_s map. However, a small level of disagreement is expected. The PySM pipeline (Thorne et al. 2017) uses the synchrotron spectral index template from Miville-Deschênes et al. (2008), smoothed to 5°. Additional small-scale fluctuations in the simulated frequency maps are introduced later through the power spectrum, as outlined in Section 3.2. The final simulated maps are smoothed to 2° in our post-processing before applying the cPILC method. These non-linear processes naturally lead to minor discrepancies between the recovered and input spectral index maps. For a more quantitative assessment, the upper panel of Figure 7 shows the histogram of the normalised deviation of the inferred spectral index, $({\beta }_s^{\mathrm{Inferred}}-{\beta }_s^{\mathrm{Input}})/{\sigma }_{{\beta }_s}^{\mathrm{stat}+\mathrm{sys}}$ (orange). If the estimation is unbiased and the uncertainties are well characterised, this distribution should follow a normal function, 𝒩(0, 1), represented by the solid black line. We find that, across a large fraction of the analysed region, the inferred β_s values are consistent with the input within 3σ confidence. Only about 13% of our analysis region exhibits lower confidence, primarily due to the reduced S/N in those regions. To further evaluate the performance, we examine the spatial correlation between the inferred and input β_s maps. The lower panel of Figure 7 presents the pixel-by-pixel correlation between the two. We find a strong correlation of 81% between the inferred and input β_s maps. Within the analysed region, the input β_s distribution has a mean of −2.96 and a standard deviation due to intrinsic scatter of 0.06, while the corresponding inferred distribution yields a mean of −2.97 and a standard deviation of 0.16.

Fig. 6. Left panel: The input βs map used in the simulation. Middle panel: The inferred βs map from the estimated synchrotron moments. Right panel: The map of statistical and systematic uncertainty of the estimated βs which are added in quadrature.

Fig. 7. Upper panel: Histogram of the normalised deviation of the derived spectral index with respect to the true input spectral index $({\beta }_s^{\mathrm{Inferred}}-{\beta }_s^{\mathrm{Input}})/{\sigma }_{{\beta }_s}^{\mathrm{stat}+\mathrm{sys}}$ (orange). The black curve is a Gaussian distribution with zero mean and unit variance. Vertical lines are the corresponding mean of the distribution. Lower panel: The two-dimensional density plot of true input and inferred βs map. The correlation between the two maps is 81%. The black dashed line correspond to 100% correlation.

5.3. Spectral index inference in presence of polarised AME

The polarisation fraction of diffuse AME remains relatively uncertain. As noted in Section 1, several previous theoretical and empirical studies have suggested that AME can be polarised at level below 1% (Draine & Hensley 2013; Hoang & Lazarian 2016; Génova-Santos et al. 2017; González-González et al. 2025). In this section, we investigate the impact of polarised diffuse AME on the estimation of the synchrotron spectral index. We perform our analysis using two sets of simulations based on the PySMa2 model, incorporating polarised AME with polarisation fractions of 0.5% and 0.1%, respectively. It is important to emphasise that we do not impose any extra constraint in our cPILC method to deproject the AME moments while recover the synchrotron moments, which would increase the residual noise level in the recovered maps and lead to greater dispersion in the inferred spectral index. Figure 8 illustrates the comparison of the normalised spectral index deviations for three cases: simulations without polarised AME (orange), with 0.5% polarised AME (magenta), and with 0.1% polarised AME (green). Our results indicate that the presence of a 0.5% polarised AME component causes the inferred synchrotron spectral index to appear slightly flatter compared to the unpolarised case. The mean and standard deviation of the spectral index in this scenario are −2.96 and 0.16, respectively, indicating that the mean value is shifted by only 0.5% relative to the unpolarised AME results. The β_s values inferred from the polarised and unpolarised AME cases are consistent within 1.5σ across all pixels. When the AME polarisation fraction is reduced to 0.1%, the mean and standard deviation of the inferred spectral index become −2.95 and 0.15, respectively, which are statistically indistinguishable from the unpolarised AME case, showing no measurable bias.

Fig. 8. Histograms of the normalised deviation of the spectral index derived from simulations of three cases: without polarised AME (orange), with 0.5% (magenta), and 0.1% (green) polarised AME. The black curve is a unit variance Gaussian distribution with zero mean. Vertical lines are the corresponding mean of the distributions.

6. Results on real data

We now proceed to apply the component separation method to the 2° smoothed maps at frequencies < 40 GHz from Planck PR3 LFI 30 GHz, WMAP K and Ka bands and QUIJOTE MFI 11 and 13 GHz channels. The results obtained with Planck PR4 data at 30 GHz are presented in Appendix A. For the QUIJOTE MFI data, the 11 GHz and 13 GHz channels are combined using weights of 55% and 45%, respectively, to avoid accounting for the correlated noise between two channels as stated in Section 4. This choice of the weights are not strictly optimal, instead they are chosen to take into account the higher S/N at 11 GHz as compared to 13 GHz. In addition to this analysis with combined MFI maps, we have also used the 11 GHz and 13 GHz maps separately in cPILC algorithm and perform subsequent analysis. The corresponding results are discussed later in this section. The constraint equations given in Eqs. (9a), (9b), and (9c) are imposed on the effective spectral response of each component, computed integrating the instrument-specific bandpass profiles (Jarosik et al. 2003; Planck Collaboration V 2014; Rubiño-Martín et al. 2023) over the spectral behaviour of respective components. Given that the WMAP 9-year map-making process incorporates bandpass shift, we adopt a similar treatment of bandpass variation in our analysis, as detailed in Bennett et al. (2013). For the combined QUIJOTE MFI data, we applied the same channel weights in the bandpass integration as those used to produce the combined maps. From the component-separated maps, we reconstruct the Stokes parameter maps corresponding to synchrotron amplitude and spectral index modulated synchrotron amplitude following Eq. (12) and Eq. (11) at N_side = 64. We display the resulting maps in Figure 9. The figure demonstrates a clear spatial anti-correlation between the zeroth and first order recovered moment maps, a feature that is further utilised in the estimation of the synchrotron spectral index.

Fig. 9. Synchrotron Stokes Qνs (first panel) and Uνs (second panel) maps at 22.8 GHz recovered from the cPILC component separation method using the QUIJOTE MFI 11 and 13 GHz combined map, WMAP K and Ka bands and Planck LFI 30 GHz. The Stokes parameters corresponding to the spectral index modulated synchrotron amplitudes, Qνsβs, Uνsβs obtained from component separation are shown in the third and fourth panels respectively. All maps are at FWHM = 2° at Nside = 64.

Subsequently, we perform a T-T correlation analysis at N_side = 32 between pairs of recovered Stokes maps corresponding to the two reconstructed synchrotron moments. To ensure the robustness of the correlation, we exclude some super-pixels from our analysis for which the Pearson correlation coefficient between the moment maps falls below 0.3. This prescription resulted in the removal of 0.5% of the total pixels from the analysis. We compute the inverse-variance-weighted mean and the uncertainty, ${\sigma }_{{\beta }_s}^{\mathrm{sys}+\mathrm{stat}}$, of the spectral indices, following Section 2.2, from the estimated spectral index for 18 rotation angles (α ∈ [0° ,85° ]). In the upper panel of Figure 10, we present the inverse-variance-weighted mean map of spectral index at N_side = 32. In the lower panel, we show the corresponding uncertainty map ${\sigma }_{{\beta }_s}^{\mathrm{sys}+\mathrm{stat}}$.

The inverse-variance weighted mean and uncertainty of the synchrotron spectral index within the Galactic plane region (|b|< 10°) is found to be ${\beta }_s^{\text{plane}}=-3.05\pm 0.01$, while at high Galactic latitudes which is dominated by the NPS for the choice of our mask, it steepens to ${\beta }_s^{\text{high-lat}}=-3.13\pm 0.02$. This behaviour is consistent with earlier findings in the literature (Fuskeland et al. 2014; Krachmalnicoff et al. 2018; Fuskeland et al. 2021), which suggest that the spectral index becomes steeper at high latitudes due to reduced depolarisation effects compared to those prevalent near the Galactic plane. It is important to note that the quoted uncertainties represent only the estimation errors arising from instrumental noise and systematics; they do not include the intrinsic dispersion of the spectral index, which we describe below. In Figure 11, we present a comparison of the spectral index distribution derived in this work (orange) with those reported by de la Hoz et al. (2023) (green) and Miville-Deschênes et al. (2008) (blue). The result in de la Hoz et al. (2023) is based on fitting a parametric model to data from QUIJOTE MFI, WMAP, and Planck, including frequency channels up to 353 GHz to improve the modelling of thermal dust emission. In contrast, our new method adopts a semi-blind strategy combined with a moment expansion approach, that makes fewer assumptions about SEDs of foregrounds and reduces reliance on detailed foreground parametrisation. Furthermore, we exclude the higher-frequency channels (ν > 40 GHz) to minimise contamination from dust emission and instead focus on the low-frequency regime (ν < 40 GHz) where synchrotron emission dominates. Nevertheless, both independent approaches yield consistent results when applied to the data. The spectral index β_s derived in Miville-Deschênes et al. (2008) is widely adopted as a reference template for the synchrotron emission modelling and is incorporated into the standard PySM models. This β_s map was estimated using 408 MHz Haslam map (Haslam et al. 1982) and WMAP 23 GHz data using a model based on a Galactic magnetic field (‘Model 4’). However, both the spectral index distributions obtained in this work and by de la Hoz et al. (2023) exhibit noticeably larger dispersion compared to that reported by Miville-Deschênes et al. (2008). To quantify this broadening, we compute the inverse-variance-weighted mean and standard deviation of the spectral index for all three β_s maps over the analysed region. The estimated mean and standard deviation due to intrinsic scatter are −2.96 and 0.06 for Miville-Deschênes et al. (2008), −3.11 and 0.21 for this work (a broadening by a factor of ∼3.5), and −3.13 and 0.16 for de la Hoz et al. (2023) (a broadening by a factor of ∼2.7), respectively. A similarly broad variation in spectral index was also reported in Weiland et al. (2022) over comparable sky regions. Additional studies (Fuskeland et al. 2014; Vidal et al. 2015; Krachmalnicoff et al. 2018; Fuskeland et al. 2021) further support the presence of significant spatial variability in the synchrotron spectral index across the sky.

We further conduct the analysis to evaluate the sensitivity of the derived synchrotron spectral index to the specific choice of QUIJOTE MFI input data. We examine how the results vary when the combined MFI map, considered as a baseline, is replaced with the individual frequency maps at 11 and 13 GHz. When we consider only 11 GHz map in the component separation, the resulting spectral index has a mean and standard deviation of −3.09 and 0.23 respectively, indicating a systematic shift of only 2.1% (δβ_s = 0.07) compared to the baseline result obtained with the combined MFI data. The corresponding spectral index and uncertainty maps are shown in the first and second panels of Figure 12, respectively. A visual comparison with Figure 10 reveals that the overall morphology remains largely consistent, indicating that the 11 GHz map alone provides a stable recovery of the synchrotron spectral structure. In contrast, when the 13 GHz map is considered in cPILC, the results show a significant deviation. The derived spectral index and uncertainty maps for this case are displayed in the third and fourth panels of Figure 12. These maps differ substantially from those obtained using either the combined MFI data or the 11 GHz-only case, particularly in high S/N regions such as the Fan region and NPS. The mean and standard deviation of the spectral index in this scenario are found to be −3.24 and 0.41, respectively, with the increased dispersion reflecting the lower S/N of the 13 GHz data. We revisited the simulations to verify the origin of the large intrinsic scatter in β_s when including the 13 GHz channel in the cPILC analysis. Our simulation results yield a mean and standard deviation of −2.95 and 0.27, respectively, representing an increase in dispersion of approximately 68% compared to the case using only the 11 GHz channel. This finding suggests that the analysis based on the 13 GHz channel alone may, in fact, degrade the reconstruction of the moment maps and consequently the reliability of the inferred spectral index due to its lower S/N. On the other hand, the analysis based on 11 GHz alone or combined maps 11 and 13 GHz give consistent results. Therefore, we adopt the joint use of both QUIJOTE MFI channels as default for analysis in the subsequent section.

Fig. 10. Synchrotron spectral index (upper panel) and uncertainty maps (lower panel) obtained from combined the QUIJOTE MFI, WMAP, and Planck data.

7. Analysis over regions

We extend the analysis to the low S/N regions of QUIJOTE MFI survey partitioning the sky into regions defined in Fuskeland et al. (2014). We multiply the QUIJOTE MFI survey mask (Rubiño-Martín et al. 2023) by the mask used in Fuskeland et al. (2014), thereby redefining the effective sky coverage within each subregion. As a result, several regions, specifically regions 11, 22, and 23 as defined in Figure 1 of Fuskeland et al. (2014) are excluded from our analysis. In addition, we apply further masking to exclude some of the brightest and most contaminated areas using P06 polarisation mask of WMAP discussed in Section 4. The final analysis mask is shown in Figure 13. We perform the spectral index estimation in each of these regions following the same methodology used in the previous section. In this analysis, we use the combined map of QUIJOTE MFI 11 and 13 GHz data in our component separation procedure. Subsequently, the T-T correlation between moments are performed over the sub-pixels of each region following same procedure as in Section 2.2. The inverse-variance-weighted mean and uncertainty of the derived spectral index values for each region are reported in Table 1, and the corresponding maps are presented in Figure 14. It is important to note that the synchrotron spectral index varies significantly within these broad regions, as seen in our earlier results. Therefore, the average spectral indices obtained by fitting the synchrotron moments are of limited interpretative value when spatial variability is ignored. This is because the underlying assumption of T–T plotting is that all pixel pairs within a region share a common spectral index, which may not hold across all regions. In the literature, several studies have conducted T–T correlation analyses using two frequency data over similar regions with Planck and WMAP data (Fuskeland et al. 2014, 2021; Watts et al. 2024). Their results are also typically scattered due to such internal variations.

Fig. 11. Distribution of the synchrotron spectral index from ‘Model 4’ of Miville-Deschênes et al. (2008) (blue), de la Hoz et al. (2023) (green) and our estimation (orange). The vertical dashed lines correspond to the weighted mean values of the respective distributions.

Fig. 12. Synchrotron spectral index and corresponding uncertainty maps obtained using the QUIJOTE MFI 11 GHz map alone in replacement of QUIJOTE MFI combined maps of 11 and 13 GHz are shown in the first and second panels, respectively. The similar plots of synchrotron spectral index and corresponding uncertainty maps when we use QUIJOTE MFI 13 GHz data alone are shown in the third and fourth panels respectively.

Fig. 13. Regions used for estimation of the spectral index defined in Fuskeland et al. (2014). The grey regions are masked since the regions are not covered by QUIJOTE MFI survey. Some of the bright sources are masked using P06 polarisation mask provided by WMAP.

Fig. 14. Upper panel: The derived spectral index over the regions. Lower panel: The uncertainty map of the derived spectral index.

We compare our results with the recent findings of the COSMOGLOBE DR1⁹ (Watts et al. 2024) and QUIJOTE (de la Hoz et al. 2023). Watts et al. (2024) performed T–T correlation analyses over comparable sky regions using COSMOGLOBE maps at K-band and 30 GHz (Watts et al. 2023). de la Hoz et al. (2023) performed the parametric fitting using QUIJOTE MFI 11, 13 GHz, WMAP 9-year K, Ka, Planck PR4 data. In Figure 15, we display the spectral index obtained in our analysis (orange), and those found in de la Hoz et al. (2023) (green) and Watts et al. (2024) (black). Since de la Hoz et al. (2023) obtained the results from pixel-by-pixel fitting the data instead of joint fitting to all pixels in the regions of Figure 13, the green points represented in Figure 15 are the inverse-variance weighted average of the spectral indices of de la Hoz et al. (2023) at all pixels enclosed in each region. The corresponding plotted error bars are the weighted standard deviation estimated over the same set of pixels, and this explains why they are typically larger. The differences in results from the three methods might be due not only to differences in the datasets used but to the different methodologies. For instance, the β_s values in parametric fitting approach of de la Hoz et al. (2023) is prior dominated especially in low S/N regions, which can produce noticeable discrepancies. Also, pixel-by-pixel fitting instead of fitting the data over entire region, could cause differences. However, we observe reasonable consistency, especially at low-Galactic latitudes. Next, we investigate the low-Galactic regions which has a large overlapping region between our mask and the mask used in Watts et al. (2024). We find that our estimated spectral indices are in very good agreement with those reported by Watts et al. (2024) in regions 13 (NPS), 14 (Galactic centre), 18, 19 (Fan region), and 20. The rest of the Galactic regions (15, 16, 17, and 21) show some level of disagreement between the two analyses. In these low-Galactic regions, we find that the spectral indices obtained in Watts et al. (2024) exhibit larger uncertainties than ours. The improvement in precision in our analysis arises from the inclusion of low-frequency data from the QUIJOTE MFI survey, which significantly enhances sensitivity to the synchrotron component. In particular, we find steeper spectral indices than Watts et al. (2024) in regions 15, 16, and 24, while in region 17, our estimate yields a slightly flatter index. At higher Galactic latitudes, however, the comparison becomes less straightforward. This is primarily due to the smaller overlap between the regions of patches used in our analysis and those in Watts et al. (2024) at high Galactic latitudes, and second, due to the low S/N data in these regions. As a result, direct comparison of the results in those patches is less reliable. To capture the overall trend, we compute the inverse-variance weighted mean and uncertainty of the spectral index across the Galactic plane (regions 15–25) and at high Galactic latitudes (regions 1–14), obtaining ${\beta }_s^{\mathrm{plane}}=-2.96\pm 0.01$ (stat + sys) and ${\beta }_s^{\mathrm{high}\text{-}\mathrm{lat}}=-3.02\pm 0.02$ (stat + sys), respectively. The value within the Galactic plane is in good agreement with previous results from Fuskeland et al. (2014) and Watts et al. (2024). At high Galactic latitudes, however, our estimate is slightly less steep than those reported in COSMOGLOBE and Fuskeland et al. (2014). This mild difference may be attributed to the different sky coverage and inclusion of low-frequency QUIJOTE MFI data in our analysis, which improves synchrotron sensitivity.

Fig. 15. Synchrotron spectral index as a function of the region number obtained in this work (orange), COSMOGLOBE DR1 (black), and de la Hoz et al. (2023) (green). The data shown in green are the inverse-variance-weighted mean of spectral index and uncertainties computed from βs and σβs maps of de la Hoz et al. (2023) over regions.

8. Summary and conclusion

We have developed a new approach based on constrained-ILC and moment expansion of CMB foregrounds that allows us to recover synchrotron moment maps, enabling the construction of a synchrotron spectral index map through their correlation. This approach is fully developed on map-based semi-blind component separation and differs from parametric fitting of a model to the data. The analysis incorporates 2° smoothed data from QUIJOTE MFI at 11 and 13 GHz, WMAP 9-year data at K and Ka bands, and the Planck PR3 data at LFI 30 GHz channel. These frequencies are dominated by synchrotron emission, making them particularly suitable for our study. We first validate the method using simulations and demonstrate that it can recover the synchrotron spectral index at high fidelity in the regions with high signal-to-noise ratio. The inferred results remain largely unbiased in the presence of polarised AME with a polarisation fraction below the current observational upper limit. Then we apply the method to real data to reconstruct the synchrotron spectral index map at the Galactic plane and the NPS region. Over this region our estimated spectral index map is consistent with the one obtained in de la Hoz et al. (2023) which used a parametric model fitting to a similar data sets. Our analysis differs from that of de la Hoz et al. (2023) in two key aspects. First, we restrict our dataset to channels below 40 GHz, whereas de la Hoz et al. (2023) employed all polarised channels up to 353 GHz to constrain dust parameters. Second, we combine the QUIJOTE MFI 11 GHz and 13 GHz maps since the noise in these two MFI channels is significantly correlated, which can compromise the performance of our semi-blind component separation technique. Instead of the combined QUIJOTE MFI maps, using only the 11 GHz data leaves the results remain approximately the same with only 2.1% drift of mean β_s. However, when we use QUIJOTE MFI 13 GHz data, the obtained spectral index gets largely dispersed due to low S/N at 13 GHz channel. We find the inverse-variance-weighted mean and dispersion of the spectral index distribution are −3.11 and 0.21 respectively with combined QUIJOTE MFI data. This variability is notably larger than that predicted by the widely used synchrotron template from Miville-Deschênes et al. (2008), indicating that this template may underestimate spatial variation in the synchrotron spectral index due to many assumptions in synchrotron modelling and significant smoothing. We find moderately steeper spectral index at high Galactic latitude compared to the Galactic plane, consistent with previous studies in the literature.

We extend our analysis across the full QUIJOTE MFI survey region by dividing the sky into 21 subregions. This analysis yields a inverse-variance-weighted mean and uncertainties of the synchrotron spectral index over all regions of ${\beta }_s^{\mathrm{all}\text{-}\mathrm{region}}=-2.99\pm 0.01$. Considering the Galactic plane and high-latitude regions separately, we find the values are ${\beta }_s^{\mathrm{plane}}=-2.96\pm 0.01$ and ${\beta }_s^{\mathrm{high}\text{-}\mathrm{lat}}=-3.02\pm 0.02$ respectively. These results are consistent with previous estimates reported in Fuskeland et al. (2014) and Watts et al. (2024), which performed T-T correlation analyses over similarly defined sky patches. We also compare our estimates on a patch-by-patch basis with the inverse-variance weighted average of the spectral indices from de la Hoz et al. (2023), and the spectral index values reported in Watts et al. (2024). Although the results from three methods are similar at some low-Galactic high S/N regions, there is a difference in high-Galactic regions due to lower sensitivity of the data and differences in the analysis techniques.

In this analysis, we assume that the synchrotron spectrum follows a simple power-law. If the true synchrotron spectrum exhibits a significant positive curvature, the second-order moment, as described in Remazeilles et al. (2021), would become relevant. However, our current study does not account for such higher-order effects. Accurate estimation of higher-order moments will require data with improved sensitivity across a larger number of low-frequency channels. Future experiments such as QUIJOTE-MFI2 (Hoyland et al. 2022; Lorenzo-Hernández et al. 2024) and LiteBIRD (LiteBIRD Collaboration 2022) will provide the necessary frequency coverage and sensitivity to better constrain the synchrotron curvature, thereby extending the applicability of this technique.

This work presents the first application of a semi-blind component separation technique for estimating maps of foreground spectral indices, offering an alternative to fully parametric model-fitting approaches. Our method mitigates bias from unwanted components by deprojecting their first few moments using partial knowledge of their SEDs, thereby enhancing robustness to model inaccuracies. The advantage of the method proposed here is the use of the moment expansion approach, which makes minimal assumptions on the SED of the foregrounds, as compared to the existing parametric methods.

Data availability

All the synchrotron moment maps and the derived spectral index maps obtained in this work are publicly available in the QUIJOTE web page.. The data in Table 1 is available in electronic form at the CDS.

Acknowledgments

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 Astrofisica de Canarias (IAC), the Instituto de Fisica 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, PGC2018-101814-B-I00, PID2019-110610RB-C21, PID2020-120514GB-I00, IACA13-3E-2336, IACA15-BE-3707, EQC2018-004918-P, PID2023-150398NB-I00, PID2023-151567NB-I00, PID2022-139223OB-C21 and PID2022-140670NA-I00, the Severo Ochoa Programs SEV-2015-0548 and CEX2019-000920-S, the Maria de Maeztu Program MDM-2017-0765, and by the Consolider-Ingenio project CSD2010-00064 (EPI: Exploring the Physics of Inflation). We acknowledge support from the ACIISI, Consejeria de Economia, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference ProID2020010108, and Red de Investigación RED2022-134715-T funded by MCIN/AEI/10.13039/501100011033. This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement number 687312 (RADIOFOREGROUNDS), and the Horizon Europe research and innovation program under GA 101135036 (RadioForegroundsPlus). MFT acknowledges support from the Enigmass+ research federation (CNRS, Université Grenoble Alpes, Université Savoie Mont-Blanc). FP acknowledges support from the MICINN under grant numbers PID2022-141915NB-C21. EMG acknowledges support from the Plan Complementario AstroHEP funded by the “European Union Next GenerationEU/PRTR” and the Government of the Autonomous Community of 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). Support for LAMBDA is provided by the NASA Office of Space Science. Some of the results in this paper have been derived using the Healpy and HEALPix packages (Górski et al. 2005; Zonca et al. 2019). The paper uses the python packages, Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020). DA thanks COSMOGLOBE project (https://www.cosmoglobe.uio.no/products) for making their spectral index map and masks publicly available.

References

Appendix A: Analysis using PR4 data at 30 GHz

In this section, we present a revised analysis incorporating the Planck PR4 data at 30 GHz, while all other frequency data remain the same as described in the main text. Figure A.2 shows the recovered moment maps. Compared to the results presented in Figure 9, the moment maps obtained using the PR4 data show no major differences relative to those derived using the PR3 data in the main analysis. Using the TT correlation between the recovered moment maps, we derive the β_s spectral index map and its corresponding standard deviation map, also shown in Figure A.2. The inverse-variance-weighted mean and dispersion of the derived spectral index map are –3.12 and 0.22, respectively. These values indicate no significant change relative to the results reported in the main text.

Fig. A.1. Synchrotron Stokes Qνs (first panel), and Uνs (second panel) maps at 22.8  GHz recovered from the cPILC component separation method using the QUIJOTE MFI 11 and 13  GHz combined map, WMAP K and Ka bands and Planck PR4 30 GHz. The Stokes parameters corresponding to the spectral index modulated synchrotron amplitudes, Qνsβs, Uνsβs obtained from component separation are shown in the third and fourth panels respectively. All maps are at FWHM = 2° at Nside = 64.

Fig. A.2. Synchrotron spectral index (left panel) and uncertainty (right panel) maps obtained from the combined QUIJOTE MFI, WMAP, and Planck PR4 data.

All Tables

All Figures

Fig. 1. Polarisation intensity, $P=\sqrt{Q^2+U^2}$, of the synchrotron zeroth moment (dash-dotted blue line) and the first moment (green) for two pivot values of βs = −3.0 (dash-dotted green line) and −2.97 (dashed green lines) at a pivot frequency, νs = 22.8 GHz. Emissions from dust, CMB, and AME (0.5% polarisation fraction) are presented by black, cyan, and magenta dash-dotted lines, respectively. The grey areas represent the frequency bands of QUIJOTE MFI bands, WMAP K and Ka-band, and Planck LFI 30 GHz.

In the text

	Fig. 2. Mask at Nside = 64 used in this analysis that retains 30% high S/N regions of the QUIJOTE MFI survey. The bright source positions are masked using P06 polarisation mask provided by WMAP.
In the text

	Fig. 3. Changes of the mean absolute deviation of the inferred spectral index with respect to the input synchrotron spectral index template, normalised with respect to the mean error of estimated βs with a different choice of the pivot value of the synchrotron spectral index.
In the text

Fig. 4. Comparison between input (upper panel) and recovered (lower panel) maps of Qνs(p), Uνs(p) and Qνs(p)βs(p), Uνs(p)βs(p) at the reference frequency of 23 GHz. All maps are shown at resolution of beam FWHM = 2°. The correlation between input and recovered maps is found to be 99% for Q and 95% for U. A similar correlation is found for Qνs(p)βs(p) and Uνs(p)βs(p) between input and recovered maps.

In the text

	Fig. 5. Correlation between recovered Stokes Q(p) (green) and U(p) (magenta) maps corresponding to zeroth and first order moments. The black dashed line is the linear fit with βs = −2.96 which corresponds to mean value of the input spectral index map.
In the text

	Fig. 6. Left panel: The input βs map used in the simulation. Middle panel: The inferred βs map from the estimated synchrotron moments. Right panel: The map of statistical and systematic uncertainty of the estimated βs which are added in quadrature.
In the text

Fig. 7. Upper panel: Histogram of the normalised deviation of the derived spectral index with respect to the true input spectral index $({\beta }_s^{\mathrm{Inferred}}-{\beta }_s^{\mathrm{Input}})/{\sigma }_{{\beta }_s}^{\mathrm{stat}+\mathrm{sys}}$ (orange). The black curve is a Gaussian distribution with zero mean and unit variance. Vertical lines are the corresponding mean of the distribution. Lower panel: The two-dimensional density plot of true input and inferred βs map. The correlation between the two maps is 81%. The black dashed line correspond to 100% correlation.

In the text

	Fig. 8. Histograms of the normalised deviation of the spectral index derived from simulations of three cases: without polarised AME (orange), with 0.5% (magenta), and 0.1% (green) polarised AME. The black curve is a unit variance Gaussian distribution with zero mean. Vertical lines are the corresponding mean of the distributions.
In the text

Fig. 9. Synchrotron Stokes Qνs (first panel) and Uνs (second panel) maps at 22.8 GHz recovered from the cPILC component separation method using the QUIJOTE MFI 11 and 13 GHz combined map, WMAP K and Ka bands and Planck LFI 30 GHz. The Stokes parameters corresponding to the spectral index modulated synchrotron amplitudes, Qνsβs, Uνsβs obtained from component separation are shown in the third and fourth panels respectively. All maps are at FWHM = 2° at Nside = 64.

In the text

	Fig. 10. Synchrotron spectral index (upper panel) and uncertainty maps (lower panel) obtained from combined the QUIJOTE MFI, WMAP, and Planck data.
In the text

	Fig. 11. Distribution of the synchrotron spectral index from ‘Model 4’ of Miville-Deschênes et al. (2008) (blue), de la Hoz et al. (2023) (green) and our estimation (orange). The vertical dashed lines correspond to the weighted mean values of the respective distributions.
In the text

Fig. 12. Synchrotron spectral index and corresponding uncertainty maps obtained using the QUIJOTE MFI 11 GHz map alone in replacement of QUIJOTE MFI combined maps of 11 and 13 GHz are shown in the first and second panels, respectively. The similar plots of synchrotron spectral index and corresponding uncertainty maps when we use QUIJOTE MFI 13 GHz data alone are shown in the third and fourth panels respectively.

In the text

	Fig. 13. Regions used for estimation of the spectral index defined in Fuskeland et al. (2014). The grey regions are masked since the regions are not covered by QUIJOTE MFI survey. Some of the bright sources are masked using P06 polarisation mask provided by WMAP.
In the text

	Fig. 14. Upper panel: The derived spectral index over the regions. Lower panel: The uncertainty map of the derived spectral index.
In the text

	Fig. 15. Synchrotron spectral index as a function of the region number obtained in this work (orange), COSMOGLOBE DR1 (black), and de la Hoz et al. (2023) (green). The data shown in green are the inverse-variance-weighted mean of spectral index and uncertainties computed from βs and σβs maps of de la Hoz et al. (2023) over regions.
In the text

Fig. A.1. Synchrotron Stokes Qνs (first panel), and Uνs (second panel) maps at 22.8  GHz recovered from the cPILC component separation method using the QUIJOTE MFI 11 and 13  GHz combined map, WMAP K and Ka bands and Planck PR4 30 GHz. The Stokes parameters corresponding to the spectral index modulated synchrotron amplitudes, Qνsβs, Uνsβs obtained from component separation are shown in the third and fourth panels respectively. All maps are at FWHM = 2° at Nside = 64.

In the text

	Fig. A.2. Synchrotron spectral index (left panel) and uncertainty (right panel) maps obtained from the combined QUIJOTE MFI, WMAP, and Planck PR4 data.
In the text