© The Authors 2023

1. Introduction

One of the primary goals of future cosmic microwave background (CMB) missions such as CMB-S4 (Abazajian et al. 2019), Simons Observatory (Ade et al. 2019), and LiteBIRD (Sugai et al. 2020) is to detect the hypothesized signal from primordial gravitational waves, or to constrain the tensor-to-scalar ratio r. In order for this goal to be reached, it is essential to characterize polarized Galactic emission at the 10−100 nK level.

Two emission mechanisms dominate the polarized Galactic signal at microwave frequencies: synchrotron emission below ∼70 GHz and thermal dust emission above ∼70 GHz (e.g., Page et al. 2007; Planck Collaboration X 2016). Recent measurements from Planck have shown that the polarization fractions of these two components can reach 40 and 20%, respectively (Planck Collaboration XII 2020), and their rich morphologies are connected by the Galactic magnetic field. Other emission mechanisms, such as free-free and anomalous microwave emission (AME), are important contributions to the total intensity observed in the same frequency range, but no significant polarization has been detected yet (e.g., Planck Collaboration XXV 2016), although a small free-free contribution is predicted on very small angular scales near the edges of H II regions (Rybicki & Lightman 1985).

Anomalous microwave emission was serendipitously detected while analyzing COBE–DMR and other data sets in the 10−60 GHz range (Kogut et al. 1996; Leitch et al. 1997; de Oliveira-Costa et al. 2004). This emission is highly correlated with Galactic thermal dust morphology and has been detected within specific clouds such as the Perseus molecular cloud (Watson et al. 2005; Planck Collaboration XXV 2016; Tibbs et al. 2010; Génova-Santos et al. 2015a), λ Orionis (Bell et al. 2019; Cepeda-Arroita et al. 2021), ρ Ophiuchus (Casassus et al. 2008; Arce-Tord et al. 2020), among others (Planck Collaboration XV 2014), and potentially even in external galaxies (Murphy et al. 2010, 2018; Battistelli et al. 2019). However, the polarized properties of this emission is still uncertain as physical models predict varying levels of potential polarized emission (Dickinson et al. 2018).

Ignoring this component during Galactic foreground removal can have a detrimental effect on polarized CMB determination. As demonstrated by Remazeilles et al. (2016), neglecting AME models with a polarization fraction as low as ∼1% can cause a nonnegligible bias on the measured tensor-to-scalar ratio for a wide range of proposed future missions.

The currently favored model for AME is spinning dust emission, in which the radiation is due to small, rotating dust grains with an electric or magnetic dipole moment (Draine & Lazarian 1998a). The total power P_ν emitted at frequency ν from a single grain rotating at frequency ω = 2πν, electric or magnetic dipole moment μ, and with angle θ between its dipole moment and rotation axis is given by

$$\begin{array}{c}\hfill P_{\nu }=\frac{2}{3}\frac{{\omega }^4{\mathrm{\mu }}^2{sin}^2\theta }{c^3}·\end{array}$$(1)

Much of the theoretical work on spinning dust emission has been devoted to predicting the distribution of grain angular velocities present in the interstellar medium, which depends on properties such as grain size and shape as well as the local environment, such as the local gas density, temperature, and ionization (Draine & Lazarian 1998a; Hoang et al. 2010; Silsbee et al. 2011). These properties are poorly constrained, so theoretical predictions for the frequency spectrum of the AME are highly uncertain.

The details of this model have been combined in the SpDust2 code to produce AME spectra for varying interstellar medium conditions (Hoang et al. 2010; Silsbee et al. 2011). In Planck Collaboration IX (2016), the Commander (Eriksen et al. 2008) analysis had to combine two SpDust2 spectra in order to get a good fit to the observed AME. While theoretical studies have attributed two component AME models to emission from the cold and warm neutral media (Hoang et al. 2010; Ysard et al. 2010), a recent study found no link between the AME spectrum and the fraction of HI in the cold neutral phase (Hensley et al. 2022), casting doubt on this interpretation. In this paper, we therefore seek constraints on AME properties without imposing strong assumptions on its frequency spectrum.

The polarization of spinning dust emission depends on how efficiently they align with the local magnetic field. Hoang et al. (2013) argue that starlight polarization features in HD 197770 and HD 147933-4 could be caused by weakly aligned polycyclic aromatic hydrocarbons (PAHs), corresponding to polarization from spinning PAHs of ≲1% for frequencies ∼20 GHz. Additionally, Hoang & Lazarian (2016) calculated that iron nanoparticles can efficiently align with the Galactic magnetic field, yielding a wide range of potential polarization fractions of their microwave spinning dust emission. On the other hand, Draine & Hensley (2016) argue that quantum suppression of energy dissipation interferes with alignment of very small grains, leading to extremely small polarization fractions at microwave frequencies.

Though the spinning dust model is favored due to intensity measurements, other models have not been definitively ruled out. A theory of emission from thermal fluctuations in magnetic dust grains has been proposed by Draine & Lazarian (1999), and subsequently revisited by Draine & Hensley (2013). They identified resonance behavior that can produced peaked, AME-like spectra dependent on the magnetic properties and shapes of the grains. Predictions for polarization fractions in such resonances range from ∼5 to 40%, in conflict with current best limits on AME polarization.

Theoretical predictions for the low-level polarized emission from AME are supported by most published observations. Using QUIJOTE data focused on molecular complex W43r, Génova-Santos et al. (2017) found rigid upper limits for polarized AME at 16.7, 22.7, and 40.6 GHz of 0.39, 0.52 and 0.22% respectively. The most recent analysis of QUIJOTE data by Poidevin et al. (2019) concluded that the polarization fraction is indeed consistent with zero within the Taurus molecular cloud. This analysis found the polarization fraction of AME to be less than 4.2% (95% confidence), assuming that all of the observed polarized emission were to come from AME. Battistelli et al. (2006) measured the total emission within the Perseus region to be polarized at $P=3.4_{-1.9}^{+1.5}\%$ at 11 GHz (95% confidence interval), while Planck Collaboration XXV (2016) derives a 2σ AME polarization percentage limit of less than 1.6% within the same region. The differing results are likely due to the handling of polarized synchrotron emission, and this emphasizes the importance of synchrotron marginalization.

Macellari et al. (2011) carried out a template based cross-correlation analysis with standard foreground templates to estimate the amplitude of the Galactic components in both intensity and polarization. Correlations were determined for synchrotron, dust, and free-free emission within the WMAP 5-year maps using Haslam et al. (1982), the WMAP 94 GHz dust prediction, and the Finkbeiner (2003) Hα maps as templates, respectively. With this method they found the polarization of both dust and free-free emission to be consistent with zero for the all-sky analysis. Concerning AME, this results in a full sky AME polarization fraction at WMAP K-band of less than 5% (95% confidence). They also highlight dust-correlated polarization fractions at a 2σ level within some regions of the sky within the K-band, while noting slightly negative dust polarizations (negative correlation between total intensity template and polarized emission) indicating a degeneracy between sky components.

In this work we add to earlier constraints in three key ways. Firstly, rather than investigating a specific region, we fit the level of polarized dust correlated emission using nearly the full sky. Secondly, we make no assumptions about the AME SED to avoid bias from an imperfect spectral fit. Instead, we make a strong spatial assumption, namely that AME is perfectly correlated (or anti-correlated) with thermal dust emission. Finally, we employ the Bayesian BEYONDPLANCK framework (BeyondPlanck Collaboration 2023) to marginalize over the uncertainty in the other Galactic foregrounds and systematic effects. This approach attacks the CMB analysis problem by fitting instrumental, astrophysical and cosmological parameters all jointly within an integrated Gibbs sampling framework. The synchrotron spectral index and amplitude are independently redetermined within the current paper, to allow joint fits with the AME parameters. Technically speaking, marginalization over systematic uncertainties is implemented by analysing an ensemble of different BEYONDPLANCK Low-Frequency Instrument (LFI) map realizations within the sampling scheme, thereby propagating low-level uncertainties from the BEYONDPLANCK mapmaking procedure. This is the first example of how systematic errors may be propagated into postproduction analyses using Monte-Carlo ensembles with LFI data.

The paper is organized as follows: Sect. 2 describes the data used in this analysis. An overview of the sky model and sampling techniques are discussed in Sects. 3 and 4 respectively, and in Sect. 5 the results of the Gibbs sampler on the BEYONDPLANCK data are discussed. Concluding remarks and outlook are reported in Sect. 6. All software is made publicly available¹.

2. Data

Because AME is only observed between 10 and 60 GHz in intensity, we restrict our data selection to the 20−70 GHz frequency range. As such, our main observations are the polarization observations from the Planck LFI (Planck Collaboration II 2020) and WMAP (Bennett et al. 2013) experiments. For LFI, we adopt the BEYONDPLANCK 30, 44, and 70 GHz frequency sky maps² (BeyondPlanck Collaboration 2023), which come in the form of an ensemble of samples, each of which represents a different realization of systematic effects for the frequency channel in question. These include, but are not limited to, corrections for gain variations (Gjerløw et al. 2023), correlated noise (Ihle et al. 2023), bandpass leakage (Svalheim et al. 2023a) and far sidelobe contamination (Galloway et al. 2023). In order to properly marginalize over these uncertainties, we draw maps from an ensemble of 1000 different map realizations.

For WMAP, we include the K- (22.8 GHz), Ka- (33.0 GHz), Q- (40.6 GHz), and V-band (60.8 GHz) channels. The W-band (94 GHz) data are however excluded due to a high level of systematic residuals (Bennett et al. 2013).

In addition to these main observations, we utilize the Planck DR4 353-GHz polarization map (Planck Collaboration Int. LVII 2020) as a tracer of polarized dust emission. We use the raw DR4 353-GHz map, with no CMB subtraction, noting that this channel is strongly dust dominated. Finally, in order to derive an estimate of the polarization fraction p_AME of AME we need an estimate of the AME intensity. For this, we adopt the Planck 2015 AME component map (Planck Collaboration X 2016), rather than the corresponding BEYONDPLANCK AME maps (Andersen et al. 2023), as the former was derived without strong spatial priors.

All maps are smoothed to a common angular resolution of 1° FWHM Gaussian beam, and discretized using the HEALPix³ pixelization format (Górski et al. 2005) with a resolution parameter of N_side = 64. We adopt the cosmological convention for the definition of the Stokes Q and U parameters, which differs from the IAU convention only in the sign of U (Górski et al. 2005). Rayleigh-Jeans brightness temperature units (μK_RJ) are utilized for the entirety of this work, and full bandpass integration for both WMAP and BEYONDPLANCK LFI maps are taken into account during the fitting procedures.

To obtain a reliable estimate of the level of detectable polarized dust emission below 60 GHz we need to mask portions of the sky which may bias our results. At the same time, it is highly desirable to include as much of the sky as possible, to maximize our signal-to-noise ratio. We therefore choose to primarily define our primary analysis mask on the central Galactic region defined in Svalheim et al. (2023b), but additionally we mask out a few compact regions with very high χ²’s, including Tau-A. A total of f_sky = 91% is included in the final analysis.

3. Sky model

3.1. Modelling considerations

In this paper, we make no assumptions regarding the AME spectral energy density (SED). Instead, we note from the Planck 2015 analysis that AME in total intensity is highly correlated with thermal dust emission (Planck Collaboration X 2016). This motivates us for now to assume that any detectable polarized AME also has the same morphology as polarized thermal dust emission as measured at 353 GHz. If the grains emitting the AME are aligned by the same magnetic field aligning those responsible for FIR polarization, and if the alignment is in the usual sense, that is the short axis parallel to the magnetic field B, then the polarization direction at AME frequencies is the same as that observed at 353 GHz (Draine & Lazarian 1998b).

The polarization fraction of dust emission is a function of grain properties (alignment, shape, size), density structure, and magnetic field geometry. As shown in Planck Collaboration XII (2020) using the Planck full-sky maps, neither grain alignment nor the dust temperature drive variations in the polarization fraction over diffuse and translucent sightlines (for column densities N_H < 8 × 10²¹ cm⁻², corresponding to 98% of the sky). The highest column densities where this conclusion breaks down lie, for the most part, within the portion of the sky which is masked out here. This conclusion is supported by Clark & Hensley (2019), who find that with H I structures alone are able to reproduce the Planck 353 GHz polarization maps at degree scales with high accuracy over most of the sky. Therefore, if the grains producing the AME and the far infrared continuum emission are aligned with the same magnetic field, we expect the polarization fraction of each to be similar. Nevertheless, it is possible that, unlike the larger grains, the alignment efficiency of nanoparticles varies substantially across the sky. In this case, the polarization fractions of the two emission mechanisms would be less correlated and AME polarization fractions in excess of those derived in this work would be permitted However, there is no strong observational evidence for this scenario and we do no consider it further in our analysis.

A known issue in AME polarization studies is the presence of polarized synchrotron emission (Planck Collaboration XXV 2016). Additionally, polarized thermal dust emission is nonnegligible even at AME frequencies. In principle, we could disentangle these distinct Galactic components by using a template for polarized synchrotron emission in a way analogous to the dust emission. However, no suitable high signal-to-noise ratio polarized synchrotron templates currently exist. As an example, one might consider using the polarized synchrotron component derived by the SMICA algorithm in Planck Collaboration IV (2020). However, this analysis explicitly assumes that the polarized thermal dust emission “vanishes at 30 GHz”, and therefore all polarized emission at 30 GHz is assigned to be exclusively synchrotron emission. This would lead to erroneously tight constraints on polarized AME emission in the current analysis. In the absence of a high signal-to-noise ratio template of polarized synchrotron emission, we adopt a simple parametric power-law model to fit the synchrotron component in each pixel, and marginalize over the free parameters.

3.2. Parametric model

We model the observed emission d_ν, p in pixel p at frequency ν as the sum of the true sky signal s_ν,  p and the noise n_ν,  p, given by

$$\begin{array}{c}\hfill d_{\nu ,p}=s_{\nu ,p}+n_{\nu ,p}.\end{array}$$(2)

The polarized sky signal, s_ν, p, is usually described by three components in the microwave regime in addition to AME, namely thermal dust emission, synchrotron emission, and CMB (Planck Collaboration XXV 2016; Planck Collaboration IV 2018; Svalheim et al. 2023b). In this analysis, however, we do not include CMB because this contribution is known to be small over the frequency range and angular scales of interest (Planck Collaboration V 2020), and including unconstrained degrees of freedom typically leads to nonphysical degeneracies. Instead, we perform dedicated sensitivity tests in Sect. 5.4 by adding simulated CMB fluctuations to the existing sky maps, and show that this has a minor impact on final results. The sky model may therefore be written out in the following form,

$$\begin{array}{c}\hfill s_{\nu ,p}=s_{\mathrm{s},\nu ,p}+s_{\mathrm{AME},\nu ,p}+s_{\mathrm{d},\nu ,p}.\end{array}$$(3)

Synchrotron emission is modeled in terms of a power-law SED with a free amplitude per Stokes parameter and pixel, $a_{\mathrm{s},p}^{Q,U}$, and a spectral index, β_s, p, per pixel, but common for Q and U,

$$\begin{array}{c}\hfill s_{\mathrm{s},\nu ,p}^{Q,U}=a_{\mathrm{s},p}^{Q,U}{\left(\frac{\nu }{{\nu }_{0,\mathrm{s}}}\right)}^{{\beta }_{\mathrm{s},p}}.\end{array}$$(4)

We choose a reference frequency of ν_0, s = 30 GHz for synchrotron emission.

For thermal dust emission, we adopt a single modified blackbody (MBB) SED,

$$\begin{array}{c}\hfill s_{\mathrm{d},\nu ,p}^{Q,U}={\mathcal{T}}_{\mathrm{d},p}^{Q,U}{\left(\frac{\nu }{{\nu }_{0,\mathrm{d}}}\right)}^{{\beta }_{\mathrm{d}}+1}\frac{{\mathrm{e}}^{h{\nu }_{0,\mathrm{d}}/k{T}_{\mathrm{d}}}-1}{{\mathrm{e}}^{h\nu /k{T}_{\mathrm{d}}}-1},\end{array}$$(5)

where h and k are the Planck and Boltzmann constants, respectively; ν_0, d is the thermal dust reference frequency, which is taken to be 353 GHz; β_d is the thermal dust spectral index; and T_d = 19.6 K is the dust temperature; and 𝒯_d,  p is the Planck 353 GHz map. We marginalize over the thermal dust spectral index using the same prior as in the main BEYONDPLANCK analysis (β_d = 1.62 ± 0.04; Svalheim et al. 2023b) by drawing a new value for β_d in each Markov chain iteration.

Finally, we assume that polarized AME is spatially perfectly correlated with polarized thermal dust emission, and we adopt the Planck DR4 353 GHz map as a spatial template for polarized AME. Allowing for a single multiplicative amplitude, a_AME,  ν, for both Stokes parameters, our model reads

$$\begin{array}{c}\hfill s_{\mathrm{AME},\nu ,p}^{Q,U}=a_{\mathrm{AME},\nu }{\mathcal{T}}_{\mathrm{d},p}^{Q,U}.\end{array}$$(6)

However, to avoid a perfect degeneracy between the spatial synchrotron template and the AME template amplitudes, a maximum of N_band − 1 template amplitudes can be fit. To break this degeneracy, we fix a_AME,  ν = 0 for the WMAP V-band (61 GHz) and LFI 70 GHz channels, where the AME-to-thermal-dust ratio in intensity is 0.033 and 0.008, respectively, in the Planck 2015 model (Planck Collaboration X 2016).

Collecting terms, we arrive at the final complete data model,

$$\begin{array}{c}\hfill s_{\nu ,p}^{Q,U}=a_{\mathrm{s},p}^{Q,U}{\left(\frac{\nu }{{\nu }_{0,\mathrm{s}}}\right)}^{{\beta }_{\mathrm{s},p}}+{\mathcal{T}}_{\mathrm{d},p}^{Q,U}({\left(\frac{\nu }{{\nu }_{0,\mathrm{d}}}\right)}^{{\beta }_{\mathrm{d}}+1}\frac{{\mathrm{e}}^{h{\nu }_{0,\mathrm{d}}/k{T}_{\mathrm{d}}}-1}{{\mathrm{e}}^{h\nu /k{T}_{\mathrm{d}}}-1}+a_{\mathrm{AME},\nu })·\end{array}$$(7)

In this model we fit $a_{\mathrm{s},p}^{Q,U}$ and β_s, p per pixel, and a_AME, ν as a single value across the sky. The synchrotron amplitude is determined independently in Stokes Q and U, while both β_s and a_AME, ν are fit jointly in Q and U. The primary goal of the following analysis is to constrain the a_AME, ν coefficients, and the other parameters are simply considered to be nuisance parameters that we marginalize over.

3.3. From AME template amplitudes to polarization fraction

The AME polarization fraction p_AME is defined as the ratio of the polarized (P_{AME,  ν,  p}) and total AME intensities (I_{AME,  ν,  p}) at each frequency ν and pixel p, i.e. P_{AME,  ν,  p}/I_{AME,  ν,  p}.

To translate the vector of template amplitudes, a_AME, to an AME polarization fraction, p_AME,  ν, we use Eq. (6) to associate polarized AME with polarized thermal dust emission,

$$\begin{array}{c}\hfill P_{\mathrm{AME},\nu ,p}=a_{\mathrm{AME},\nu }\sqrt{{{\mathcal{T}}_{\mathrm{d},p}^Q}^2+{{\mathcal{T}}_{\mathrm{d},p}^U}^2}=a_{\mathrm{AME},\nu }{P}_{353,p}.\end{array}$$(8)

Following the results from Planck 2015, we can describe I_{AME,  ν,  p} as the product of the frequency dependence f_AME(ν) and the AME amplitude map A_AME,  p, yielding

$$\begin{array}{cc}\hfill p_{\mathrm{AME},\nu ,p}& =a_{\mathrm{AME},\nu }{P}_{353,p}/I_{\mathrm{AME},\nu ,p}\hfill \end{array}$$(9)

$$\begin{array}{cc}& =p_{353,p}\frac{a_{\mathrm{AME},\nu }}{f_{\mathrm{AME}}(\nu )}\frac{I_{353,p}}{A_{\mathrm{AME},p}},\hfill \end{array}$$(10)

where p_353,  p is equivalently defined as the polarization fraction at 353 GHz, i.e. p_353,  p = P_353,  p/I_353,  p. We estimate the maximum AME polarization fraction $p_{\mathrm{AME},\nu }^{\max }$, corresponding to the magnetic field in the plane of the sky as

$$\begin{array}{c}\hfill p_{\mathrm{AME},\nu }^{\max }=p_{353}^{\max }\frac{a_{\mathrm{AME},\nu }}{f_{\mathrm{AME}}(\nu )}\langle \frac{I_{353}}{A_{\mathrm{AME}}}\rangle ,\end{array}$$(11)

where $p_{353}^{\max }$ is the maximum dust polarization fraction at 353 GHz and $\langle \frac{I_{353}}{A_{\mathrm{AME}}}\rangle$ is the ratio of the 353 GHz dust emission to 22.8 GHz AME emission, assumed to be constant over the sky in analogy with our template analysis for polarization. We adopt $p_{353}^{\max }=22\%$ (Planck Collaboration XII 2020), and $⟨I_{353}/A_{\mathrm{AME}}⟩=2.5_{-0.7}^{+0.2}$ and f_AME(ν) are taken from the Planck 2015 data products (Planck Collaboration X 2016). This allows us to estimate $p_{\mathrm{AME},\nu }^{\max }$ directly from our fit amplitudes a_AME,  ν.

We note that Eq. (10) provides a direct linear estimate of the AME polarization fraction in terms of the fitted AME template amplitudes, which themselves are fitted linearly to the raw data. As such, Eq. (10) supports simple Gaussian error propagation without explicit noise debiasing, which otherwise often is a problem for polarization fraction estimates with low signal-to-noise data; readers can refer to Montier et al. (2015), for example. With the current formulation, this issue is circumvented through linear association with high signal-to-noise thermal dust emission estimates that themselves have been noise debiased.

4. Sampling algorithms

4.1. Posterior distribution and Gibbs sampling

Following the general approach in the BEYONDPLANCK environment, we employ Bayesian sampling methods to estimate the various free parameters in our model. Let us define the set of all free parameters in Eq. (7) as ω = {a_AME, a_s, β_s, β_d}. The appropriate posterior distribution is then given by Bayes’ theorem,

$$\begin{array}{c}\hfill P(\mathit{\omega }\mid \mathit{d})=\frac{P(\mathit{d}\mid \mathit{\omega })P(\mathit{\omega })}{P(\mathit{d})}\propto \mathcal{L}(\mathit{\omega })P(\mathit{\omega }),\end{array}$$(12)

where P(d ∣ ω)≡ℒ(ω) is the likelihood function, P(ω) is the prior, and P(d) is a normalization factor which we disregard here as it is independent of the parameters ω.

In order to map out the probability distribution of our multidimensional parameter space in a computationally efficient way, we adopt Gibbs sampling (e.g., Gelman et al. 2003). In the Gibbs sampling framework we cycle through each parameter independently, drawing each parameter from a distribution conditional on the values of the other parameters; see, e.g., BeyondPlanck Collaboration (2023) for details of the much larger BEYONDPLANCK Gibbs sampler. In the present work, the Gibbs sampler takes the following form,

$$\begin{array}{cc}& ({\mathit{a}}_{\mathrm{AME}},{\mathit{a}}_{\mathrm{s}})\leftarrow P(a_{\mathrm{AME}},{\mathit{a}}_{\mathrm{s}}|\mathit{d},{\beta }_{\mathrm{s}},{\beta }_{\mathrm{d}})\hfill \end{array}$$(13)

$$\begin{array}{cc}& {\beta }_{\mathrm{s}}\leftarrow P({\beta }_{\mathrm{s}}|\mathit{d},{\mathit{a}}_{\mathrm{AME}},{\mathit{a}}_{\mathrm{s}},{\beta }_{\mathrm{d}})\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\beta }_{\mathrm{d}}\leftarrow P({\beta }_{\mathrm{d}}).\hfill \end{array}$$(15)

The two first steps will be detailed below, while the third step just denotes sampling the thermal dust spectral index from its Gaussian prior, as discussed in Sect. 3.2.

4.2. Amplitude sampling

To sample from the conditional distribution for the amplitudes in Eq. (13), we note from Eq. (7) that the signal model is linear in both a_AME and a_s. The data model in Eq. (2) may therefore be written in the following compact vector form,

$$\begin{array}{c}\hfill d_{\nu ,p}={\mathsf{T}}_{\nu }\mathit{a}+s_{\mathrm{d},\nu }+n_{\nu ,p},\end{array}$$(16)

where

$$\begin{array}{cc}& \mathit{a}\equiv {\left[\begin{array}{cccccc}{a}_{\mathrm{s},0}& \cdots & a_{\mathrm{s},{\mathrm{N}}_{\mathrm{pix}}-1}& a_{\mathrm{AME},1}& \cdots & a_{\mathrm{AME},{\mathrm{N}}_{\mathrm{bands}}-1}\end{array}\right]}^T\hfill \\ \hfill & {\mathsf{T}}_{\nu }\equiv \left[\begin{array}{cccccc}{(\frac{\nu }{{\nu }_{0,\mathrm{s}}})}^{{\beta }_{\mathrm{s},0}}& 0& 0& {\delta }_{\nu ,1}{\mathcal{T}}_{\mathrm{d},0}& \cdots & {\delta }_{\nu ,n-1}{\mathcal{T}}_{\mathrm{d},0}\\ 0& \ddots & 0& ⋮& \cdots & ⋮\\ 0& 0& {(\frac{\nu }{{\nu }_{0,\mathrm{s}}})}^{{\beta }_{\mathrm{s},N_{\mathrm{pix}}-1}}& {\delta }_{\nu ,1}{\mathcal{T}}_{\mathrm{d},N_{\mathrm{pix}}-1}& \cdots & {\delta }_{\nu ,n-1}{\mathcal{T}}_{\mathrm{d},N_{\mathrm{pix}}-1}\end{array}\right].\hfill \end{array}$$

In these expressions, indices for Stokes Q and U parameters are suppressed, but we note that all frequency or component maps should be considered as stacked vectors, a = {a^Q, a^U}, and matrices are generalized accordingly.

With this notation, d_ν − T_νa − s_d, ν = n_ν, which we assume to be Gaussian distributed with a covariance matrix N_ν. Therefore, the likelihood may be written in terms of a standard multivariate Gaussian,

$$\begin{array}{c}\hfill \mathcal{L}(\omega )\propto {\mathrm{e}}^{-\frac{1}{2}{\sum }_{\nu }{({\mathit{d}}_{\nu }-{\mathsf{T}}_{\nu }\mathit{a}-{\mathit{s}}_{\mathrm{d},\nu })}^t{\mathsf{N}}_{\nu }^{-1}({\mathit{d}}_{\nu }-{\mathsf{T}}_{\nu }\mathit{a}-{\mathit{s}}_{\mathrm{d},\nu })}\equiv {\mathrm{e}}^{-\frac{1}{2}{\chi }^2(\mathit{a},\beta )}.\end{array}$$(17)

When interpreted as a function of a only, this is a strict Gaussian distribution from which a sample may be drawn through the following set of linear equations (see Appendix A of BeyondPlanck Collaboration 2023 for a pedagogical derivation),

$$\begin{array}{c}\hfill {\displaystyle \sum _{\nu }\left({\mathsf{T}}_{\nu ,j}^T{\mathsf{N}}_{\nu }^{-1}{\mathsf{T}}_{\nu ,j}\right)\mathit{a}=\sum _{\nu }({\mathsf{T}}_{\nu ,j}^T{\mathsf{N}}_{\nu }^{-1}({\mathit{d}}_{\nu }-{\mathit{s}}_{\mathrm{d},\nu })+{\mathsf{T}}_{\nu ,j}^T{\mathsf{N}}_{\nu }^{-1/2}{\mathit{\eta }}_j).}\end{array}$$

Here η_j is a vector of standard Gaussian random variates, N(0, 1). This equation may be solved efficiently through standard Conjugate Gradient techniques (Shewchuk 1994).

4.3. Sampling the synchrotron spectral index

Sampling the synchrotron spectral index has been an active topic of analysis throughout the history of both WMAP and Planck. Both the WMAP and Planck Collaborations note difficulties in independently determining β_s in polarization and therefore adopted the determination from measurements in total intensity (Bennett et al. 2013; Planck Collaboration IV 2020). Dunkley et al. (2009) assumed a prior of β_s = −3.0 ± 0.3 from Rybicki & Lightman (1985), finding a mean index of −3.03 ± 0.04 in regions of the sky with high signal-to-noise. More recent analyses including S-PASS data find spatial variations in the spectral index ranging between β_s ∼ −2.8 at low Galactic latitudes, and β_s ∼ −3.1 at higher Galactic latitudes, though both are consistent with β_s ∼ −3.1. Additionally, recent work with C-BASS data shows the average synchrotron spectral index between 4.76 and 22.8 GHz to be β_s = −3.1 ± 0.02 (Harper et al. 2022). Here we explore a variety of priors for with −3.3 < ⟨β_s⟩<  − 3.0, and assume a nominal prior of β_s = −3.1 ± 0.1, following the results of Svalheim et al. (2023b).

To sample from Eq. (14), we employ a Metropolis Markov chain Monte Carlo sampler (Metropolis et al. 1953; BeyondPlanck Collaboration 2023). Specifically, let β_j represent the jth sample of β_s in a Markov chain (implicitly suppressing the synchrotron label for readability), and let T(β_j + 1|β_j) be a symmetric stochastic proposal probability distribution, i.e., T(β_j + 1|β_j) = T(β_j|β_j + 1); we adopt a simple Gaussian distribution with a mean of β_j and a (tunable) standard deviation of σ.

The Metropolis algorithm is then given by the following steps: First, the chain is initialized at some parameter value β₀, which can chosen as a good estimate to the parameter β is as some random value. Secondly, a random sample is drawn from the proposal distribution, as β_j + 1 ← T(β_j + 1 ∣ β_j), and define ω_j = {a_AME, a_s, β_j, β_d}, while keeping all other parameters fixed. In order to decide if the sample should be accepted or rejected, the Metropolis acceptance probability q is computed,

$$\begin{array}{c}\hfill q=\frac{\mathcal{L}({\omega }_{j+1})}{\mathcal{L}({\omega }_j)}\frac{P({\omega }_{j+1})}{P({\omega }_j)},\end{array}$$(18)

where ℒ(ω) is given by Eq. (17) and P(ω) denotes a Gaussian prior. Once q is calculated, a random number, η, is drawn from a uniform distribution U[0, 1], and the proposed ω_j + 1 is accepted if η < q. If the sample is rejected, set ω_j + 1 = ω_j. This procedure is repeated until convergence.

An intuitive observation regarding the Metropolis method is that all samples with a higher posterior probability are accepted, while some samples with a lower posterior probability are accepted. By accepting with a probability given by the ratio of the sample posteriors, the density of samples within a given parameter volume is proportional to the underlying probability density. In this work the proposal probability distribution is Gaussian with a width σ determined by the standard deviation of the Gaussian prior distribution. An ensemble of Gaussian priors is used in this analysis to explore the effect of the prior on the final results.

5. Results

We are now ready to present the results derived by applying the Gibbs sampler described in Sect. 4 to the data summarized in Sect. 2. Unless otherwise noted, the default synchrotron spectral index prior is β_s = −3.1 ± 0.1. For each chain, 500 full Gibbs samples are generated. Between each main sample, we recall that the LFI frequency maps are replaced by different BEYONDPLANCK samples to marginalize over systematic effects, as discussed in Sect. 2.

5.1. Goodness-of-fit

Before examining the details of each emission component, we consider the overall goodness-of-fit of our data model. Figure 1 shows mean residual maps for each frequency map on the form r_ν = d_ν − T_νa − s_d, ν for both LFI and WMAP. The gray region indicates the analysis mask discussed in Sect. 2. Generally we observe an excellent fit in terms of Galactic signal suppression; the Galactic plane is only barely visible at the level of a few microkelvin in the LFI 30 GHz and WMAP K-band. However, at higher latitudes there are clear large-scale residuals at the level of ≲3 μK for LFI and ≲5 μK for WMAP with a clear instrumental origin; for a detailed discussion of these structures, we refer the interested reader to Jarosik et al. (2007), BeyondPlanck Collaboration (2023), Svalheim et al. (2023b), and Watts et al. (2023).

Fig. 1. Posterior mean residual maps (rν = dν − sν) in Stokes Q and U for (from top to bottom) the BEYONDPLANCK LFI 30, 44, and 70 GHz maps and the WMAP K, Ka, Q, and V bands. All maps are smoothed to a common resolution of 3° FWHM.

It is worth noting that the LFI residual maps shown in Fig. 1 appear significantly different from the corresponding maps in the main BEYONDPLANCK analysis. The reason for this is, as discussed earlier, that the WMAP K-band channel is excluded from the latter analysis. And from this plot we can intuitively understand why: The high signal-to-noise of the WMAP K-band heavily influences the fit of the Planck 30 GHz band, and potential uncontrolled systematic effects in the K-band may therefore compromise the LFI data themselves. We also note that the 30 GHz and K-band residuals are morphologically similar, but with opposite signs, indicating good agreement in terms of the Galactic signal, but with significantly different large-scale scanning-induced systematics. For the main BEYONDPLANCK analysis, which focuses on LFI, it is therefore clear why the K-band is omitted–while in the current paper, we want to establish a joint estimate between LFI and WMAP, and we want to maximize our signal-to-noise ratio, and we therefore include all data.

Figure 2 shows a corresponding χ² map, which has been normalized by the number of bands. At this level, we observe very little discernible structure, except for small excesses very near the mask, indicating that the adopted mask indeed performs well.

Fig. 2. Mean χ2/Nband of all Gibbs iterations, averaged over Stokes Q and U.

5.2. Synchrotron component

Next we examine the derived synchrotron posterior component maps. Figure 3 shows the posterior mean map for the synchrotron spectral index, and we recall that the prior in this case is β_s = −3.1 ± 0.1. Overall we see that β_s is very close to the prior mean over large portions of the sky, though we can clearly see the impact from the likelihood in some high signal-to-noise regions of the sky. In particular, we note significant deviations from the prior along the Northern Galactic Spur, in the Fan Region, and along the Galactic Plane, in agreement with Svalheim et al. (2023b). However, by far the most important conclusion to be drawn from this figure is the fact that even the combination of LFI and WMAP has very little constraining power compared to the prior with respect to the spectral index of synchrotron emission, and all results that depends on this will necessarily be sensitive to the assumed prior. In the following, we therefore consider the synchrotron index prior mean as a free hyperparameter, and all main results will be plotted as a function of this parameter.

Fig. 3. Posterior mean polarized synchrotron spectral index map.

Next, the synchrotron Stokes Q and U amplitude components are summarized in the form of posterior mean (top panels) and standard deviation (bottom panels) maps in Fig. 4. Figure 5 shows difference maps between the synchrotron amplitude map derived in this paper and those generated by the BEYONDPLANCK (top row) and Planck DR4 (bottom row) pipelines. We find excellent agreement with the different pipelines to around the 3 μK_RJ level. The most dominant differences are systematic in nature. In particular, we note that the main difference between the BEYONDPLANCK and current analysis is the addition of WMAP K-band, while the difference between the current analysis and the Planck DR4 additionally highlights the difference between the end-to-end Bayesian global approach (BeyondPlanck Collaboration 2023) and the classic frequentist single-channel approach (Planck Collaboration Int. LVII 2020).

Fig. 4. Posterior mean (top) and standard deviation (bottom) maps of the synchrotron amplitude at 30 GHz in Stokes Q (left) and U (right).

Fig. 5. Difference between the synchrotron solution derived here and the solution derived in other works. Top: Difference between the BEYONDPLANCK synchrotron solution (Svalheim et al. 2023b) and the synchrotron component derived here. Bottom: Difference between Planck DR4 and the synchrotron component derived here. We note a difference along the Galactic plane, which is likely due to the strong degeneracy between dust and synchrotron emissions in this region. The left and right columns correspond to Stokes Q and U, smoothed with a 3° FWHM Gaussian beam.

5.3. Constraints on AME polarization

Finally, we are ready to present the main results of the current paper, namely constraints on polarized AME emission, and we start the linear amplitude parameters, a_AME. As an example, these are shown in Fig. 6 for a synchrotron prior of β_s = −3.1 ± 0.1. For convenience, we normalize a_AME with respect to the 353 GHz channel, and we plot the sum of the thermal dust emission and AME, with different AME polarization fractions marked as colored curves, ranging between 0 and 10%.

Fig. 6. Fitted AME amplitudes for a synchrotron prior of βs ∼ N(−3.1, 0.1), with various models for different levels of $p_{\mathrm{AME}}^{\max }$ overplotted as colored curves. The red points are taken from Planck Collaboration XI (2020).

As discussed in Sect. 3.2, AME amplitudes are not fit at LFI 70 GHz and WMAP V-band as these bands are assumed to be fully described by synchrotron and thermal dust emission. However, Planck Collaboration XI (2020) performed a similar thermal dust emission fit for LFI 70 and 44 GHz, and found good agreement with the SED at 70 GHz. Their 44 GHz result fits shows a slight excess, while in this work we find that the template amplitude fit for BEYONDPLANCK 44 GHz (and WMAP Q-band) is also consistent with the thermal dust SED alone, at least for a synchrotron prior of β_s = −3.1 ± 0.1. Also for the lower frequency bands we find results consistent with nonpolarized AME, but with large uncertainties.

To constrain the corresponding AME polarization fraction, we map out the corresponding posterior distribution defined by the change-of-variables in Eq. (10). We create a grid of potential polarized AME SEDs for different polarization fractions, s_ν(p_AME), and calculate a log-likelihood given by

$$\begin{array}{c}\hfill -2ln\mathcal{L}(p_{\mathrm{AME}})={\displaystyle \sum _{\nu }{\left(\frac{{({\widehat{a}}_{\nu }-s_{\nu }(p_{\mathrm{AME}}))}^2}{{\sigma }_{a,\nu }^2}\right)}^2.}\end{array}$$(19)

Uncertainties in the thermal dust SED are taken into account by averaging over random draws from the MBB parameters as described in Sect. 3. Slices of the full posterior distribution for each of the β_s priors are shown in Fig. 7, showing the 95% confidence interval for $p_{\mathrm{AME}}^{\max }$, as well as the posterior peak. The nominal case (shown in the form of SED constraints in Fig. 6) is shown in Fig. 8 in orange. In this particular case, we find that the posterior peaks below 0% with an upper limit of $p_{\mathrm{AME}}^{\max }<1.0\%$ at the 95% confidence level.

Fig. 7. 95% confidence intervals for $p_{\mathrm{AME}}^{\max }$ as a function of the βs prior mean.

However, as shown in Sect. 5.2, the current data set has very little constraining power with respect to β_s. We therefore repeat the above calculations for a grid of β_s, allowing the prior mean to range between β = −3.3 and β = −3.0. The results from these analyses are summarized in Fig. 7. Here we see that for any β_s ≲ −3.1, the distributions are consistent with a vanishing AME polarization fraction, while for flatter indices a nominal positive detection emerges. The most likely explanation for this is that synchrotron emission is oversubtracted by the flat spectral index, and this may be countered by adding a spurious polarized AME component.

Figure 7 summarizes this in the form of the 95% confidence region as a function of the β_s prior mean. Here we find that the typical upper limit is $p_{\mathrm{AME}}^{\max }\lesssim 1.0\%$ at 95% confidence for β_s ≲ −3.1, while for β_s = −3.0, there is a nominal detection of polarized AME with p_AME = 2.5 ± 1.0% at 95% confidence.

5.4. Impact of CMB fluctuations

To determine the validity of our assumption to ignore the CMB component in Sect. 3.2, we run our analysis on the same set of data, but with an additional CMB component added to each frequency map. This allows us to estimate the affect the CMB component has on our $p_{\mathrm{AME}}^{\max }$ constraints.

The CMB component is generated using healpy’s (Zonca et al. 2019) synfast routine, using the Planck best-fit ΛCDM power spectrum⁴ (Planck Collaboration V 2020). The results are shown in Fig. 8, where we show the $p_{\mathrm{AME}}^{\max }$ posterior slices of the fiducial analysis in orange and the fiducial-plus-CMB analysis in cyan. We see that an additional set of CMB fluctuations increases the upper limits on $p_{\mathrm{AME}}^{\max }$ by 10 − 15%. As a result, the main limits quoted in paper are likely overestimated, and therefore conservative, as the addition of unmodeled CMB fluctuations lead to slightly weaker constraints.

Fig. 8. Comparison of likelihood slices for $p_{\mathrm{AME}}^{\max }$ for a synchrotron prior of βs ∼ N(−3.1, 0.1) with (blue) and without (orange) additional CMB fluctuations.

6. Summary

We have presented the first constraints on large-scale polarized AME derived based on the combination of Planck LFI and WMAP polarization measurements. These constraints are derived within the context of the BEYONDPLANCK framework, which allows for detailed propagation of systematic errors; the current work is an explicit example of how the BEYONDPLANCK Monte Carlo frequency map ensembles may be used for error propagation in external analyses. However, we do note that only Planck LFI data are currently processed through this machinery, while traditionally processed maps are employed for WMAP. The effect of this mismatch is seen in various residual and difference maps, clearly demonstrating the presence of instrumental systematic effects in one or both experiments. Work is currently on-going to reprocess the WMAP data jointly with LFI (Watts et al. 2023), but for now we note that the two experiments agree very well in terms of Galactic estimates, even if their instrumental residuals differ. We also note that the χ² of our fits is excellent, which indicates that systematic errors are small compared to the noise level of each experiment.

We find that the LFI and WMAP data only have limited constraining power with respect to the spectral index of polarized synchrotron emission, and any constraint on the AME polarization fraction is significantly affected by the choice of synchrotron spectral index prior: The current data set is simply not able to robustly and simultaneously constrain polarized AME and synchrotron emission. As a result, we choose to report a limit on the AME polarization fraction that is explicitly β_s-dependent. Specifically, for β_s ≲ −3.1, we find an upper limit of $p_{\mathrm{AME}}^{\max }\lesssim 0.6\%$ at 95% confidence, while for β_s ≳ −3.1 there is formally a positive detection of polarized AME power.

The maximum polarization fraction of 3.5% (for β_s = −3.0) derived here is in good agreement with and more stringent than the large-scale Macellari et al. (2011) limit of p_AME < 5%. Our upper limit also agrees with more stringent upper limits placed in individual molecular cloud complexes, such as those from Battistelli et al. (2006), Génova-Santos et al. (2017), Poidevin et al. (2019), though we emphasize that the AME may be significantly less polarized than the theoretical maximum in these regions based on the magnetic field geometry and level of coherence within the beam and along the line of sight. The nondetection of polarization is in agreement with spinning dust emission theory, driven either by PAHs or nanosilicate grains, but we do not have the sensitivity to discriminate between models with modest levels of polarization and models which predict no polarization. AME driven by emission mechanisms which have polarization fractions of greater than a few percent, such as spinning iron grains (Hoang et al. 2016), are excluded.

We do note that both Planck 2018 (Planck Collaboration V 2020) and BEYONDPLANCK (BeyondPlanck Collaboration 2023; Svalheim et al. 2023b) prefer a steep spectral index of β_s ≈ −3.3 at high Galactic latitudes. However, those analyses consider only LFI data, and are therefore susceptible to even stronger degeneracies than what is observed in the current analysis, which also includes WMAP data. It is conceivable that the steep spectral index they observe is partially caused by polarized AME. Another complication is the possibility of curvature in the synchrotron SED, which has not been addressed in the current analysis.

The only way to actually break these degeneracies is through additional high signal-to-noise observations, both by ground-based low-frequency experiments such as C-BASS (Jew et al. 2019), S-PASS (Carretti et al. 2019), QUIJOTE (Génova-Santos et al. 2015b), and by next-generation large-scale experiments such as Simons Observatory (Ade et al. 2019) and CMB-S4 (Abazajian et al. 2019). Indeed, we would argue that high-sensitivity characterization of the low-frequency polarized foreground should be a primary objective for near-term CMB experiments, as the complexity of this frequency range can have dramatic consequences for next-generation CMB B-mode satellite experiments.

This is illustrated in Fig. 9, which provides an overview of the polarized microwave sky from 1 to 500 GHz. Green and red bands show synchrotron and thermal dust emission as constrained by BEYONDPLANCK (Svalheim et al. 2023b), while the orange curve shows the conservative upper limit on polarized AME for a synchrotron spectral index prior of β_s = −3.1 ± 0.1, namely $p_{\mathrm{AME}}^{\max }\lesssim 1.0\%$. Clearly, the presence of a polarized AME component with an amplitude this level would be critically important for any future B-mode experiments. In sum, the current analysis serves as a useful reminder about how little is still known about low-frequency polarized foregrounds, even after Planck and WMAP, and more data are desperately needed.

Fig. 9. Polarization spectrum in the microwave frequency range. In addition to the traditional spectrum which only includes thermal dust and synchrotron emission, an AME (spinning dust) component is added here with a maximum polarization fraction of 1%. The thermal dust and synchrotron spectra are determined by results from BEYONDPLANCK LFI analysis (Svalheim et al. 2023b). The cyan and orange vertical bands correspond to the WMAP, and Planck bands respectively.

Acknowledgments

We thank Prof. Pedro Ferreira and Dr. Charles Lawrence for useful suggestions, comments and discussions. We also 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 776282 (COMPET-4; BEYONDPLANCK), 772253 (ERC; BITS2COSMOLOGY), and 819478 (ERC; COSMOGLOBE). In addition, the collaboration acknowledges support from ESA; ASI and INAF (Italy); NASA and DoE (USA); Tekes, Academy of Finland (grant no. 295113), CSC, and Magnus Ehrnrooth foundation (Finland); RCN (Norway; grant nos. 263011, 274990); and PRACE (EU).

References

All Figures

	Fig. 1. Posterior mean residual maps (rν = dν − sν) in Stokes Q and U for (from top to bottom) the BEYONDPLANCK LFI 30, 44, and 70 GHz maps and the WMAP K, Ka, Q, and V bands. All maps are smoothed to a common resolution of 3° FWHM.
In the text

	Fig. 2. Mean χ2/Nband of all Gibbs iterations, averaged over Stokes Q and U.
In the text

	Fig. 3. Posterior mean polarized synchrotron spectral index map.
In the text

	Fig. 4. Posterior mean (top) and standard deviation (bottom) maps of the synchrotron amplitude at 30 GHz in Stokes Q (left) and U (right).
In the text

Fig. 5. Difference between the synchrotron solution derived here and the solution derived in other works. Top: Difference between the BEYONDPLANCK synchrotron solution (Svalheim et al. 2023b) and the synchrotron component derived here. Bottom: Difference between Planck DR4 and the synchrotron component derived here. We note a difference along the Galactic plane, which is likely due to the strong degeneracy between dust and synchrotron emissions in this region. The left and right columns correspond to Stokes Q and U, smoothed with a 3° FWHM Gaussian beam.

In the text

	Fig. 6. Fitted AME amplitudes for a synchrotron prior of βs ∼ N(−3.1, 0.1), with various models for different levels of $p_{\mathrm{AME}}^{\max }$ overplotted as colored curves. The red points are taken from Planck Collaboration XI (2020).
In the text

	Fig. 7. 95% confidence intervals for $p_{\mathrm{AME}}^{\max }$ as a function of the βs prior mean.
In the text

	Fig. 8. Comparison of likelihood slices for $p_{\mathrm{AME}}^{\max }$ for a synchrotron prior of βs ∼ N(−3.1, 0.1) with (blue) and without (orange) additional CMB fluctuations.
In the text

Fig. 9. Polarization spectrum in the microwave frequency range. In addition to the traditional spectrum which only includes thermal dust and synchrotron emission, an AME (spinning dust) component is added here with a maximum polarization fraction of 1%. The thermal dust and synchrotron spectra are determined by results from BEYONDPLANCK LFI analysis (Svalheim et al. 2023b). The cyan and orange vertical bands correspond to the WMAP, and Planck bands respectively.

In the text