© The Authors 2024

1. Introduction

Understanding the polarization of the cosmic microwave background (CMB) will be a primary focus of observational cosmology in the coming decades. There has been phenomenal observational success over the past few decades, from the satellite-based COBE (Smoot et al. 1992; Mather et al. 1994; Hauser et al. 1998), Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2013), and Planck (Planck Collaboration I 2020) experiments, to the many sub-orbital experiments, including ACT (Madhavacheril et al. 2024), BICEP/Keck (Ade et al. 2021), CLASS (Eimer et al. 2024), QUIJOTE (Rubiño-Martín et al. 2023), SPIDER (SPIDER Collaboration 2022), SPT (Carlstrom et al. 2011), and many others. Future experiments, including those conducted with the LiteBIRD satellite (LiteBIRD Collaboration 2023), Simons Observatory (Ade et al. 2019), and CMB-S4 (Abazajian et al. 2019), will create the most sensitive maps of the polarized sky yet, resulting in stringent constraints on primordial gravitational waves.

The dominant astrophysical sources for polarized radiation in the microwave frequency range are synchrotron emission from relativistic electrons moving through the magnetic field of the Milky Way and thermal dust emission from vibrating dust grains. In order to map the valuable CMB fluctuations, both of these contaminants must be characterized and subtracted at high precision (Planck Collaboration X 2016). As shown by both Planck and BICEP2/Keck (Planck Collaboration II 2020; BICEP2/Keck Array and Planck Collaborations 2015), uncertainties on polarization-based cosmological constraints have been limited not only by instrumental sensitivity, but by incomplete knowledge of the sky itself. Uncertainty in the sky model has been mitigated by designing experiments with broad frequency coverage, such as WMAP (Bennett et al. 2013) and Planck (Planck Collaboration I 2020), or analyzing maps from different experiments jointly (e.g., Bennett et al. 1996, 2013, Planck Collaboration X 2016, and BICEP2/Keck Array and Planck Collaborations 2015, among many others). A major impediment to joint analysis is the difficulty of combining individually analyzed datasets with different survey strategies and incompletely characterized systematics. To maximize scientific throughput, one must either design an experiment that can characterize every relevant observable on its own or jointly analyze different datasets in the same joint framework. In practice, joint analysis has therefore until now usually been performed by combining datasets at the likelihood level.

The BEYONDPLANCK project performed a conceptually different form of joint analysis by directly combining the Planck Low Frequency Instrument (LFI) time-ordered data (TOD; BeyondPlanck Collaboration 2023) with external pixel-based data in the same analysis pipeline, including the Haslam 408 MHz map (Haslam et al. 1982), WMAP Ka–V bands, Planck 353 GHz in polarization, and Planck 857 GHz in intensity. Planck Collaboration II (2020) found that the LFI detector gain solution depended on the assumed polarization and intensity of the sky. To break this circular dependence, the BEYONDPLANCK framework solved for the intrinsic sky signal and instrumental parameters iteratively, providing an accurate model of the entire system with full error propagation. By leveraging the external WMAP, Haslam, and Planck High Frequency Instrument (HFI) data, BEYONDPLANCK was able to create Planck LFI maps of cosmological quality at all frequency channels, while simultaneously generating a robust model of the foreground sky.

The COSMOGLOBE project¹ has now generalized the BEYONDPLANCK analysis by also processing WMAP K–W-band data as TOD. This framework, fully described by Watts et al. (2023a,b), not only improves the quality of the WMAP maps themselves, but provides the most robust full-sky model of low frequency polarized emission available to date. In particular, poorly measured modes due to transmission imbalance (Jarosik et al. 2007) have been effectively marginalized over, in large part due to the use of a global sky model that was not available for the fiducial WMAP analyses. The removal of poorly measured modes alone will improve Galactic modeling for years to come, as the current state-of-the-art models use WMAP K-band (23 GHz) maps as polarized synchrotron templates for modeling synchrotron without any mitigation of the poorly measured modes (Delabrouille et al. 2013; Thorne et al. 2017; Zonca et al. 2021).

To use the synchrotron model, we need robust estimates of the spectral behavior. Fundamental properties of synchrotron emission, such as the spectral index as a function of the position on the sky, spatial de-correlation, and even the functional form of the spectral energy distribution (SED), are yet to be fully characterized observationally. Modern attempts have been frustrated by the lack of high signal-to-noise data; de Belsunce et al. (2022), for example, find spectral indices ranging from −5 to −1, which is discrepant with predictions from direct measurements of the cosmic ray energy spectrum of −2 to −3 (Rybicki & Lightman 1985; Orlando & Strong 2013; Neronov et al. 2017), largely due to fitting too many parameters to data with too low signal-to-noise ratios. Such analyses can be greatly improved by combining complementary datasets.

While the combination of as many datasets as possible would help constrain the large-scale properties of polarized synchrotron, long-standing discrepancies between existing datasets complicate this. As shown by, for example, Planck Collaboration X (2016) and Weiland et al. (2018), there are discrepancies in the polarization measurements of WMAP and Planck, partially due to known instrumental effects (Bennett et al. 2013; Planck Collaboration II 2020). BeyondPlanck Collaboration (2023) largely resolved these issues for LFI, while Watts et al. (2023b) removed the poorly measured modes in WMAP. In the first COSMOGLOBE data release (DR1), of which this paper is a part, Watts et al. (2023a) demonstrate that the polarized maps resulting from these two experiments can be free of systematics below the noise level if both experiments are jointly recalibrated. At the same time, the new maps must be validated by a variety of analysis methods. As noted in, for example, Weiland et al. (2022), simply using a different analysis pipeline can result in different estimates of the underlying frequency map. As such, we took care to appropriately marginalize over instrumental effects in our spectral analysis, and attempt to find the physical reason for the differences between different pipelines.

With these longstanding discrepancies addressed, the main goal of the current paper is to derive the highest signal-to-noise ratio model of polarized synchrotron emission from Planck and WMAP data published to date, including both amplitude and spectral index parameters. This work is organized as follows. In Sect. 2 we provide a brief overview of the data products and the COSMOGLOBE data processing. Section 3 discusses the polarization amplitude of the COSMOGLOBE synchrotron model, its overall uncertainty, and the power spectrum properties compared to external products. Section 4 estimates the synchrotron spectral index, employing both the T–T plot method and Gibbs sampling using Commander, which is followed by discussions and conclusions in Sect. 5.

2. Data products

The data products used in this paper are the WMAP and Planck LFI maps, with most of the statistical weight coming from 23–33 GHz, or the WMAP K and Ka bands and LFI’s 30 GHz band. These frequencies are low enough that we can treat them as synchrotron tracers and hence ignore thermal dust and CMB emission, but not so low that we need to take into account effects like Faraday rotation, as in S-PASS data (Krachmalnicoff et al. 2018; Fuskeland et al. 2021). We describe the legacy data products and the data products from COSMOGLOBE DR1 in the following subsections.

2.1. WMAP and Planck legacy products

WMAP was a NASA-funded satellite mission that observed from August 2001 to August 2010, designed to characterize the microwave sky well enough to measure the primary CMB anisotropies across the full sky down to a resolution of 13′ full width at half maximum (FWHM). Using a differential scanning strategy inspired by COBE Differential Microwave Radiometer (DMR), WMAP produced maps of the sky at 23 (K), 33 (Ka), 41 (Q), 61 (V), and 94 GHz (W) in both polarization and total intensity (Bennett et al. 2013), with angular resolutions of 53′ at 23 GHz to 13′ at 94 GHz. The maps are available on the LAMBDA website².

The Planck LFI produced 30, 44, and 70 GHz maps in both intensity and polarization, while the HFI produced 100, 143, 217, and 353 GHz maps in polarization and intensity, and 545 and 857 GHz maps in intensity alone. The LFI data have less white noise and higher angular resolutions (of 30′, 20′, and 13′ FWHM for 30, 44, and 70 GHz) than WMAP. In contrast to WMAP, the LFI measurements used a single horn, and the Planck scanning strategy followed rings closely aligned with ecliptic meridians. The Planck legacy datasets, PR3 (Planck Collaboration I 2020) and PR4 (Planck Collaboration Int LVII 2020), are both publicly available on the Planck Legacy Archive (PLA)³.

2.2. COSMOGLOBE products

A main goal of COSMOGLOBE is to perform joint end-to-end analyses on multiple datasets, preferably beginning from raw TOD. An important advantage of such end-to-end processing is that it offers a robust path to breaking internal degeneracies within and between different datasets, in general reducing the magnitude of systematic effects. The analysis described by BeyondPlanck Collaboration (2023) and Watts et al. (2023a) were performed on raw TOD, producing cosmological parameters using a Bayesian Gibbs sampler called Commander3 (Galloway et al. 2023).

In this framework, we produced a full sky model and set of instrumental parameters for each Gibbs sample, and the set of all such samples allows for a thorough characterization of the dependence of low-level instrumental parameters on the sky model. The sky model includes all relevant components in the WMAP and LFI frequency range, specifically the CMB, synchrotron, thermal dust, free-free emission, anomalous microwave emission, and radio point sources, the first three of which we modeled as polarized. The full products from this analysis and individual maps are available on the COSMOGLOBE website⁴.

As described by Watts et al. (2023a), COSMOGLOBE DR1 includes 500 end-to-end Gibbs samples; we performed our basic analysis on each of these samples individually and then formed posterior summary statistics from the full ensemble. This allowed us to fully marginalize over the low-level systematic parameters, quantifying the extent to which instrumental processing propagates to the synchrotron spectral index determination.

3. Polarized synchrotron amplitude

In this section, we give an overview of the polarized synchrotron amplitude properties. In Sect. 3.1 we focus on the map-based properties, and compare with independently processed results in Sect. 3.2. In Sect. 3.3 we evaluate the power spectrum of the polarized synchrotron map and compute the B-to-E ratio.

3.1. Polarization amplitude maps

We start by comparing the polarized spectral amplitude (as defined by $P=\sqrt{Q^2+U^2}$) from the COSMOGLOBE synchrotron map with the official WMAP K-band and Planck 30 GHz maps in Fig. 1. The noise level of the COSMOGLOBE synchrotron map is already at a visual level clearly lower than those of the K-band and 30 GHz maps, even taking into account the scaling of the K-band map to 30 GHz. This is particularly visible in the leftmost column of Fig. 1, an inset of a high Galactic latitude, low signal region.

Fig. 1. Polarized intensity and white noise levels of the WMAP K band (top), Planck 30 GHz band (middle), and synchrotron amplitude from the COSMOGLOBE Gibbs chain (bottom), all evaluated at 30 GHz with a resolution of 72′. All maps have been scaled to a common frequency of 30 GHz. The leftmost column is a 10° wide square centered on the low signal-to-noise region highlighted by a red square in the middle column’s panels. The middle column shows the polarized total amplitude, P. The rightmost column shows the rms noise for the K-band and 30 GHz frequency maps and the COSMOGLOBE DR1 synchrotron amplitude map.

To quantify the noise improvement in the synchrotron maps over pure templates based on either K-band or 30 GHz maps, we compare the posterior standard deviation of the COSMOGLOBE synchrotron map with the frequency maps’ rms values per pixel. To do so, it is important to take into account the fact that the COSMOGLOBE synchrotron map has effectively been convolved with a 30′ FWHM Gaussian smoothing prior, which suppresses small-scale fluctuations, resulting in a 72′ resolution. To compare the K-band and 30 GHz noise levels with this product, we therefore simulate noise realizations, smooth with a corresponding beam, and scale the maps to the default reference frequency of 30 GHz assuming a polarized synchrotron spectral index β_s = −3.1, consistent with the post-processed synchrotron map. We display these smoothed rms maps in the right column of Fig. 1.

At a FWHM of 2°, we find that the mean rms’s for K-band and 30 GHz maps are 4.8 μK and 4.7 μK, respectively, compared to the mean value of 3.4 μK for the COSMOGLOBE DR1 synchrotron map⁵; this is consistent with adding the scaled maps in quadrature. While this may seem obvious on its face, the combination of WMAP K-band and Planck 30 GHz maps is not straightforward due to instrumental effects that remain in the maps, notably poorly measured modes in WMAP (Bennett et al. 2013; Weiland et al. 2018) and gain uncertainty in Planck (Planck Collaboration II 2020). Although the calculated white noise level from adding COSMOGLOBEK-band and 30 GHz maps are not significantly lower than with the PR4 and WMAP9 products, this combination of data yields a synchrotron map with reduced instrumental effects. These effects remain in the official WMAP9 (Bennett et al. 2013) and Planck PR3/PR4 (Planck Collaboration II 2020; Planck Collaboration Int LVII 2020) maps, requiring these effects to be carefully modeled and projected out before being used for Galactic science and cosmological analyses. However, the end-to-end BEYONDPLANCK (BeyondPlanck Collaboration 2023) products effectively removed the gain uncertainty modes from the Planck LFI maps, and the COSMOGLOBE DR1 results (Watts et al. 2023a) are free from the poorly measured modes. As shown in Fig. 46 of Watts et al. (2023a), the Planck 30 GHz and the WMAP K-band maps are consistent with each other at the 10 μK level. This is slightly larger than the expected white noise level of 3.4 μK, and can be explained by un-modeled deviations from β_s = −3.1 or systematics. In general, the COSMOGLOBE frequency maps have similar white noise levels to the official products, but with lower systematics and better recovery of large angular scales. Our polarized synchrotron map, derived from consistent datasets, has a white noise level 29% lower than both the K-band map and the 30 GHz map.

In Fig. 2 we plot the ratio between the COSMOGLOBE and BEYONDPLANCK (BeyondPlanck Collaboration 2023) synchrotron rms maps. Both these products are created using the Commander3 pipeline and similar data selection, but COSMOGLOBE benefits from additionally using the K-band data. The COSMOGLOBE synchrotron map has a mean posterior rms of 3.4 μK, nearly a factor of two improvement over the BEYONDPLANCK map. As expected, the regions with the lowest ratios correspond to the deep WMAP observations. The synchrotron in the Galactic plane has a higher signal-to-noise due to the high signal-to-noise K-band data, which stabilizes the relative calibration in the 30 GHz data. These two effects combine to give a smaller posterior standard deviation in the Galactic plane, as there is less variation in the component separation step between Gibbs samples. Regions with nearly identical rms include regions with the highest Planck depth, such as the ecliptic poles, and regions less deeply observed by WMAP, corresponding to planet crossings and artifacts from the processing mask. Notably, stripes corresponding the Planck scan strategy also show improvement with respect to BEYONDPLANCK. This is due to the interaction between the sky model and LFI’s instrumental parameters – with the high signal-to-noise K-band data, the sky model becomes more stable, and LFI’s relative gain solution becomes better determined.

Fig. 2. Ratio between the COSMOGLOBE and BEYONDPLANCK polarized synchrotron amplitude standard deviations per HEALPixNside = 1024 pixel. The posterior standard deviation is computed over all Gibbs samples for both the COSMOGLOBE and BEYONDPLANCK chains.

3.2. Comparison with independent datasets

Next, we compared the COSMOGLOBE polarized synchrotron map with the frequency maps produced in the main DR1 chain, paying special attention to maps with the highest polarized synchrotron signal-to-noise ratio; K, Ka, 30 GHz, and 44 GHz. Using the polarized synchrotron model generated in the main COSMOGLOBE DR1 chain, we evaluated the synchrotron emission at each frequency and subtracted this from frequency maps produced by various different processing pipelines. We evaluated the COSMOGLOBE sky model employing the cosmoglobe Python package, which uses the full bandpass information of each instrument⁶. Because the COSMOGLOBE sky model is fit using WMAP data with mitigated poorly measured modes, the model will necessarily provide a better fit to the COSMOGLOBEK-band maps than WMAP9 K-band maps. Had the model been fit using the WMAP9 maps, the WMAP9 K-band residuals would be more consistent with white noise while causing larger residuals in the WMAP9 and COSMOGLOBEKa-band maps.

The first and third rows of Fig. 3 show the WMAP9 and PR3 residuals with respect to the COSMOGLOBE sky model. These residuals match previously documented observational effects in both experiments. In particular, the WMAP maps show the poorly measured modes due to transmission imbalance in the differential horns (Jarosik et al. 2007; Bennett et al. 2013), while the PR3 differences are mostly due to relative gain errors (Planck Collaboration II 2020; Planck Collaboration Int LVII 2020). In contrast, as seen in the second and fourth rows of Fig. 3, the synchrotron model matches the COSMOGLOBE DR1 frequency maps within 5 μK across the sky, with few observational artifacts.

Fig. 3. Frequency map minus sky model with respect to the COSMOGLOBE DR1 sky model evaluated at 5°. COSMOGLOBE maps are labeled CG, Planck 2018 maps are labeled PR3, and the legacy WMAP9 maps are labeled WMAP. Each of the WMAP9, Planck PR3, and COSMOGLOBE maps have had the COSMOGLOBE DR1 sky model evaluated at their respective bandpass and the resolution subtracted.

Most of the COSMOGLOBE residuals are associated with the Galactic plane and diffuse structures uncorrelated with the WMAP and Planck observation strategies. The residuals in the WMAP9 K and Ka maps include imprints of the Galaxy due to relative gain differences and poorly measured modes. Neither of these imprints are present in the COSMOGLOBEK and Ka map residuals, although the removal of the poorly measured modes has made the existence of a positive excess in the Ka Stokes Q map more apparent. This residual strongly resembles the bandpass correction in the WMAP K band, suggesting that an improved modeling of bandpass corrections can further improve the Ka-band residuals. The LFI residuals, while much improved, still show trace residuals, especially near the Galactic center, that are somewhat correlated with the gain correction templates, but not at a level that high Galactic latitude features can be identified.

3.3. Power spectra

Next, we considered the angular power spectrum of polarized synchrotron emission evaluated at 30 GHz. To estimate the power spectra without the positive-definite bias induced by noise fluctuations, we performed a Commander3 run using half-mission splits with odd-numbered scans and even-numbered scans being analyzed in runs labeled HM1 and HM2, similar to the DetSet approach used by Planck Collaboration Int XXX (2016). Each of these chains were performed using the same data as in the main COSMOGLOBE chain, with 200 samples each⁷. The highest quality of similar half-mission splits that are publicly available are from the Planck PR3 analysis, as discussed by Planck Collaboration IV (2020). We therefore computed power spectra from the PR3 results and compared them directly with the COSMOGLOBE HM splits.

We performed power spectrum estimation using a so-called pseudo-spectrum code NaMaster (LSST Dark Energy Science Collaboration 2019), coupled with the Planck 2018 common polarization mask (f_sky = 0.78) and 1° apodization. To quantify the uncertainty, we took the cross spectrum for each pair of Gibbs samples from HM1 and HM2, respectively, and report the 68% confidence intervals on this posterior. The power spectra are displayed in Fig. 4, and the standard deviation is computed using the within-bin variance of each bin, with the posterior standard deviation of the Gibbs chain added in quadrature for the COSMOGLOBE spectra. Other than the very lowest and very highest multipole bins, there is good per-multipole agreement between both the PR3 and COSMOGLOBE DR1 spectra. For comparison, we plotted the E-mode power spectrum predicted by the Planck 2018 cosmological parameters.

Fig. 4. Half-mission cross spectra for polarized synchrotron emission as evaluated for Planck PR3 (blue) and COSMOGLOBE DR1 (black). Filled circles correspond to E modes, while empty circles correspond to B modes. The thick solid and dashed lines are the best-fit E- and B-mode power law fits to synchrotron, and the thin black line is the Λ cold dark matter prediction for E modes. The lowest bin is not included in either fit.

Following Planck Collaboration IV (2020), we performed power law fits to the power spectra of the form ${\mathcal{D}}_{\ell }^{\mathcal{EE}/\mathcal{BB}}={\mathcal{A}}_{\mathrm{s}}^{\mathcal{EE}/\mathcal{BB}}{(\ell /\mathcal{80})}^{\alpha }$, using multipoles ℓ ∈ [4, 140]. The 68% confidence intervals for each quantity, including the $A_{\mathrm{s}}^{\mathit{BB}}/A_{\mathrm{s}}^{\mathit{EE}}$ ratio, are reported in Table 1. The primary differences between the fits to the two datasets are ∼2 σ discrepancies in the $A_{\mathrm{s}}^{\mathit{BB}}$ and ${\alpha }_{\mathrm{s}}^{\mathit{EE}}$ fits, while all others are consistent within ≃0.5 σ. The primary drivers of these differences are lower ${\mathcal{D}}_{\mathcal{l}}^{BB}$ and higher ${\mathcal{D}}_{\mathcal{l}}^{EE}$ in the highest bins, indicating potential bias in the PR3 spectra.

A question of special interest is the ratio of synchrotron B-mode to E-mode power, which has been consistently noted to be less than one (Page et al. 2007; Planck Collaboration X 2016; Planck Collaboration IV 2020; Krachmalnicoff et al. 2018; Martire et al. 2022; Rubiño-Martín et al. 2023; Eimer et al. 2024). There is also a potential dependence of this ratio on the mask choice, as suggested by results in Krachmalnicoff et al. (2018), Martire et al. (2022), Rubiño-Martín et al. (2023), and Eimer et al. (2024). The physical mechanism responsible for the ratio being less than one and for dependence on sky mask has been discussed in the context of Galactic magnetic fields and polarized thermal dust (e.g., Kandel et al. 2017, Hensley et al. 2022, Vacher et al. 2023), and similar mechanisms are likely to be relevant for synchrotron polarization. In Table 1, we report $A_{\mathrm{s}}^{\mathit{BB}}/A_{\mathrm{s}}^{\mathit{EE}}$ of 0.39 for the COSMOGLOBE splits and 0.47 for PR3. In this respect, we note that our value of $A_{\mathrm{s}}^{\mathit{BB}}/A_{\mathrm{s}}^{\mathit{EE}}$ for Planck 2018 is higher than that reported by Planck Collaboration IV (2020) of 0.34, despite using a nearly identical methodology. The origin of this discrepancy is not yet understood, but it must necessarily be associated with details in the power spectrum estimation algorithm. For a quantitative sense of the stability of these results, changing the minimum multipole from ℓ = 4 to ℓ = 2 alters the $A_{\mathrm{s}}^{\mathit{BB}}/A_{\mathrm{s}}^{\mathit{EE}}$ ratios from 0.39 to 0.40 and 0.46 to 0.47 for COSMOGLOBE and PR3, respectively, a 0.5σ shift. As the focus of this analysis is to check consistency with the values with the PR3 results, we leave the exploration of the power spectrum ratio’s dependence on the mask and multipole range to further work.

4. Polarized synchrotron spectral indices

Next, we turned our attention to the SED of polarized synchrotron emission, which we in this section will assume scales as a perfect power law in frequency, ν^β_s. The determination of β_s across the sky has already been studied in detail using several different data combinations and methods. Although the small-scale details vary, nearly every analysis has found β_s ≃ −2.8 in the Galactic plane and β_s ≃ −3.3 in high Galactic latitudes (Fuskeland et al. 2014, 2021; Krachmalnicoff et al. 2018; Weiland et al. 2022), with the exception of QUIJOTE (Rubiño-Martín et al. 2023; de la Hoz et al. 2023), who find a slightly steeper spectral index along the Galactic plane. Both Fuskeland et al. (2014) and Weiland et al. (2022) report oscillations with Galactic longitude close to the Galactic plane, but high-latitude regions variations are more difficult to determine, and tend to depend on the specific dataset chosen and the analysis method chosen.

To probe the stability of our results with respect to algorithmic details, we performed two different analyses, namely a so-called T–T plot analysis (see Sect. 5) and then a Gibbs sampling analysis using Commander1 (Eriksen et al. 2008; see Sect. 4.2)⁸. Section 5 focuses on pairs of channels to better isolate potential un-modeled systematic effects, while Sect. 4.2 uses all WMAP and LFI channels plus Planck 353 GHz to maximize the joint statistical weight of all available data.

Fig. 5. Spectral index regions as defined by Fuskeland et al. (2014). The most prominent point sources are masked out and shown in gray circular areas.

While the COSMOGLOBE pipeline does produce samples of β_s as part of the main Gibbs chain, the full marginal signal-to-noise ratio per pixel with respect to β_s is generally very low due to degeneracies with respect to gain and bandpass uncertainties (Gjerløw et al. 2023; Svalheim et al. 2023b), and this makes it susceptible to systematic uncertainties. In preliminary runs where β_s was sampled with weak priors, we found that errors in low-level parameters such as gain and bandpass shift could be mimicked by β_s variation in high Galactic latitude low signal-to-noise regions. In high signal-to-noise regions, the spectral index was well determined, albeit with a small systematic uncertainty, which is discussed further in Sect. 4.2. To avoid long correlation lengths with unphysical parameters, the default COSMOGLOBE DR1 processing adopts an informative β_s ∼ 𝒩(−3.1.5,0.05) prior sampling (Watts et al. 2023a). This choice has direct implications for the estimates of β_s presented in this section, which are derived from the COSMOGLOBE frequency maps. To at least partially quantify this effect, we considered various different prior choices (see Sect. 4.2), and we accounted for shifts resulting from different prior means into our final uncertainties.

4.1. T–T plot analysis

The T–T plot analysis can be used for pairs of datasets in a model-independent way, allowing for us to probe the difference between different datasets without making strong assumptions on the underlying physical model. This method can also easily be adapted to probe the dependence on polarization angle, and it can be used to study the orientation dependence of β_s with respect to un-modeled instrumental effects, as discussed by Wehus et al. (2013).

4.1.1. Method

Following Fuskeland et al. (2014, 2021), we applied linear regression via the T–T plot method, in which spectral indices can be estimated over extended regions with approximately constant spectral indices. Here we used the regions labeled in Fig. 5⁹. In this approach, our data model reads

$$\begin{array}{c}\hfill {\mathit{m}}_{\nu }={\mathit{m}}_{{\nu }_0}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_{\mathrm{s}}}+c_{\nu }+{\mathit{n}}_{\nu },\end{array}$$(1)

where m_ν is a spatially varying amplitude map at frequency ν; ν₀ is a reference frequency; β_s is a power law index fit jointly for Q and U; c_ν is the spatially constant offset per band; and n_ν represents Gaussian noise. In the case of noiseless data with no offset, the spectral index may be estimated from any two frequency maps, ν₁ and ν₂, using a simple ratio,

$$\begin{array}{c}\hfill \frac{m_{{\nu }_1,p}}{m_{{\nu }_2,p}}={\left(\frac{{\nu }_1}{{\nu }_2}\right)}^{{\beta }_{\mathrm{s},p}}\Rightarrow {\beta }_{\mathrm{s},p}=\frac{ln(m_{{\nu }_1,p}/m_{{\nu }_2,p})}{ln({\nu }_1/{\nu }_2)}.\end{array}$$(2)

In the case where one map is much noisier than the other, the standard T–T plot method involves performing a linear regression m_ν1 = am_ν2+ b, and associates β_s with lna/ln(ν₁/ν₂). More care must be taken when the noise amplitudes in both maps are comparable to each other, so we adopted the effective variance method of Orear (1982) as implemented by Fuskeland et al. (2014).

To obtain robust results that are independent of the orientation of the coordinate system, we marginalized the spectral index over polarization angle as explained in Fuskeland et al. (2014, 2021). Here we produced separate datasets considering angles, α, between 0° and 85° in steps of 5°, and conducted the T–T plot analysis on all of these datasets. In some figures we report all the resulting spectral indices, and as a final result per region we took the inverse variance weighted mean.

As described by Fuskeland et al. (2021), two sources of systematic uncertainties are added in quadrature to the statistical uncertainty. The first is a bootstrap uncertainty where we randomly drew 10 000 new datasets from the original dataset and adopted the standard deviation of the resulting spectral index as a systematic uncertainty. The second effect is a systematic uncertainty that takes into account the variation in β_s over rotation angle, [max(β_s, α)−min(β_s, α)]/2. For the COSMOGLOBE analyses in this paper we have a whole suite of maps represented by the individual samples instead of just one mean map. Here the standard deviation of the spectral indices of the samples is also added in quadrature to represent an additional systematic uncertainty. This enables us to have a better propagation of uncertainties from the maps to the final spectral indices. The total uncertainty of the spectral index in a region is calculated as the minimum of the uncertainties in each rotation angle, and for the final spectral index as the minimum of all regions, noting that the uncertainty should never increase when we add more data.

For the T–T plot analysis, we focused exclusively on bands between 23 and 33 GHz. Following Fuskeland et al. (2014), we used the WMAP K- and Ka-band Stokes Q and U parameter maps at 23 GHz and 33 GHz. The respective effective frequencies used are 22.45 GHz and 32.64 GHz for a spectral index of −2.7 (Page et al. 2003). We note that the values for the effective frequencies when the spectral index varies between −2.8 and −3.2 only yields a sub-percent effect on β_s. The maps originally at a HEALPix¹⁰ pixelization of N_side = 512 are downgraded to N_side = 64 and smoothed to a common resolution of 1° FWHM. The Planck data products used are the 30 GHz Stokes Q and U maps, with an effective frequency of 28.4 GHz. This value is when bandpass shifts are applied, the nominal value without bandpass correction is 28.1 GHz. The difference between these two effective frequencies leads to a difference of 4.5% in the spectral index. Formally, the uncertainty of 0.3 GHz in the LFI 30 bandpass reported in Svalheim et al. (2023b) and Watts et al. (2023a) should be propagated, but the effective error of order 0.01 in β_s is negligible compared to the instrumental and Monte Carlo uncertainties, and is therefore ignored in this analysis. Both the BEYONDPLANCK and COSMOGLOBE products are natively at N_side = 512, while for the Planck PR4, it is 1024.

4.1.2. Results

We show the full set of T–T scatter plots between WMAP K-band and LFI 30 GHz maps in Fig. 6. Overall, the two sets of point clouds appear to be in reasonable good agreement at this visual level. Perhaps the most notable features are the fact that while the best-fit lines are generally overlapping in most low Galactic latitude regions, the COSMOGLOBE slopes are generally a bit shallower than the official maps at high Galactic regions. The most discrepant case is region 22 along the Galactic plane, in which the best-fit line for the PR3 maps is poorly aligned with the diffuse cloud, indicating that a substantial fraction of the pixels have anomalous values. This is typically indicative of strong systematic effects, for instance incomplete temperature-to-polarization leakage.

Fig. 6. T–T plots for Stokes Q and U maps of the COSMOGLOBE DR1 K- band versus the COSMOGLOBE DR1 30 GHz (black) and the WMAP9 K- band versus Planck PR3 30 GHz pixels (red) for all regions. The horizontal (solid and dotted) lines indicate the corresponding inverse variance weighted values of the spectral index, averaged over rotation angle, and in the COSMOGLOBE case also weighted over samples. The corresponding βs values can be found in Fig. 8.

The information in Fig. 6 is significantly compressed in Fig. 7, which shows the best-fit spectral index as a function of polarization angle (for full discussion of this type of plots, see, e.g., Wehus et al. 2013 and Fuskeland et al. 2021). Variations in these curves could in principle be indicative of true β_s dependence on polarization angle within a region, but very large variations are typically associated with uncorrected instrumental systematic errors. Comparing the COSMOGLOBE (black curves) and official maps (red curves), we generally see that the former show less internal variations than the latter, suggesting greater consistency between the Stokes parameters. For region 22 we see that the official maps fall outside the plotted range. As noted above, we inverse-variance weighed these results when reporting final co-added spectral indices in the following, and we included the Stokes’ variation in the final uncertainties.

Fig. 7. Synchrotron spectral index as a function of rotation angle, computed using a T–T plot between the COSMOGLOBE DR1 K-band and the COSMOGLOBE DR1 30 GHz pixels (black) compared to the spectral index using the WMAP9 K-band and Planck PR3 30 GHz (red) pixels for all regions. The horizontal (solid and dotted) lines indicates the corresponding inverse variance weighted values of the spectral index, averaged over rotation angle, and in the COSMOGLOBE case also samples.

Figure 8 shows final spectral index estimates of this type for four different generations of WMAP K-band and Planck 30 GHz data in the form of regionalized sky maps. The top three maps show results derived from the official WMAP K-band map versus Planck PR3, PR4, and BEYONDPLANCK 30 GHz. The bottom map shows the results using the improved maps from the current COSMOGLOBE analysis. While all cases are inverse variance weighted over rotation angle, the COSMOGLOBE case is additionally averaged over posterior samples, accounting for low-level instrumental uncertainties. Some of the brightest point sources have been masked out using gray.

Fig. 8. Spatial variation in the synchrotron spectral index, computed using T–T plots between the (from top to bottom) WMAP9 K-band and Planck PR3 30 GHz, WMAP9 K-band and Planck PR4 30 GHz, WMAP9 K-band and BEYONDPLANCK 30 GHz, and COSMOGLOBEK-band and COSMOGLOBE 30 GHz maps. The spectral index is inverse variance weighted over rotation angle, and in the COSMOGLOBE case also samples.

We see in this figure an overall narrowing in the range of spectral spectral indices, from regions with extreme values in the top map represented by dark red and blue values, with colors gradually fading in the lower figures. The range in the plot is quite wide for spectral indices, and goes from −5 to −1. This improvement in the value of the spectral index is especially prominent in the high Galactic regions, as well as region number 22 along the Galactic plane. Using the COSMOGLOBE DR1 data, we find a full-sky inverse variance weighted mean of the spectral index of β_s = −3.07 ± 0.07, while the mean of the high latitude regions 1–14 is β_s = −3.31 ± 0.07.

The same results are also summarized in Fig. 9, with the spectral indices plotted as a function of region number. The different datasets are shown in different colors, and with the same order as the four maps in Fig. 8. Here we have also included one sigma uncertainties as error bars as described in Sect. 4.1.1. The magnitudes of the error bars are set by the reported white noise level in each dataset except for COSMOGLOBE, which also takes into account systematic uncertainties through sample averaging. We also plot the estimated spectral index range for high Galactic regions, in our case represented by regions 1–12, as derived by the Planck PR3 LFI likelihood analysis (see Fig. 22 of Planck Collaboration V 2020). With each iteration of the data processing, the region estimates draw closer to this result.

Fig. 9. Synchrotron spectral index as a function of region number, computed using T–T plots between the WMAP9 K-band and Planck PR3 30 GHz (red), WMAP9 K-band and Planck PR4 30 GHz (blue), WMAP9 K-band and BEYONDPLANCK 30 GHz (purple), and COSMOGLOBEK-band and COSMOGLOBE 30 GHz (black) maps. The spectral index is inverse variance weighted over rotation angles, and samples. The horizontal line in the high latitude regions corresponds to the estimated spectral index values from the Planck PR3 likelihood analysis (Planck Collaboration V 2020).

When interpreting the variations between different data generations seen in Fig. 9, it is useful to overlay the region boundaries on the raw frequency map differences between the COSMOGLOBE and official WMAP9 maps (Watts et al. 2023a); these are shown in Fig. 10 for the K-band and Planck 30 GHz maps. As discussed by Watts et al. (2023a), at high Galactic latitudes these large-scale patterns are largely dominated by different gain (for Planck) and transmission imbalance (for WMAP) estimates. The most extreme outlier in Fig. 9 is region 4, close to the North Galactic Pole, for which the official maps yield a mean spectral index of about β_s ≈ −5. Inspecting Fig. 10, we see that that region has a strong eastward gradient from positive to negative in the Planck PR3 30 GHz difference, while the opposite is true for the WMAP K band. The PR4 30 GHz map also contains such a gradient within this region with a smaller amplitude. Similar considerations hold for many other regions as well.

Fig. 10. Differences between the COSMOGLOBE and official WMAP 23 GHz (top row) and between COSMOGLOBE and Planck PR3 30 GHz (bottom row) maps, with spectral index regions overlaid. Columns show Stokes Q and U parameters.

Next, we considered similar T–T plots for the K-band and Ka-band maps. For the WMAP data, this analysis was already performed by Fuskeland et al. (2014), so here we only present results using the new COSMOGLOBE data. First, as a consistency check, in Fig. 11 we plot the spectral indices for all 24 regions in the form of COSMOGLOBEK/30 GHz versus COSMOGLOBEK/Ka. We see that for the regions in the north and south Galactic spurs, and along the Galactic plane (regions 13–24), there is a good agreement between the spectral indices obtained by the two pairs of datasets. Regions 1–12, all of which are high-latitude regions with low signal-to-noise, are consistent with the two data combinations. The most prominent outliers, regions 6 and 12, are the regions with the lowest signal-to-noise ratio, and are thus more prone to noise fluctuations.

Fig. 11. Synchrotron spectral index computed using the T–T plot method with the COSMOGLOBE DR1 K-band (23 GHz) and 30 GHz data versus COSMOGLOBE DR1 K-band and Ka-band (33 GHz) data.

In Fig. 12 the value of the spectral indices is shown for the 24 regions. We see here that the colors are even fainter than the COSMOGLOBEK-band versus 30 GHz values (the bottom panel in Fig. 8), meaning there are fewer outliers with respect to a standard value in the range around −3. This is especially visible in the regions with lowest signal to noise, like regions 6, 8 and 12, at high latitude. Taking the inverse variance weighted mean of the spectral indices of all 24 regions, we get β_s = −2.95 ± 0.07, while for the high latitude regions 1–14 we get −3.20 ± 0.10.

Fig. 12. Spatial variation in the synchrotron spectral index, computed using the T–T plot method between the COSMOGLOBE DR1 K- and Ka-band maps. The spectral index is inverse variance weighted over the rotation angle and samples.

4.2. Parametric component estimation

We used Commander1¹¹ to perform pixel-based component separation on maps smoothed to a common resolution of 5° and at N_side = 64. For this analysis, we generated 100 pixel-based Commander1 Gibbs samples for each of the 500 COSMOGLOBE DR1 Gibbs samples produced by Commander3. This allowed us to decompose the pure statistical error assuming white noise alone for each of the 100 Commander1 samples, while changes between each main DR1 sample show the effects of low-level instrumental processing. In these analyses, we used the same polarized bands as in the main COSMOGLOBE DR1 analysis, namely WMAP, Planck LFI, and Planck 353 GHz from PR4. We used a prior β_s ∼ 𝒩(−3.1,0.1) for the synchrotron spectral index, which is looser than the algorithmic prior used for the main DR1 processing, but still tight enough to ensure physical values at high Galactic latitudes. For thermal dust emission, we adopted the dust temperature map, T_d, from Planck Collaboration X (2016), coupled with a constant spectral index of β_d = 1.55; minor variations in the thermal dust model has only a small impact on the synchrotron spectral index. The dust amplitude was fit as a free parameter conditioned on these values.

Overall, we considered the same data model as in COSMOGLOBE DR1 but allowed for a spatially varying β_s in two different forms. First, we allowed β_s to vary over the same regions as in the T–T analysis. We find that for such large regions, there is only a small effect due to the prior. For instance, when considering very different priors of 𝒩(−3.5, 0.1) and 𝒩(−2.7, 0.1), we find a maximum difference of Δβ_s = 0.3 in the lowest signal-to-noise region; generally it is much smaller than this. We find that a flat prior of U(−5, −1) provided reasonable results, but used the priors with means of −3.5 and −2.7 to assess the potential impact of additional external data. We then added the resulting differences in quadrature to the uncertainty due to statistical and systematic errors to account for prior ignorance. We did this by running cases with priors centered at −2.7, −3.1, and −3.5, and used the differences between the −2.7 and −3.5 cases as an estimate of prior-based uncertainty. All results are quoted with the prior centered at −3.1.

The results from this analysis are summarized in Fig. 13 for the nominal priors, and here we also compare our results with corresponding results derived by QUIJOTE (de la Hoz et al. 2023) and CLASS (Eimer et al. 2024). We chose these datasets because they are to date the most sensitive independent datasets probing a frequency range that is comparable to WMAP and LFI; for comparison S-PASS is much more strongly affected by Faraday rotation (Krachmalnicoff et al. 2018; Fuskeland et al. 2021). The β_s maps delivered by the two projects are pixelized at N_side = 64 and 32, respectively, with associated uncertainty maps taking into account the expected instrumental noise levels. We took an inverse-weighted average of the respective maps and report the weighted standard error within each region, displayed in Fig. 13. We only display regions for which 75% of the pixels lie within the region in question to better report the agreement between different experiments.

Fig. 13. Spectral index estimates from COSMOGLOBE using T–T plots with K and Ka maps and Gibbs sampling using all frequencies, and from QUIJOTE and CLASS using T–T plots. Only regions with ≥75% of the pixels within a given region are displayed. Red points and black points are both calculated using COSMOGLOBE frequency maps, but are obtained with different techniques.

In general, the uncertainties in our parametric analysis are smaller than each of the other published results, each for slightly different reasons. First, the T–T analyses will inherently have less constraining power than a full likelihood analysis, as this approach only uses two frequency channels via a linear regression, so the uncertainty is determined by the noise level in each frequency channel and the inherent variation within a given sky region. Beyond that, the T–T plot in this paper marginalizes over dependence on the polarization angle α, and by design accentuates systematic effects, such as beam ellipticity. This is best understood when inspecting the effect of beam ellipticity on a point source, as in, for example, Wehus et al. (2013). By convolving a point source with an elliptical beam, spurious polarization is created, leading to an artificial variation as a function of the polarization angle, increasing the error bar of β_s. Finally, the QUIJOTE analysis is most similar in data choice (MFI 11/13 GHz, WMAP9 K/Ka, Planck PR4) and methodology (B-SeCRET; de la Hoz et al. 2022), but still yields higher uncertainty than the Commander1 spectral index region analysis. This is most likely due to different spatial resolution and modeling choices. In particular, de la Hoz et al. (2023) sampled for β_d and T_d with priors 𝒩(1.55, 0.1) and 𝒩(21, 3) while using a relatively wide prior on β_s of 𝒩(−3.1, 0.3). Another modeling difference is that the Commander1 analysis presented here is performed at 5° resolution versus the 2° resolution per-pixel analysis by de la Hoz et al. (2023), although this is unlikely to affect the final results.

At high Galactic latitudes (regions 1–12), there is good agreement between each of the treatments, and all values are consistent with a single constant value. In the north and south polar spur (regions 13 and 14) there is excellent agreement between all of the pipelines, while along the Galactic plane (regions 15–24), there are mild differences between the values. This is almost entirely due to the differing coverage between CLASS, QUIJOTE, and the COSMOGLOBE maps. When restricting the analysis to use pixel regions that are identical within all experiments, we find that all spectral index estimates are consistent within 1σ.

To more finely probe the spatial variation in β_s, we performed a second Commander1 analysis with identical data and model choices, except β_s was now allowed to vary within each N_side = 16 HEALPix pixel. At higher resolutions than this, the spectral index is prior dominated at high Galactic latitudes. In Fig. 14, we display the mean of all of the Gibbs samples, along with the standard deviation evaluated per COSMOGLOBE DR1 Gibbs sample, ${\sigma }_{{\beta }_{\text{s}}}^{\text{stat}}$, and the standard deviation over all COSMOGLOBE DR1 samples and Commander1 samples, ${\sigma }_{{\beta }_{\text{s}}}^{\text{sys+stat}}$. Thus, ${\sigma }_{{\beta }_{\text{s}}}^{\text{stat}}$ is the standard deviation when the input maps themselves are static, and ${\sigma }_{{\beta }_{\text{s}}}^{\text{sys+stat}}$ includes variations in the frequency maps themselves, corresponding to underlying instrumental effects, including bandpass, gain, noise characterization, and baseline estimation. The uncertainty due to white noise alone traces the high signal-to-noise regions of the polarized synchrotron, especially the prominent loops and spurs and the Fan region. At high Galactic latitudes, the standard deviation is 0.1, indicating that the posterior uncertainty is limited by the prior. To further quantify this, we performed additional Commander1 runs with priors of 𝒩(−3.2, 0.1) and 𝒩(−3.0, 0.1), and adopt the deviation from the mean in the fiducial analysis as an estimate of the extent of prior domination, ${\sigma }_{{\beta }_{\text{s}}}^{\text{prior}}\equiv (\langle {\beta }_{-3.0}\rangle -\langle {\beta }_{-3.2}\rangle )/2$. Adding this in quadrature gives the total uncertainty, as shown in the bottom panel of Fig. 14.

Fig. 14. COSMOGLOBE DR1 βs constraints. First: Mean spectral index over all samples. Second: Standard deviation over a single COSMOGLOBE DR1 sample. Third: Standard deviation over all Commander1 and COSMOGLOBE DR1 samples. Fourth: Difference due to prior choice added in quadrature.

The uncertainty across the entire Gibbs chain does not merely trace high signal-to-noise regions, and in fact there are variations that exceed the prior surrounding the Galactic center. These variations are primarily due to bandpass and gain uncertainties. Despite the relative increase in noise level when including instrumental effects, the Fan region and the Galactic loops are well constrained by the data. The Cygnus region is quite prominent in these maps, which is somewhat unexpected given its low polarized amplitude.

As a final quality check, we display the residuals with respect to the Commander1 sky model in Fig. 15, as well as the total scaled and normalized χ². In all bands but K, 30 GHz, and Ka, there are no visible artifacts due to instrumental uncertainty, with each map showing fluctuations consistent with the estimated white noise level calculated in the DR1 processing. The only exceptions are the Galactic plane and the ecliptic plane, regions that require more detailed Galactic modeling or have higher instrumental noise. In addition, many of the residuals visible in Fig. 3 have been reduced, demonstrating that polarized synchrotron spectral index variation provides meaningful improvements to the sky model fit. It is possible that there are residual systematics absorbed into variations in the spectral index, but any such effects must be smaller than our reported uncertainties. One exception to this is the Cygnus region, which in addition to having a high β_s with low uncertainty, has high residuals in the K band and each of the LFI bands. This suggests the presence of incompletely modeled temperature-to-polarization leakage.

Fig. 15. Residual maps and normalized χ2 in units of σ for the Commander1 analysis of the COSMOGLOBE DR1 data. Panels are organized by LFI channels (rows 1–2), WMAP (rows 3–7), and the total reduced χ2 (row 8).

The most salient remaining residuals are in K, 30 GHz, and Ka. Specifically, K-band and 30 GHz maps are anticorrelated surrounding the Galactic center, indicating tension between these two high signal-to-noise datasets. This could be due to either a genuine mis-modeling of the sky or incompletely modeled instrumental parameters. In particular, the signature is reminiscent of bandpass leakage corrections, which are shown, for example, in Fig. 9 of Svalheim et al. (2023b).

A more persistent residual is found in the Stokes QKa-band map. This feature has appeared in several different analyses, for example Fig. 4 of Svalheim et al. (2023a) and Fig. 8 of Weiland et al. (2022), but was not as clear without full removal of the poorly measured modes in the final map. The lack of this feature in the corresponding Stokes U map suggests that the effect is not a true Galactic effect, and is in some way due to instrumental processing, or un-modeled systematics; potential sources of this signal are discussed below.

5. Discussion and conclusion

We have presented a new state-of-the-art model of polarized synchrotron using the COSMOGLOBE DR1 data products, for the first time combining all Planck LFI and WMAP frequency maps and thereby leveraging the full statistical power of these data. The polarized synchrotron map as delivered by Watts et al. (2023a) has an effective white noise level of 3.4 μK at 2° FWHM, a 29% improvement over the white noise levels of the WMAP K-band and LFI 30 GHz maps. This model characterizes the polarized sky signal to within 5 μK in all bands. We have also reproduced the previously reported B-to-E-mode ratio in the polarized synchrotron power spectrum. Furthermore, we have shown that physically reasonable spectral indices can be recovered using a variety of methods, all of which are consistent with the two methods presented in this paper and with previously published results from CLASS and QUIJOTE, confirming variation in β_s along the Galactic plane and steepening at high Galactic latitudes.

Our improved processing of the WMAP and LFI data in conjunction with improved polarized synchrotron modeling has allowed us to dig deeper into the underlying differences between the two datasets. We have shown agreement between the β_s values derived with the T–T plot using K/30 GHz and K/Ka data combinations, further evidence that these datasets are now consistent with each other. The remaining unexplained differences between the sky maps and our derived model are primarily in the Galactic center, and have morphologies consistent with bandpass errors and SED complexities.

The least well-understood residual is in the Ka band, specifically Stokes Q. Of all the known instrumental effects in both WMAP and Planck LFI data, the only one that morphologically resembles this residual is the bandpass correction in the K band, which was previously presented in Fig. 12 of Watts et al. (2023a). The bandpass correction term depends solely on the laboratory measurements presented by Jarosik et al. (2003). The high statistical weight of the K band for determining polarized synchrotron could easily lead to an under-subtraction in the Ka-band maps. Conversely, the nominal bandpass correction to the Ka band itself is negligible, as all of the Ka-band radiometers were reported to have nearly identical bandpasses. A full accounting of this effect, including sampling of bandpass differences in the WMAP TOD not performed by Watts et al. 2023a, will be performed in future work.

The estimation of β_s is difficult precisely because the quantity is dependent on differences between different frequencies, which can often exacerbate systematic effects and processing choices. The differences between a sky model and datasets is best determined in the time streams – much of the improvement demonstrated in this analysis would not have been possible through a purely map-based analysis. In addition to improving the data quality, this approach allowed for a natural end-to-end propagation of errors, and demonstrated that some of the brightest regions of the sky do not in fact have well-determined spectral indices, despite their high signal-to-instrumental noise ratio. As increasingly stringent uncertainty constraints are becoming necessary to measure a nonzero tensor-to-scalar ratio, the use of information from different experiments with complementary observation strategies will increasingly become a necessity. Future joint analyses, including of QUIJOTE and CLASS data, will continue to improve our knowledge of the polarized sky.

Acknowledgments

We thank Dr. Ricardo Génova-Santos, whose suggestions improved this work. We thank the entire Planck and WMAP teams for invaluable support and discussions, and for their dedicated efforts through several decades without which this work would not be possible. The current work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement numbers 819478 (ERC; COSMOGLOBE), 772253 (ERC; BITS2COSMOLOGY), and 101007633 (Marie Skłodowska-Curie, CMB-INFLATE). In addition, the collaboration acknowledges support from RCN (Norway; grant no. 274990). The research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. Some of the results in this paper have been derived using the healpy and HEALPix http://healpix.sf.net packages (Górski et al. 2005; Zonca et al. 2019). This work made use of Astropy: http://www.astropy.org a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

References

All Tables

All Figures

Fig. 1. Polarized intensity and white noise levels of the WMAP K band (top), Planck 30 GHz band (middle), and synchrotron amplitude from the COSMOGLOBE Gibbs chain (bottom), all evaluated at 30 GHz with a resolution of 72′. All maps have been scaled to a common frequency of 30 GHz. The leftmost column is a 10° wide square centered on the low signal-to-noise region highlighted by a red square in the middle column’s panels. The middle column shows the polarized total amplitude, P. The rightmost column shows the rms noise for the K-band and 30 GHz frequency maps and the COSMOGLOBE DR1 synchrotron amplitude map.

In the text

	Fig. 2. Ratio between the COSMOGLOBE and BEYONDPLANCK polarized synchrotron amplitude standard deviations per HEALPixNside = 1024 pixel. The posterior standard deviation is computed over all Gibbs samples for both the COSMOGLOBE and BEYONDPLANCK chains.
In the text

	Fig. 3. Frequency map minus sky model with respect to the COSMOGLOBE DR1 sky model evaluated at 5°. COSMOGLOBE maps are labeled CG, Planck 2018 maps are labeled PR3, and the legacy WMAP9 maps are labeled WMAP. Each of the WMAP9, Planck PR3, and COSMOGLOBE maps have had the COSMOGLOBE DR1 sky model evaluated at their respective bandpass and the resolution subtracted.
In the text

Fig. 4. Half-mission cross spectra for polarized synchrotron emission as evaluated for Planck PR3 (blue) and COSMOGLOBE DR1 (black). Filled circles correspond to E modes, while empty circles correspond to B modes. The thick solid and dashed lines are the best-fit E- and B-mode power law fits to synchrotron, and the thin black line is the Λ cold dark matter prediction for E modes. The lowest bin is not included in either fit.

In the text

	Fig. 5. Spectral index regions as defined by Fuskeland et al. (2014). The most prominent point sources are masked out and shown in gray circular areas.
In the text

Fig. 6. T–T plots for Stokes Q and U maps of the COSMOGLOBE DR1 K- band versus the COSMOGLOBE DR1 30 GHz (black) and the WMAP9 K- band versus Planck PR3 30 GHz pixels (red) for all regions. The horizontal (solid and dotted) lines indicate the corresponding inverse variance weighted values of the spectral index, averaged over rotation angle, and in the COSMOGLOBE case also weighted over samples. The corresponding βs values can be found in Fig. 8.

In the text

Fig. 7. Synchrotron spectral index as a function of rotation angle, computed using a T–T plot between the COSMOGLOBE DR1 K-band and the COSMOGLOBE DR1 30 GHz pixels (black) compared to the spectral index using the WMAP9 K-band and Planck PR3 30 GHz (red) pixels for all regions. The horizontal (solid and dotted) lines indicates the corresponding inverse variance weighted values of the spectral index, averaged over rotation angle, and in the COSMOGLOBE case also samples.

In the text

	Fig. 8. Spatial variation in the synchrotron spectral index, computed using T–T plots between the (from top to bottom) WMAP9 K-band and Planck PR3 30 GHz, WMAP9 K-band and Planck PR4 30 GHz, WMAP9 K-band and BEYONDPLANCK 30 GHz, and COSMOGLOBEK-band and COSMOGLOBE 30 GHz maps. The spectral index is inverse variance weighted over rotation angle, and in the COSMOGLOBE case also samples.
In the text

Fig. 9. Synchrotron spectral index as a function of region number, computed using T–T plots between the WMAP9 K-band and Planck PR3 30 GHz (red), WMAP9 K-band and Planck PR4 30 GHz (blue), WMAP9 K-band and BEYONDPLANCK 30 GHz (purple), and COSMOGLOBEK-band and COSMOGLOBE 30 GHz (black) maps. The spectral index is inverse variance weighted over rotation angles, and samples. The horizontal line in the high latitude regions corresponds to the estimated spectral index values from the Planck PR3 likelihood analysis (Planck Collaboration V 2020).

In the text

	Fig. 10. Differences between the COSMOGLOBE and official WMAP 23 GHz (top row) and between COSMOGLOBE and Planck PR3 30 GHz (bottom row) maps, with spectral index regions overlaid. Columns show Stokes Q and U parameters.
In the text

	Fig. 11. Synchrotron spectral index computed using the T–T plot method with the COSMOGLOBE DR1 K-band (23 GHz) and 30 GHz data versus COSMOGLOBE DR1 K-band and Ka-band (33 GHz) data.
In the text

	Fig. 12. Spatial variation in the synchrotron spectral index, computed using the T–T plot method between the COSMOGLOBE DR1 K- and Ka-band maps. The spectral index is inverse variance weighted over the rotation angle and samples.
In the text

	Fig. 13. Spectral index estimates from COSMOGLOBE using T–T plots with K and Ka maps and Gibbs sampling using all frequencies, and from QUIJOTE and CLASS using T–T plots. Only regions with ≥75% of the pixels within a given region are displayed. Red points and black points are both calculated using COSMOGLOBE frequency maps, but are obtained with different techniques.
In the text

	Fig. 14. COSMOGLOBE DR1 βs constraints. First: Mean spectral index over all samples. Second: Standard deviation over a single COSMOGLOBE DR1 sample. Third: Standard deviation over all Commander1 and COSMOGLOBE DR1 samples. Fourth: Difference due to prior choice added in quadrature.
In the text

	Fig. 15. Residual maps and normalized χ2 in units of σ for the Commander1 analysis of the COSMOGLOBE DR1 data. Panels are organized by LFI channels (rows 1–2), WMAP (rows 3–7), and the total reduced χ2 (row 8).
In the text