DR1 results. I. Improved Wilkinson Microwave Anisotropy Probe maps through Bayesian end-to-end analysis

Watts et al.

We present  Data Release 1, which implements the first joint analysis of WMAP and Planck LFI time-ordered data, processed within a single Bayesian end-to-end framework. This framework builds directly on a similar analysis of the LFI measurements by the  collaboration, and approaches the CMB analysis challenge through Gibbs sampling of a global posterior distribution, simultaneously accounting for calibration, mapmaking, and component separation. The computational cost of producing one complete +LFI Gibbs sample is 812 CPU-hr, of which 603 CPU-hrs are spent on  low-level processing; this demonstrates that end-to-end Bayesian analysis of the  data is computationally feasible. We find that our  posterior mean temperature sky maps and CMB temperature power spectrum are largely consistent with the official WMAP9 results. Perhaps the most notable difference is that our CMB dipole amplitude is $3366.2\pm1.4\,\mathrm{\mu K}$, which is 11 higher than the  estimate and $2.5\,\sigma$ higher than ; however, it is in perfect agreement with the HFI-dominated  PR4 result. In contrast, our  polarization maps differ more notably from the  results, and in general exhibit significantly lower large-scale residuals. We attribute this to a better constrained gain and transmission imbalance model. It is particularly noteworthy that the -band polarization sky map, which was excluded from the official  cosmological analysis, for the first time appears visually consistent with the -band sky map. Similarly, the long standing discrepancy between the  -band and LFI 30 GHz maps is finally resolved, and the difference between the two maps appears consistent with instrumental noise at high Galactic latitudes. Relatedly, these updated maps allow us for the first time to combine  and LFI polarization data into a single coherent model of large-scale polarized synchrotron emission. Still, we identify a few issues that require additional work, including 1) low-level noise modeling, 2) large-scale temperature residuals at the 1–2 level; and 3) a strong degeneracy between the absolute $K$-band calibration and the dipole of the anomalous microwave emission component. We conclude that leveraging the complementary strengths of  and LFI has allowed the mitigation of both experiments' weaknesses, and resulted in new state-of-the-art  sky maps. All maps and the associated code are made publicly available through the  web page.

§ INTRODUCTION

The discovery of the cosmic microwave background (CMB) by <cit.> marked a paradigm shift in the field of cosmology, providing direct evidence that the Universe was once much hotter than it is today, effectively ruling out the steady-state theory of the universe <cit.>. This discovery spurred a series of ground-breaking cosmological experiments, including the Nobel Prize-winning measurements by -FIRAS that confirmed the blackbody nature of the CMB <cit.> and -DMR that measured temperature variations from the primordial gravitational field <cit.>.

The NASA-funded Wilkinson Microwave Anisotropy Probe (; ) mission was launched a decade after -DMR, and mapped the microwave sky with 45 times higher sensitivity and 33 times higher angular resolution, thereby revolutionizing our understanding of early universe physics <cit.>.  In addition, the 3-year measurements presented by <cit.> included the first ever detection of large-scale polarization in the CMB, opening a new window into the process of cosmic reionization. As quantified by <cit.>, the permissible parameter space volume for a standard  model was decreased by a factor of 68,000 by , and the best pre- determination of the age of the universe was $t_0<14\,\mathrm{Gyr}$ from Boomerang <cit.>, with best-fit values of 9–11 Gyr; the latter values in apparent contradiction with direct measurements of the oldest globular clusters <cit.>.

The ESA-led  satellite <cit.> was developed concurrently with , and their operation lifetimes partially overlapped, with  observing from 2001–2011 and  from 2009–2013. 's stated goal was to fully characterize the primary CMB temperature fluctuations from recombination, as well as to characterize the polarized microwave sky on large angular scales.  Overall, 's raw CMB sensitivity was an order of magnitude higher than 's, and its angular resolution more than twice higher. Today,  represents the state-of-the-art in terms of full-sky microwave sky measurements.

comprised two independent experiments, namely the Low Frequency Instrument (LFI; ) and High Frequency Instrument (HFI; ), respectively. The LFI detectors were based on HEMT (high electron mobility transistor) amplifiers, spanning three frequency channels between 30 and 70 GHz, while the HFI detectors were based on spider-web and polarization sensitive bolometers, and spanned six frequency channels between 100 and 857 GHz. For comparison,  was also HEMT-based, with comparable sensitivity to LFI alone, and spanned five frequencies between 23 and 94 GHz. At the same time, the two experiments implemented very different scanning strategies, and as a result they are highly complementary and synergistic; together they provide a clearer view of the low-frequency microwave sky than either can alone.

Towards the end of the  analysis phase it became clear that the interplay between instrument calibration and astrophysical component separation was a main limiting factor in terms of systematic effects for high signal-to-noise measurements <cit.>. Specifically, in order to calibrate the instrument to sufficient precision, it became apparent that it was necessary to know the true sky to a comparably high precision – but to know the sky, it was also necessary to know the instrumental calibration. The data analysis was thus fundamentally circular and global in nature. The final official  LFI analysis performed four complete iterations between calibration and component separation <cit.>, aiming to probe this degeneracy. However, it was recognized that this was not sufficient to reach full convergence, and this sub-optimality led to the  project <cit.>, which aimed to perform thousands of complete analysis cycles, as opposed to just a handful. This framework was implemented using the  <cit.> code, a CMB Gibbs sampler that performs integrated high-level and low-level parameter estimation in a single integrated framework. This analysis demonstrated the feasibility of end-to-end CMB analysis through Gibbs sampling analysis, while at the same time it provided the highest-quality LFI maps to date.

Rather than simply probing the degeneracy between instrument calibration and component separation, a better solution is to actually break it. The optimal approach to do so is by jointly analyzing complementary datasets, each of which provide key information regarding the full system. This insight led to the [<https://cosmoglobe.uio.no>] initiative, which is an Open Source and community-wide effort that aims to derive a single joint model of the radio, microwave, and sub-millimeter sky by combining all available state-of-the-art experiments. An obvious first extension of the LFI-oriented  project is to analyze the  measurements in the same framework. Indeed, already as part of the  suite of papers, <cit.> integrated  -band time-ordered data (TOD) into the  framework, calibrated off of the  sky model.

In this paper, we present the first end-to-end Bayesian analysis of the full  TOD, processed within the  framework. As such, this paper also presents the first ever joint analysis of two major CMB experiments (LFI and ) at the lowest possible level, and it therefore constitutes a major milestone of the  initiative. We refer to the current products as  Data Release 1 (CG1), and the scientific results from this are described in a series of four papers. The current paper gives a detailed discussion of data processing methods, instrumental parameters, frequency maps, and preliminary astrophysical results, while updated constraints on anomalous microwave emission and polarized synchrotron emission are presented by <cit.> and <cit.>, respectively. <cit.> use these new products to provide new constraints on cosmic birefringence. In the future, many more datasets and astrophysical components will be added to this framework, gradually providing stronger and stronger constraints on both the true astrophysical sky and the instrumental calibration of all previous experiments.

The rest of this paper is organized as follows. In Sect. <ref>, we provide a brief review of the Bayesian end-to-end statistical framework used in this work, before describing the underlying data and computational expenses in Sect. <ref>. The main results, as expressed by the global posterior distribution, are described in Sects. <ref>–<ref>, summarizing instrumental parameters, frequency sky maps, and preliminary astrophysical results, respectively. In Sect. <ref> we address unresolved issues that should be further analyzed in future work. We conclude in Sect. <ref>, and lay a path forward for the  project.

§ END-TO-END BAYESIAN CMB ANALYSIS

The general computational analysis framework used in this work has been described in detail by <cit.> and <cit.> and references therein. In this section, we give a brief summary of the main points, and emphasize in particular the differences with respect to earlier work.

§.§ Official  instrument model and analysis pipeline

The main goal of the current paper is to perform a similar analysis as <cit.> did for  LFI, but this time including the  in terms of time-ordered data, and thereby solve some of the long-standing unresolved issues with the official maps, in particular related to poorly constrained large-scale polarization modes. Before presenting our algorithm, however, it is useful to briefly review the official  instrument model and analysis pipeline, which improved gradually over a total of five data releases, often referred to as the 1-, 3-, 5-, 7-, and 9-year data releases, respectively. Unless otherwise noted, we will refer to the final 9-year results <cit.>, and denote these as . A concise summary of the  mission, data processing, and results is available in <cit.>. The full data archive can be found on LAMBDA.[<https://lambda.gsfc.nasa.gov/product/wmap/dr5/m_products.html>]

The  satellite carried twenty differential polarization-sensitive
radiometers, grouped into ten differencing assemblies (DAs), where one
was sensitive to the difference in signal at one polarization
orientation and the other sensitive to the orthogonal
polarization. In total, of the ten DAs there were:
one -band (23 GHz), one -band (33 GHz), two -bands (41 GHz),
two -bands (61 GHz), and four -bands (94 GHz). Each radiometer comprised two detector diodes, which each recorded a science sample every $1.536/N_\mathrm{obs}$ seconds, where $N_\mathrm{obs}$ is 12, 12, 15, 20 and 30 for , , , , and , respectively. The raw data are recorded as 16-bit integers with units du (digital unit).

The  bandpasses were measured pre-launch on the ground, sweeping a signal source through 201 frequencies and recording the output <cit.>. The bandpass responses available on LAMBDA have not been updated since the initial data release. However, as noted by <cit.>, there has been an observed drift in the center frequency of , , , and -band corresponding to a $\sim$$\,0.1\,\%$ decrease over time. In practice, this did not affect the  data processing because each year was mapped separately and co-added afterwards. An effective frequency calculator was delivered in the DR5 release as part of the IDL library to mitigate this effect during astrophysical analyses.[<https://lambda.gsfc.nasa.gov/product/wmap/dr5/m_sw.html>]

The beams were characterized in the form of maps, with separate products for the central portion of the beam pattern and the far sidelobes. The main beam and near sidelobes were characterized using a combination of physical optics codes and observations of Jupiter for each horn separately. The maps of Jupiter were then combined with the best-fit parameters from physical optics codes to create a map of the beam response <cit.>.

Far sidelobes were estimated using a combination of laboratory measurements and Moon data taken during the mission <cit.>, as well as a physical optics model described by <cit.>. To remove the far sidelobe  in the TOD, an estimate was calculated by convolving the intensity map and the orbital dipole signal with the measured sidelobe signal <cit.>. Although the sidelobe pickup was modeled by <cit.>, it was determined that the results were small enough to be neglected and have not been explicitly reported in any of the subsequent  data releases.

The  pointing solution was determined using the boresight vectors of individual feedhorns in spacecraft coordinates, in combination with on-board star trackers. Thermal flexure of the tracking structure introduced small pointing errors, as discussed by <cit.>. Using the temperature variation measured by onboard thermistors, the pointing solution was corrected using a model that returns angular deviation per kelvin. The residual pointing errors were computed using observations of Jupiter and Saturn, and the reported upper limit was estimated to be 10 <cit.>.

The  data were calibrated by jointly estimating the time-dependent gains, $\g$, and baselines, $\bv$, as described by <cit.>, <cit.>, and <cit.>.
The TOD were initially modeled as having
constant gain and baseline for a 1–24 hour period, with parameters that were fit to the orbital
dipole assuming $T_0$ from <cit.> and a map made from a previous
iteration of the mapmaking procedure. Once the gain and baseline solution had
converged, the data were fit to a parametric form of the radiometer response
as a function of housekeeping data, given in Appendix A of <cit.>.

had two primary mirrors positioned on opposite sides of the vertical satellite axis, tilted approximately $19.5^\circ$ towards the Solar shield. Essentially, when horn A was pointed at pixel $p_\A$, horn B was pointed at a pixel $p_\B$ approximately $141^\circ$ away <cit.>. The incoming radiation was differenced in the electronics before being deposited on the detectors, recording radiation proportional to the observed maps $\boldsymbol m$ at their respective pixels, $\boldsymbol m_{p_\A}-\boldsymbol m_{p_\B}$ and $\boldsymbol m_{p_\B}-\boldsymbol m_{p_\A}$ <cit.>. Each radiometer had a partner that observed the same pixels with sensitivity to the orthogonal polarization direction. Taking all these effects into account, the total data model for a single radiometer is given by

\begin{align}
d_t^\mathrm{imbal}&\propto (1+x_\mathrm{im})T_{p_\A}-(1-x_\mathrm{im})T_{p_\B}
\\
\end{align}

where $T_{p_\A}$ and $T_{p_\B}$ are the A- and B-side antenna temperatures, and $x_\mathrm{im}$ is the differential optical pickup between horns A and B. This effect is taken into account during mapmaking. However, inaccuracies in the determination of $x_\mathrm{im}$ will yield a spurious polarization component, and create artificial imbalance modes due to coupling with the sky signal, in particular with the bright Solar CMB dipole <cit.>. The  transmission imbalance factors were fit to the Solar dipole in TOD space, accounting for both common and differential modes <cit.>.

Data were flagged and masked before the final mapmaking step. In particular, station-keeping maneuvers, solar flares, and unscheduled events caused certain data to be unusable – the full catalog of these events is listed in Table 1.8 of <cit.>. In addition, data were masked depending on the channel frequency and the planet itself, with the full list of exclusion radii enumerated in Table 4 of <cit.>.

To create the sky maps $\boldsymbol m$, the calibrated data were put into the asymmetric mapmaking equation,

\begin{equation}
\label{eq:asymmetric_mapmaking}
\mathsf P_\mathrm{am}^T\mathsf N^{-1}\mathsf P\boldsymbol m=\mathsf P^T_\mathrm{am}\mathsf N^{-1}\boldsymbol d,
\end{equation}

where $\mathsf N$ is the noise covariance matrix, and the pointing matrix $\mathsf P$ is implicitly defined for each datastream, $\boldsymbol d_1$ and $\boldsymbol d_2$ sensitive to different polarization orientations.
The asymmetric mapmaking matrix, $\mathsf P_\mathrm{am}$, was used because, as noted by <cit.>, large signals observed in one beam could leak into the solution for the pixel observed by the other beam, leading to incorrect signals in the final map. The asymmetric mapmaking solution is defined by only updating the matrix multiplication for beam A when beam A is in a high emission region and beam B is not, and vice versa. <cit.> also identified that these effects are pronounced when one horn is crossing a large temperature gradient, leading to excesses $140^\circ$ away from the Galactic center if an appropriate processing mask is not used.
For each side $\A/\B$, the maps are defined as a function of the Stokes parameters $T_{\A/\B}$, $Q_{\A/\B}$, and $U_{\A/\B}$, with polarization angle $\gamma_{\A/\B}$, such that

\begin{align}
\boldsymbol d_1&=\mathsf P_1\boldsymbol m
\nonumber
\\
\nonumber
\\
\label{eq:mapeq1}
\end{align}

\begin{align}
\boldsymbol d_2&=\mathsf P_2\boldsymbol m
\nonumber
\\
\nonumber
\\
\label{eq:mapeq2}
\end{align}

In this formalism, $S_{\A/\B}$ acts as an extra Stokes parameter that absorbs the effects of differing bandpass responses between radiometers $\boldsymbol d_1$ and $\boldsymbol d_2$ <cit.>.

An accurate noise model was necessary both to perform the maximum likelihood mapmaking and for the evaluation of the dense time-space inverse noise covariance matrix  $\mathsf N^{-1}$. The  team defined this in the form of a time domain autocorrelation function that was estimated separately for each year of data. This was then Fourier transformed, inverted, and inverse Fourier transformed to create an effective inverse noise operator $N_{tt'}^{-1}$.
Finally, to create the sky maps themselves, the  team processed the data one year at a time, producing maps by solving Eq. (<ref>) using the iterative Bi-Conjugate Gradient Stabilized Method <cit.>.

§.§  instrument model

Time-ordered data segment for the 113 diode. From top to bottom, the panels show 1) raw uncalibrated TOD $\boldsymbol d$; 2) sky signal $\boldsymbol s_\mathrm{sky}$; 3) calibrated correlated noise $\boldsymbol n_\mathrm{corr}$; 4) orbital CMB dipole signal $\boldsymbol s_\mathrm{orb}$; 5) sidelobe correction $\boldsymbol s_\mathrm{sl}$; 6) bandpass leakage correction $\boldsymbol s_\mathrm{leak}$; and 7) residual TOD, $\boldsymbol d_\mathrm{res}=(\boldsymbol d -\boldsymbol n_\mathrm{corr}-\boldsymbol b)/g-\boldsymbol s_\mathrm{sky}-\boldsymbol s_\mathrm{orb}-\boldsymbol s_\mathrm{leak} -\boldsymbol s_\mathrm{sl}$, in units of $\sigma_0[\mathrm{du}]$ for this TOD segment. Note that the vertical range varies and units vary from panel to panel.

A fundamental difference between the  and  analysis pipelines (and those of most other CMB experiments) is that while the  pipeline models each channel in isolation, the  framework simultaneously considers all data, both internally within , and also from all other sources, and most notably from  LFI. The main advantage of such a global approach is significantly reduced parameter degeneracies, as data from observations with different frequency coverages and instrumental designs break the same degeneracies. For this approach to be computationally tractable, one must establish a global parametric model that simultaneously accounts for both the astrophysical sky and all relevant instruments. For the current +LFI oriented analysis, we adopt the following expression <cit.>,

\begin{equation}
\label{eq:model}
\boldsymbol d =\mathsf G\mathsf P[\mathsf B^\mathrm{symm}\mathsf M\boldsymbol a+\mathsf B^\mathrm{4\pi}(\boldsymbol s^\mathrm{orb}
+\boldsymbol s^\mathrm{fsl})] + \boldsymbol s^\mathrm{inst}+ \boldsymbol n^\mathrm{corr}+\boldsymbol n^\mathrm w,
\end{equation}

where $\mathsf G$ is the time-dependent gain in the form of the matrix $\mathrm{diag}(g_t)$; $\mathsf P$ is the $n_p\times n_\mathrm{TOD}$ pointing matrix, where $n_p$ is the number of pixels and $n_\mathrm{TOD}$ length of the TOD;
$\mathsf B^\mathrm{symm}$ and $\mathsf B^{4\pi}$ are the symmetrized and full asymmetric beam, respectively; $\mathsf M$ is the mixing matrix between a given sky component $c$ with spectral energy distribution $f_c(\nu/\nu_{0,c})$ and a detector $i$ with bandpass $\tau_i(\nu)$, given by

\begin{equation}
\mathsf M_{ci}=\int d\nu\,\tau_i(\nu)f_c(\nu/\nu_{c,0}).
\end{equation}

(In practice, $\M$ also accounts for unit conversion, but this is suppressed for readability in this expression; see  for further details.)
The maps $\boldsymbol a$ represent the Stokes parameters for each astrophysical component, while $\boldsymbol s^\mathrm{orb}$ is the orbital dipole induced by the motion of the telescope with respect to the Sun, and $\boldsymbol s^\mathrm{fsl}$ is the time-dependent far sidelobe signal. Following <cit.>, we model the correlated noise component $\boldsymbol n^\mathrm{corr}$ in terms of a $1/f$ power spectral density (PSD), which explicitly takes the form $P_n(f) = \sigma^2_0 (1 + (f/f_\mathrm{knee})^\alpha)$, where $\sigma_0$ denotes the white noise amplitude, $f_\mathrm{knee}$ is the so-called $1/f$ knee frequency, and $\alpha$ is a free power law slope. For notational purposes, we denote the set of all correlated noise parameters by $\boldsymbol\xi_n = \{\sigma_0, f_{\mathrm{knee}}, \alpha\}$. We note that this model represents a significant approximation compared to the more flexible  autocorrelation model, as the actual  noise is known to be colored  at high temporal frequencies <cit.>. The main impact of this approximation is a worse-than-expected $\chi^2$ goodness of fit statistic. However, measured in absolute noise levels the effect is very small, and has very little if any impact on the final science results; for further discussion of this approximation, see Sect. <ref>.

The term $\boldsymbol s^\mathrm{inst}$ denotes any instrument-specific terms that might be required for a given experiment. For instance, for LFI it is used to model the 1 Hz spike contribution due to electronic cross-talk. For , we use it for first-order baseline corrections, and set $s^\mathit{WMAP}_t = b_0 + b_1\,t$, where $b_0$ and $b_1$ represent the mean and slope of the baselines over the data segment in question. We note that while the  team fitted a single constant baseline over either 1- or 24-hour periods, our data segments are typically several days long (corresponding to a number of samples chosen to optimize Fourier transforms). A natural question is therefore whether nonlinear baseline variations could induce artifacts. In this regard, it is important to note that the correlated noise component effectively acts as a single-sample baseline correction that can absorb by far most such nonlinearities, as long as their total effect on the power spectrum does not exceed that imposed by the $1/f$ model. In practice, this is a very mild constraint. At the same time, visual inspection of $\n^{\mathrm{corr}}$ projected into sky maps provides a very powerful check on any potential baseline residuals, which will appear as correlated stripes aligned with the  scanning path; for the full set of correlated noise maps derived for all ten  DAs, see Fig. <ref> in Appendix <ref>. Such maps have been used to identify and mitigate modeling errors several times in the course of this analysis. In sum, it is important to note that the  model allows for a more flexible baseline behaviour than the  pipeline, even though the dedicated baseline parameters themselves apply to relatively long timescales.

A third notable difference between the  and  data models concerns bandpass mismatch. While the  pipeline simply projects out any bandpass difference from the polarization maps by solving for the spurious $S$ maps, we model it explicitly through the use of the global astrophysical sky model <cit.>. Explicitly, the expected calibrated sky signal for diode $i$ is given by

\begin{equation}
m_{p,i}=\mathsf B_{p,p'}\sum_c\mathsf M_{c,i}a^c_{p'}+n_{i,p}^\mathrm w.
\end{equation}

Since $\mathsf M_{c,i}$ encodes the bandpass response of every detector $i$ to every sky component $c$, the detector-specific maps, $\boldsymbol m_i$, will each be slightly different depending on their bandpass $\tau_i$. Therefore, before averaging different detectors together, we estimate the average over all detectors in a given frequency channel $\boldsymbol m\equiv \langle \boldsymbol m_i\rangle$, and subtract it directly in the timestream;

\begin{equation}
\delta s_{t,i}^\mathrm{leak}=\mathsf P_{t,p}^i\mathsf B_{p,p'}^i\left(\boldsymbol m_{i,p'}-\boldsymbol m_{p'}\right).
\end{equation}

This leakage term uses the expected bandpass response to remove the expected component that deviates from the mean in the timestream, directly reducing polarization contamination.

To build intuition regarding this model, we plot in Fig. <ref> both the TOD and the individual model components for an arbitrarily selected ten-minute segment for the 's 113 diode. The uncalibrated data, $\boldsymbol d_\mathrm{raw}$, are displayed in the top panel, with the sky signal $\boldsymbol s_\mathrm{sky}=\mathsf P\mathsf B^\mathrm{symm}\mathsf M\boldsymbol a$ plotted directly underneath. The next four panels show the correlated noise realization $\boldsymbol n_\mathrm{corr}$, the orbital dipole $\boldsymbol s_\mathrm{orb}$, the far sidelobe contribution $\boldsymbol s_\mathrm{sl}$, and the bandpass leakage $\boldsymbol s_\mathrm{leak}$. Finally, we also plot the time-ordered residual for this segment of data, obtained by subtracting the model from the raw data, in units of the estimated white noise level.

§.§ Sky model

Following <cit.>, we assume that the sky across the frequency range of interest can be modeled as a linear combination of CMB fluctuations ($\vec{a}_{\mathrm{CMB}}$), synchrotron ($\vec{a}_{\mathrm{s}}$), free-free emission ($\vec{a}_{\mathrm{ff}}$), anomalous microwave emission (AME; $\vec{a}_{\mathrm{AME}}$), thermal dust ($\vec{a}_{\mathrm{d}}$), and radio point sources ($\vec{a}_{j,\mathrm{src}}$). Explicitly, we assume that the astrophysical sky (in units of brightness temperature) may be modeled as follows,

\begin{align}
\vec{s}_{\mathrm{RJ}} &= \left(\vec{a}_{\mathrm{CMB}}+\vec{a}_{\mathrm{quad}}(\nu)\right) \frac{x^2 \e^x}{(\e^x -1)^2}+\label{eq:cmb_astsky}\\
&+ \vec{a}_{\mathrm{s}} \left(\frac{\nu}{\nuzeros}\right)^{\bsynch} + \label{eq:synch_astsky}\\
&+ \vec{a}_{\mathrm{ff}} \left(\frac{\nuzeroff}{\nu}\right)^2 \frac{g_{\mathrm{ff}}(\nu;\Te) }{g_{\mathrm{ff}}(\nuzeroff;\Te)} +\label{eq:ff_astsky}\\
&+ \vec{a}_{\mathrm{AME}} \e^{\beta_{\mathrm{AME}}(\nu-\nu_{0,\mathrm{AME}})}+\label{eq:ame_astsky}  \\
&+ \vec{a}_{\mathrm{d}} \left(\frac{\nu}{\nuzerod}\right)^{\bdust+1} \frac{\e^{h\nuzerod/\kB\Tdust}-1}{\e^{h\nu/\kB\Tdust}-1}+ \label{eq:dust_astsky}\\
&+ U_{\mathrm{mJy}} \sum_{j=1}^{N_{\mathrm{src}}} \vec{a}_{j,\mathrm{src}} \left(\frac{\nu}{\nuzerosrc}\right)^{\alpha_{j,\mathrm{src}}-2}, \label{eq:sum_ptsrc}
\end{align}

where $x=h\nu/kT_{\mathrm{CMB}}$; $\nu_{0,c}$ is a reference frequency for component $c$; $\beta_{\mathrm{s}}$ is a power-law index for synchrotron emission (which may take different values for temperature and polarization); $T_e$ is the electron temperature, and $g_{\mathrm{ff}}$ is the so-called Gaunt factor <cit.>; $\beta_{\mathrm{AME}}$ is an exponential scale factor for AME emission (see below); $\beta_{\mathrm{d}}$ and $T_{\mathrm{d}}$ are the emissivity and temperature parameters for a single modified blackbody thermal dust model; $\alpha_{j,\mathrm{src}}$ is the spectral index of point source $j$ relative to the same source catalog as used by <cit.>; and $U_{\mathrm{mJy}}$ is the conversion factor between flux density (in millijansky) and brightness temperature (in $\mathrm{K_{RJ}}$) for the channel in question. Finally, $\a_{\mathrm{quad}}$ accounts for a relativistic quadrupole correction due to the Sun's motion through space <cit.>.

In general, this model is nearly identical to the one adopted by <cit.>. However, there is one notable exception, namely the spectral energy density (SED) for the AME component, $\boldsymbol s_0^\mathrm{sd}(\nu)$. In this work, we adopt a simple exponential function for this component, as for instance proposed by <cit.>, and this is notably different from the  model <cit.> that was used in the  analysis. The motivation for this modification is discussed in detail by <cit.>. First and foremost, the current combination of  and LFI data appears to prefer a higher AME amplitude at frequencies between 40 and 60 GHz than can easily be supported by . This was first noted by <cit.>, who solved this issue by introducing a second independent AME component. For the original  analysis, on the other hand, this excess was not statistically significant, simply because that analysis did not include the powerful  -band data. In the current analysis, the excess is obvious. The observation that a simple one-parameter exponential model fits the data as well as the complicated multi-parameter model of <cit.> is a novel result from the current work. Indeed, it performs about as well as the commonly used log-normal model derived by <cit.>, which also has one extra parameter. By virtue of having fewer degrees of freedom than any of the previous models, we adopt the exponential model.

§.§ Priors and poorly measured modes

Dependence of AME amplitude evaluated at 22 GHz on the absolute calibration. Each map comes from the fifth iteration of a dedicated  run that fixed $g_0$ while letting all other TOD parameters be fit. The values of $g_0=1.175$ and $g_0=1.187$ represent $6\sigma$ draws from the prior distribution with mean $1.181$ and standard deviation $0.001$. Extreme outliers were chosen to illustrate this effect. The dipole visible in the top and bottom panels is aligned perfectly with the Solar dipole, and is directly due to variations in the -band absolute calibration.

The model described in Sects. <ref> and <ref> is prone to several degeneracies, allowing for unphysical solutions to be explored in the Gibbs chain. Such unphysical degeneracies are highly undesirable for two main reasons. First, they increase the statistical uncertainties on most (if not all) other important parameters in the model – sometimes to the point that the target quantity is rendered entirely unmeasurable. Secondly, and perhaps even more importantly, the data model described above is known to be a (sometimes crude) approximation to the real observations, and there will invariably be modeling errors. Degeneracies tend to amplify their impact, in the sense that any unconstrained parameters will typically be used to fit such small modeling errors. For both these reasons, it is preferable to impose either informative or algorithmic priors on the unconstrained parameters, rather than to leave them entirely unconstrained in the model.

An important example of an algorithmic prior is the foreground smoothing prior used by <cit.> and <cit.>, which dictates that astrophysical foregrounds must be smooth on small angular scales. This is justified by noting that the angular spectrum on large and intermediate scales typically falls as a power-law in multipole space; extrapolating this into the noise dominated regime prevents the overall foreground model from becoming degenerate at small scales.

Correspondingly, important examples of informative priors are the use of HFI constraints on the thermal dust SED parameters, $\beta_{\mathrm{d}}$ and $T_{\mathrm{d}}$ in . Because that analysis only included the highest HFI frequency channel, they had very little constraining power on the thermal dust SED. Rather than trying to fit these directly from LFI  alone, they instead imposed informative Gaussian priors on each of these parameters, as derived from the HFI observations <cit.>.

Unless otherwise noted, we adopt the same algorithmic and informative priors as <cit.>. However, there are three notable exceptions, as detailed below. All of these are dictated either by the fact that we include the  -band channel (which has a strong impact on the low-frequency foreground model), or by the fact that we now process the  data in the time domain, and therefore are subject to the same degeneracies as the official  low-level pipeline; degeneracies that were solved with either implicit or explicit priors in the original analysis.

First and foremost, and discussed further in Sect. <ref>, we observe a very strong degeneracy between the absolute calibration of the -band channel and the dipole of the AME map. This is not unexpected, considering -band is by far the strongest channel in terms of AME signal-to-noise ratio, exceeding that of LFI 30 GHz by about 50 % (see Sect. <ref> for details). Effectively, a small variation in the absolute gain may be countered by subtracting the corresponding CMB Solar dipole variation from the AME map, and end up with a nearly identical total $\chi^2$; the orbital CMB dipole is not bright enough at 23 GHz relative to AME emission to break this degeneracy on its own.

This is illustrated in Fig. <ref>, which shows the AME amplitude map as derived for three different values of the mean -band gain, $g_0$, namely 1.175, 1.181, and $1.184\,\mathrm{du\,mK^{-1}}$; the extreme values differ only by 0.5 %. All of these three values appear equally acceptable from a pure $\chi^2$ point-of-view, relative to the noise level and modeling errors of these data. At the same time, it is clear from visual inspection that only the middle value actually makes physical sense, given what we know about the structure of the Milky Way. For this reason, we apply a Gaussian prior on the absolute -band gain of $g_0 \sim \mathcal N(1.181, 0.001^2)$ to regularize this issue. Thus, the extreme panels in Fig. <ref> represent $\pm6\,\sigma$ outliers, respectively, and will appear in our Markov chains with a frequency of about 1-in-$10^9$.

It is reasonable to ask why the  pipeline produced sensible results without applying such a prior during their calibration procedure. We posit that the answer is due to the main difference between the two approaches. While  attempts to fit a single overall parametric model to all data at once, the  pipeline calibrated each channel independently by co-adding data from one channel into a map, subtracting that map from the TOD, fitting the gain to the orbital dipole, and iterating until the solution became stable. An advantage of the single-channel approach is that the solution is independent of the assumed sky model. However, a disadvantage is that it is impossible to break any potential inherent degeneracies; it cannot be combined with external observations in any meaningful way. One important example of this regarding the  data is a strong degeneracy between the transmission imbalance factors and the polarized sky signal; it is exceedingly difficult to break this degeneracy using data from only one DA alone, and the resulting errors will propagate to most other aspects in the analysis. In the global approach, on the other hand, the polarization modes that are poorly measured by  alone are well measured by  and vice-versa, resulting in an overall better constrained fit.

Second, as reported by <cit.> for the  analysis, another important degeneracy in the current global model concerns the spectral index of polarized synchrotron emission versus the time-variable detector gain; when fitting both the polarized synchrotron amplitude and calibration freely without priors, the synchrotron spectral index at high Galactic latitudes tend to be biased toward unreasonably flat values, $\beta_\mathrm s \lesssim -2.5$, which was likely due to a low level of unmodeled systematics, for instance temperature-to-polarization leakage, rather than true polarized synchrotron emission. In turn, this resulted in a contaminated CMB sky map with a strong synchrotron morphology. To break this degeneracy, <cit.> chose to marginalize the high-latitude synchrotron spectral index over a Gaussian prior of $\mathcal N(-3.30,0.1^2)$, informed by <cit.>, rather than estimate it from the data themselves. We observe the same degeneracy, and the introduction of the -band data is not sufficient to break it on its own. For this reason, we choose to apply the same informative prior.

Third and finally, we also marginalize over the AME scale index with a prior of $\beta_{\mathrm{AME}}\sim\mathcal N(3.56, 0.1^2)$. The parameters of these priors were determined by running a grid over $\beta_{\mathrm{AME}}$, and identifying the range that resulted in reasonable residuals near the Galactic plane, similar to that shown in Fig. <ref> for the absolute calibration of -band. We note that this prior should in principle be replaced with direct $\chi^2$-based posterior optimization, combined with a properly tailored analysis mask. However, the recent release of the QUIJOTE data <cit.>, which covers the 11–19 GHz frequency range, suggests that the entire AME model should be revisited in a future joint +LFI+QUIJOTE analysis. We therefore leave detailed prior and SED optimization to that work. For further information regarding AME modeling with the current dataset, we refer the interested reader to <cit.>.

In sum, we impose strong priors on all foreground spectral indices, with parameters that are informed by the requirement of obtaining physically meaningful component maps. These strong priors imply that the foreground spectral parameters sampled within the chain do not carry independent significance in the traditional posterior sense, but are in practice only nuisance parameters used to marginalize over externally defined uncertainties.

§.§ Posterior distribution and Gibbs sampling

As shown by <cit.>, this joint parametric description of the instrumental effects and sky allows us to write down a total model for the data, $\boldsymbol d=\boldsymbol s^\mathrm{tot}(\boldsymbol\omega)+\boldsymbol n^\mathrm w$, where $\boldsymbol s^\mathrm{tot}$ encompasses all of the terms in Eq. (<ref>) except for the white noise term. Assuming that all instrumental effects have been modeled adequately, and that the white noise is Gaussian distributed, the data should then also be Gaussian distributed with a mean of $\boldsymbol s^\mathrm{tot}(\boldsymbol\omega)$ and variance $\boldsymbol \sigma_0^2$. In general, the likelihood reads

\begin{equation}
P(\boldsymbol d\mid\boldsymbol\omega)\propto\exp\left(-\frac12\sum_t\frac{(d_t-s^\mathrm{tot}_t(\boldsymbol\omega))^2}{\sigma_0^2}
\right).
\end{equation}

If $\boldsymbol d\sim\mathcal N(\boldsymbol s^\mathrm{tot},\boldsymbol\sigma_0^2)$ is the correct model for the data, the argument of the exponent is proportional to a $\chi^2$-distribution with $n_\mathrm{TOD}$ degrees of freedom, where $n_\mathrm{TOD}$ number of datapoints within a given datastream. In the limit of large $n$, a $\chi^2$ distribution is well-approximated by a Gaussian with mean $n$ and variance $2n$. Therefore we define and use the reduced normalized $\chi^2$ statistic,

\begin{equation}
\chi^2\equiv \frac{\sum_t((d_t-s_t^\mathrm{tot})/\sigma_0)^2 - n_\mathrm{TOD}}{\sqrt{2n_\mathrm{TOD}}},
\label{eq:chisq}
\end{equation}

which is approximately drawn from the standard normal distribution $\mathcal N(0,1)$.

Following <cit.>, the  Gibbs chain for this analysis is given by

s^inst     ←P(s^inst                 |,     s^inst,     ,     x_im,     ,     ξ_n,          β,     a,     C_ℓ)

←P(                   |,     s^inst,     ,     x_im,     ,     ξ_n,          β,     a,     C_ℓ)

x_im     ←P(x_im     |,     s^inst,     ,     x_im,     ,     ξ_n,          β,     a,     C_ℓ)

←P(            |,     s^inst,     ,     x_im,     ,     ξ_n,          β,     a,     C_ℓ)

ξ_n     ←P(ξ_n                |,     s^inst,     ,     x_im,     ,     ξ_n,          β,     a,     C_ℓ)

β    ←P(β)

a     ←P(a                |,     s^inst,     ,     x_im,     ,     ξ_n,          β,     a,     C_ℓ)

C_ℓ     ←P(C_ℓ             |a)

with each step requiring its own dedicated sampling algorithm. The  pipeline is designed so that results of each Gibbs sample can be easily passed to each other, and that the internal calculations of each step do not directly depend on the inner workings of each other, which greatly increases modularity of the code.

§.§ Sampling algorithms

Before we discuss the results of this Gibbs chain as applied to the  LFI and  data, we summarize the TOD processing steps in this section. Each step of the Gibbs chain requires its own conditional distribution sampling algorithm. In Sect. <ref> we review the sampling algorithms implemented in the  suite of papers, while Sects. <ref>–<ref> provide an overview of the -specific processing steps.

§.§.§ Review of sampling algorithms

Most of the techniques required for  data analysis have already been  described in the  project and implemented in . This section includes a summary of the algorithms that were used previously for the analysis of LFI data. In each of these cases, every part of the model not explicitly mentioned is held fixed unless specified otherwise.

Noise estimation and calibration are described by <cit.> and <cit.>, respectively. As noted in those works, these two steps are strongly correlated, simply because the timestream

\begin{equation}
\boldsymbol d_i=\boldsymbol g_i\boldsymbol s_{i}^\mathrm{tot}+\boldsymbol n_{i}^\mathrm{corr}+\boldsymbol n_{i}^\mathrm{w}
\end{equation}

may be almost equally well fit by two solutions defined schematically by

\begin{equation}
\boldsymbol g'=\boldsymbol g\boldsymbol s^\mathrm{tot}/(\boldsymbol s^\mathrm{tot})'
\end{equation}

\begin{equation}
{(\boldsymbol n^\mathrm{corr})'=\boldsymbol n^\mathrm{corr}+\boldsymbol g\boldsymbol s^\mathrm{tot}+\boldsymbol g'(\boldsymbol s^\mathrm{tot})'};
\end{equation}
the only thing that breaks this degeneracy is the noise PSD, which is a relatively loose constraint. A Gibbs sampler is not very effective for nearly degenerate distributions, and we therefore instead define a joint sampling step for the correlated noise and gain. In practice, this is done by first drawing the calibration from its marginal distribution with respect to $\boldsymbol n^\mathrm{corr}$, and then drawing $\boldsymbol n^\mathrm{corr}$ from its conditional distribution with respect to $\boldsymbol g$,

←P(    |,              ξ_n)

←P(    |,     ,         ξ_n).

One can see that this is a valid sample from the joint distribution from the definition of a conditional distribution, $P(\g,\ncorr\mid\boldsymbol\omega)=P(\ncorr\mid\g,\boldsymbol\omega)P(\g\mid\boldsymbol\omega)$. In practice, this simply means that when sampling for $\g$, the covariance matrix $\mathsf N=\mathsf N_\textrm{w}+\mathsf N_\textrm{corr}$ must be used, rather than just $\mathsf N_\mathrm{w}$.

models the gain at each timestream $t$ for a detector $i$ as

\begin{equation}
g_{t,i}=g_0+\Delta g_i+\delta g_{q,i},
\end{equation}

where $q$ labels the time interval for which we assume the gain is constant, typically a single scan. In order to sample the gain, we write down a generative model for the TOD,

\begin{equation}
\dv_i=\g_{i}\s_{i}^\mathrm{tot} +\n_{i}^\mathrm{tot}\sim\mathcal N(\g_i\s_{i}^\mathrm{tot},\mathsf N_i).
\end{equation}

Since $\dv_i$ is given as a linear combination of the fixed signal and the gains, a random sample of the gain can be drawn by solving[See, e.g., Appendix A.2 of <cit.> for a derivation of this result.]

\begin{equation}
\label{eq:gain_gauss_samp}
[(\s_i^\mathrm{tot})^T\mathsf N_i^{-1}\s_i^\mathrm{tot}]\g_i
\end{equation}

where $\boldsymbol\eta\sim\mathcal N(0,1)$ is a vector of standard normal variables.
Note that $\mathsf N_i$ depends implicitly on the noise PSD $\boldsymbol \xi_n$, while the fluctuations corresponding to $\ncorr$ are properly downweighted by $\mathsf N_{\mathrm{corr},i}$.
As detailed by <cit.>,  samples $g_0$, $\Delta g_i$, and $\delta g_{q,i}$ in separate steps. Specifically, the absolute calibration $g_0$, for the CMB-dominated channels, is fitted using only the orbital dipole, while the relative calibrations, $\Delta g_i$, exploits the full sky signal. The same is true for the time-dependent gain fluctuations, $\delta g_{q,i}$, and in this case an additional smoothness prior is applied through an effective Wiener filter. The Gibbs chain is formally broken by fitting the absolute gain $g_0$ to the orbital dipole alone, as opposed to the full sky signal. However, this makes the sampling more robust with respect to unmodeled systematic effects, somewhat analogous to applying a confidence mask when estimating the CMB power spectrum.

The correlated noise sampling, described by <cit.>, follows a similar procedure, except this now conditions upon the previous gain estimate, which is sampled immediately before the correlated noise component in the code. Similar to the gain case, we can write a generative model for the data,

\begin{equation}
\dv_i=\g_i\s_i^\mathrm{tot}+\n_i^\mathrm{corr}+\n^\mathrm{w}_i
\sim\mathcal N(\g_i\s_i^\mathrm{tot},\N_{\mathrm{corr},i}+\N_{\mathrm{w},i}).
\end{equation}

Given fixed $\boldsymbol r_i=\dv_i-\g_i\s_i^\mathrm{tot}$, we can again write a sampling equation,

\begin{equation}
\label{eq:ncorr_gauss_samp}
\n^\mathrm{corr}_i
=\N^{-1}_{\mathrm{w},i}\boldsymbol r_i
\end{equation}

This gives a sample of the underlying correlated noise.

To sample the correlated noise parameters, we assume that the correlated noise is drawn from a correlated Gaussian and from the conditional posterior distribution,

\begin{equation}
\end{equation}

where $P(\boldsymbol\xi _n)$ is a flat prior over the PSD parameters.
The simplest and most commonly used parametrization for correlated noise is given by

\begin{equation}
\mathsf N_\mathrm{corr}(f)=\sigma_0^2\left(\frac f{f_\mathrm{knee}}\right)^\alpha.
\end{equation}

This can in principal be modified, and for  LFI a Gaussian log-normal bump was added at a late stage in the  analysis. Rather than sampling for $\sigma_0$, we effectively fix the white noise level to the noise level at the highest frequency, e.g.,

\begin{equation}
\sigma_0^2\equiv\frac{\mathrm{Var}(r_{t+1}-r_t)}2,
\end{equation}

where $t$ and $t+1$ are consecutive time samples, and ${\boldsymbol r\equiv\boldsymbol d-\boldsymbol g\boldsymbol s^\mathrm{tot}
-\boldsymbol n^\mathrm{corr}}$. In practice, this makes $\sigma_0$ a deterministic function of the sampled sky and gain parameters. The parameters $\alpha$ and $f_\mathrm{knee}$ are not linear in the data, but they can be sampled efficiently using a standard inversion sampler (see, e.g., Appendix A.3 of  or Chapter 7.3.2 of  for further details). In practice, this requires computing the posterior over a linear grid one parameter at a time.

Once the instrumental parameters have been sampled,  computes the calibrated TOD for each band,

\begin{equation}
+s_{t,i}^\mathrm{fsl}+\delta s_{t,i}^\mathrm{leak}+{s}_{t,i}^\mathrm{inst}\right)
\label{eq:caltod}
\end{equation}

where $\s^\mathrm{orb}$ is the orbital dipole <cit.>, $\s^\mathrm{fsl}$ is the far sidelobe timestream <cit.>, $\delta\s^\mathrm{leak}$ is the bandpass leakage <cit.>, and $\s^\mathrm{inst}$ is some instrumental-specific contribution, e.g., the 1 Hz electronic spike for LFI. With a correlated noise realization removed, one can perform simple binned mapmaking, weighting each pixel by the white noise amplitude.

§.§.§ Differential mapmaking

The first additional algorithm that needs to be added to  in order to process  TOD data is support for differential mapmaking <cit.>. After calibration and correction for instrumental effects, the TOD can be modeled as

\begin{equation}
\boldsymbol d=\mathsf P\boldsymbol m+\boldsymbol n^\mathrm{w},
\end{equation}

\begin{equation}
\boldsymbol m=\mathsf B^\mathrm{symm}\mathsf M\boldsymbol a
\end{equation}

is the expected map for each detector after removing the orbital dipole, far sidelobe, baseline, and a realization of correlated noise. The differential pointing strategy can be represented in matrix form as

\begin{align}
\label{eq:wmap_pointing}
\mathsf P_{tp}&=
\\
\end{align}

where $p_\A$ and $p_\B$ are the time-dependent pointings for each DA. Note that this is equivalent to Eqs. (<ref>) and (<ref>), except without the spurious component $S_{\A/\B}$ as the bandpass mismatch is explicitly subtracted beforehand in Eq. (<ref>). The maximum likelihood map can now in principle be derived using the usual mapmaking equation,

\begin{equation}
\label{eq:mapmapking_eqn1}
\mathsf P^T\mathsf N^{-1}\mathsf P\boldsymbol m=\mathsf P^T\mathsf N^{-1}\boldsymbol d.
\end{equation}

For a single-horn experiment, i.e.,  LFI, this reduces to a $3\times3$ matrix that can be inverted for each pixel independently. For the pointing matrix in Eq. (<ref>), this is no longer possible, as there is inherently coupling between horns A and B in the timestreams. The $3N_\mathrm{pix}\times3N_\mathrm{pix}$ matrix can in principle be solved using an iterative algorithm, e.g., preconditioned conjugate gradients <cit.>.

<cit.> identified an issue where a large difference in the sky temperature values at pixel A versus pixel B induced artifacts in the mapmaking procedure. We adopt the procedure first described by <cit.> where only the pixel in a bright region, defined by a small processing mask <cit.> is accumulated, thus modifying the mapmaking equation to

\begin{equation}
\mathsf P^T_\mathrm{am}\mathsf N^{-1}\mathsf P\boldsymbol m
=\mathsf P^T_\mathrm{am}\mathsf N^{-1}\boldsymbol d.
\label{eq:bicg_stab}
\end{equation}

This equation can be solved using the BiCG-STAB algorithm for a non-symmetric matrix $\mathsf A$ where $\mathsf A\boldsymbol x=\boldsymbol b$. We apply a preconditioner $\mathsf M$ by numerically inverting the same problem with $N_\mathrm{side}=16$ maps and applying a diagonal noise matrix. Numerically, we define convergence as when the residual $\boldsymbol r\equiv\boldsymbol b-\mathsf A\boldsymbol x$ satisfies $\boldsymbol r^T\mathsf M^{-1}\boldsymbol r/\boldsymbol b^T\mathsf M^{-1}\boldsymbol b<10^{-10}$, which typically takes about 20 iterations for producing frequency maps.

The full noise covariance matrix $\mathsf N_{pp'}$ is given by the inverse of $\mathsf P^T_\mathrm{am}\mathsf N^{-1}\mathsf P$, where the diagonals $\mathsf N_{pp}$ are the white noise variance for each Stokes parameter. An additional quantity that was computed in  but not delivered in the final products is the covariance of the Stokes parameters within a single pixel, $\sigma_{QU,p}$. We find that the correlation between Stokes parameters, $\rho_{QU}\equiv\sigma_{QU}/\sqrt{\sigma_{QQ}\sigma_{UU}}$, is of order 0.5 for the  DAs, as shown in Fig. <ref>. For  LFI, the 30 and 70 GHz channels have $|\rho_{QU}|\sim0.1$, while the 44 GHz correlations are notably higher with $|\rho_{QU}|\sim0.5$. The reason for this difference is that the 44 GHz channel has three horns. Two of those are aligned with the scanning direction in the focal plane, and have polarization angles that are rotated by $45^{\circ}$ with respect to each other. Together those two horns therefore disentangle polarization information very efficiently. The third horn, however, does not have a corresponding partner, and relies only on satellite precession to recover individual Stokes parameters. For comparison, all 30 and 70 GHz horns have partners aligned with the scanning direction. In the current work, we have implemented support for the full $3\times 3$ noise matrices, including $\sigma_{QU},$ for component separation and map-based $\chi^2$ calculations for both  and LFI.

Cross-correlation $\rho_{QU}$ per pixel for (top:)  and (bottom:) LFI 30 GHz.

§.§.§ Baseline sampling

The data model adopted by <cit.> can be written in raw du as

\begin{equation}
\data = \mathsf{GPBM}\boldsymbol a+\boldsymbol n+\boldsymbol b,
\end{equation}

where $\boldsymbol b$ is the instrumental baseline and $\boldsymbol n$ is the
total instrumental noise. As noted above,  divides
the noise into $\n=\n^\mathrm w+\n^\mathrm{corr}$, a white noise term and a
correlated noise term. Because the white noise is by definition uncorrelated in time, it does not have any
correlations between adjacent pixels, so that any pixel-pixel covariance should
be fully described by realizations of the $\n^\mathrm{corr}$ timestream.

estimates the baseline using the full estimate of the current sky
model, $\boldsymbol r=\boldsymbol d-\boldsymbol g\boldsymbol s^\mathrm{tot}=\boldsymbol
b+\boldsymbol n$. Modeling $\boldsymbol b=b_0+b_1\Delta t$, we solve for $b_0$
and $b_1$ using linear regression in each timestream while masking out samples
that lie within the processing mask. Strictly speaking, this is breaking the
Gibbs chain, as we are not formally sampling $b_0$ and $b_1$ for each TOD
chunk. In practice, baseline estimation uncertainty propagates to correlated
noise realizations and PSD parameters, as discussed below.

The approach detailed by <cit.> and the
implementation differ mainly in two ways. First, the
assumed stable timescales are different – the initial
baseline is estimated over one hour timescales, and assumed to
be an actual constant, whereas  assumes constant values
through the entire time chunk, which is 3–7 days depending on the
band in question, but allows a linear term in the baseline. Second,
the two methods differ in how they treat nonlinear residuals in the
first-order baseline model.  As noted by <cit.>,
residual baseline variations manifest as correlated noise stripes in
the final maps, and  solves this using a time domain filter,
downweighting the data based on the noise characterization. This
is similar to the  approach, which
accounts for this as part of the correlated noise component. The main
advantages of the latter is that it allows for proper error
propagation at all angular scales without the use of a dense
pixel-pixel noise covariance, and provides a convenient
means for inspecting the residuals visually by binning the correlated
noise into a sky map.

§.§.§ Transmission imbalance estimation

Transmission imbalance, the differential power transmission of the optics and waveguide components between horns A and B, can be parameterized as

\begin{equation}
\end{equation}

This can be decomposed into a differential (d) and common-mode (c) signal such that

\begin{equation}
d_{t,i}=g_{t,i}[s_{t,i}^\mathrm d+x_{\mathrm{im},i}s_{t,i}^\mathrm c]+n_t.
\end{equation}

In this form, the imbalance parameters can be estimated by drawing Gaussian samples from the standard mean and standard deviation over the entire mission. To draw samples for $x_{\mathrm{im},i}$, we construct a sampling routine analogous to the gain estimation of Eq. (<ref>) and correlated noise estimation of  Eq. (<ref>), with $\boldsymbol r=\boldsymbol d-\boldsymbol g\boldsymbol s^\mathrm d$,

\begin{equation}
[(\boldsymbol g\boldsymbol s^\mathrm c)^T\mathsf N^{-1}\boldsymbol g\boldsymbol s^\mathrm c]x_\mathrm{im}
=(\boldsymbol g\boldsymbol s^\mathrm c)^T\mathsf N^{-1}\boldsymbol r+(\boldsymbol g\boldsymbol s^\mathrm c)^T\mathsf N^{-1/2}\boldsymbol\eta,
\end{equation}

cross-correlating the common-mode signal with $\boldsymbol r$ with appropriate weights and adding a Gaussian random variable with the correct weighting. Note that we are marginalizing over the correlated noise here by using $\N=\N_\mathrm{w}+\N_\mathrm{corr}$. This mitigates any baseline drifts being erroneously attributed to the common-mode signal and biasing the estimate of $x_\mathrm{im}$.

The  procedure, described by <cit.>, fit for common-mode and differential coefficients along with a cubic baseline over 10 precession periods at a time, corresponding to 10 hours of observation. The mean and uncertainty were then calculated by averaging and taking the standard deviation of these values. This approach has the benefit of allowing for the tracking of possible transmission imbalance variation throughout the mission. However, none of the  suite of papers have found evidence for this, and it has not arisen in our analysis, so we model this as an effect whose value is constant throughout the mission.

§ DATA AND DATA PROCESSING

flagging statistics for each DA. The second column indicates the fraction of data that are removed by the official  flags, while the third column indicates the fraction that is additionally discarded in the current processing for computational reasons. The fourth column indicates the total fraction of data actually used to generate the final maps.

Band     Flagged (%)     Discarded (%)     Used (%)

1.72     0.87     97.4

1.64     0.88     97.5

1      1.84     0.84     96.5

2      1.62     0.81     97.6

1      1.62     1.10     97.3

2      1.61     1.01     97.4

1      1.76     1.03     97.2

2      1.60     0.81     97.6

3      1.61     0.87     97.5

4      1.60     0.81     97.6

Computational resources required for end-to-end
processing. All times correspond to CPU hours, and all data volumes are reported in GB. Reported times are
averaged over more than 100 samples, and vary by $\lesssim\,5\,\%$ from sample to
sample. Note that the average cost per sample takes into account the undersampling of 70, , and .

to 5.0cm#1em

5Data volume2pt
10pt Compressed TOD volume     **86    *178    **597    **13    **12    **15    **15    **19    **18    **26   **26    **26    **26    1 0532pt
5Processing time (cost per run)10pt TOD initialization/IO time                       1.8    2.5    7.8    0.7    0.6    0.8    0.7    0.9    0.8    1.3    1.3    1.0    0.9    *21.110pt Other initialization                                                                                *14.610pt Total initialization                                                                                *35.72pt
5Gibbs sampling steps (cost per sample)10pt Huffman decompression                                1.2    2.2    23.2    0.8    0.9    1.1    1.1    1.5    1.4    2.0    2.0    2.0    2.0    *41.410pt Array allocation                                     0.4    0.9    51.6    1.3    1.3    1.5    1.5    3.1    3.3    4.0    3.8    4.0    4.0    *80.710pt TOD projection ($\P$ operation)                      0.9    2.0     12.3    6.1    7.1    8.7    8.9    11.4    11.3    15.9    15.8    15.7    15.8    131.910pt Sidelobe evaluation                                  1.2    2.6     9.5    3.0    3.5    4.1    4.2    5.5    5.4    7.8    7.7    7.7    7.5    *69.710pt Orbital dipole                                       0.9    2.0    9.0    1.2    1.5    1.8    1.9    2.6    2.5    3.8    3.8    3.8    3.8    *38.610pt Gain sampling                                        0.6    0.9    2.2    1.3    1.3    0.8    0.8    1.3    1.3    1.2    1.2    1.2    1.2    *15.310pt 1 Hz spike sampling                                 0.3    0.4    1.9                                            **2.710pt Correlated noise sampling                            2.1    4.3    24.8    2.7    2.9    3.7    3.8    6.2    5.4    7.7    7.4    6.9    8.3    *86.410pt Correlated noise PSD sampling                        5.0    6.2    1.6    0.3    0.3    0.3    0.3    0.5    0.5    0.7    0.6    0.6    0.7    *17.610pt TOD binning ($\P^t$ operation)                       0.1    0.1    10.5    0.8    0.8    1.0    1.0    1.7    1.6    2.4    2.4    2.4    2.4    *27.210pt Mapmaking                                                        9.2    9.7    13.1    12.7    21.7    20.2    35.4    34.9    36.1    39.3    232.310pt MPI load-balancing                          1.2    1.7    9.2    2.2    2.0    2.2    2.1    3.6    3.3    4.8    4.6    4.5    4.6   *46.010pt Sum of other TOD processing                          0.7    1.6    13.1    0.1    0.2    0.5    0.4    0.7    0.8    0.9    1.0    0.9    1.2   *22.110pt TOD processing cost per sample                 14.6    24.9    169.7     28.8    31.5    38.7    38.7    59.8    57.0    86.6    85.2    85.8    90.8    812.12pt
10pt Amplitude sampling                                                    *16.210pt Spectral index sampling                                                    *32.12pt
10pt Average cost per sample                                                                      418.94pt5pt

We describe the public  data products in Sect. <ref>, then describe the treatment we apply to make them compatible with  in Sect. <ref>. Finally, we describe the computational requirements in Sect. <ref>.

§.§ Publicly available  products

The full  dataset is hosted at the Legacy Archive for Microwave Background Data Analysis (LAMBDA).[<https://lambda.gsfc.nasa.gov/product/wmap/dr5/m_products.html>] In addition to the primary scientific products, e.g., cosmological parameters, CMB power spectra and anisotropy maps and frequency maps, the time-ordered data (TOD) can be downloaded, both in uncalibrated and calibrated form.[<https://lambda.gsfc.nasa.gov/product/wmap/dr5/tod_info.html>] In principle, thanks to these data and the explanatory supplements <cit.>, the entire data analysis pipeline can be reproduced from uncalibrated TOD to frequency maps.

For this analysis, we keep certain instrumental parameters fixed to the reported values. For example, we have made no attempts to rederive the pointing solutions, re-estimate the main beam response and far sidelobe pickup, or recover data that were flagged in the  event log. These and other analyses, such as estimating the bandpass shift over the course of the mission, are certainly possible within the larger Gibbs sampling framework. However, in this work we limit ourselves to recalibrating the TOD, estimating their noise properties, and applying bandpass corrections to the data before mapmaking.

§.§ TOD preprocessing and data selection

The full nine-year  archive spans from 10 August 2001 to 10 August 2010, with the raw uncalibrated data comprising 626 GB. A little over 1 % of the data were lost or rejected due to incomplete satellite telemetry, thermal disturbances, spacecraft anomalies, and station-keeping maneuvers, with an extra 0.1 % rejected due to planet flagging <cit.>.
The final results reported by <cit.> included roughly 98.4 % of the total data volume.
A full accounting of all data cuts can be found in Table 1.8 of <cit.>. In this analysis we flag the same data indicated in the fiducial  analysis, and use the same planet exclusion radii.

As shown by <cit.>, a large fraction of 's computational time is spent performing Fast Fourier Transforms (FFTs) on individual scans. Rather than truncating datastreams to have lengths equal to “magic numbers” for which  <cit.> performs efficiently, as was done in the  analysis,
we redistribute the data into scans of length $2^N$, where $N=22$ for –, $N=23$ for –. This yields scans with lengths of 6.21 days for - and -band, 4.97 days for -band, 7.46 days for -band, and 4.97 days for -band.[Note that scans with equal $n_\mathrm{TOD}$ cover different lengths of time due to the different sampling rate for each frequency.]
These datastream lengths are short enough to be processed quickly and distributed efficiently across multiple processors, while being long enough to properly characterize the noise properties of the timestreams, whose $f_\mathrm{knee}$ values are on the order $1\,\mathrm{mHz}$. Most importantly,  performs fastest when the datastream is of length $2^N$.

When redistributing the data, timestreams of length $2^N$ were interrupted by events logged in Table 1.8 of <cit.>.
When we encountered these events, interrupted TOD segments were appended to the previous TOD, in most cases creating TODs with lengths $>2^N$. We found that events of length $<2^N$ were too short to accurately estimate the noise PSD parameters. This criterion led us to discard these otherwise useful data. In addition, when $>10\,\%$ of the TOD are flagged, the large number of gaps in the data makes the solution of Eq. (<ref>) for $\ncorr$ computationally more expensive. Given that data near many large gaps are more likely to have unmodeled effects than stable data, and they are more expensive to process, we chose to remove these from the analysis. Together, these two effects led to $\simeq1\%$ of the data to be discarded. We summarize the full flagging statistics for our maps in Table <ref>. In total, the  maps use about 1 % less data than the WMAP9 official products. The total difference in data volume can be entirely accounted for by the cuts described in this paragraph.

§.§ Computational resources and future plans

A key motivation of the current analysis is to evaluate whether it is feasible to perform a joint analysis of two datasets simultaneously, each with its own particular processing requirements and algorithmic treatment. One of the results from <cit.> was that most of the data processing procedures for   and  LFI overlapped, with the notable exception of mapmaking. While the algorithmic requirements have been discussed in Sect. <ref>, we have not yet quantified the requirements in terms of RAM and CPU hours. In Table <ref>, we enumerate the RAM requirements and CPU time for each sampling step using a single AMD EPYC 7H12, 2.6 GHz cluster node with 128 cores and 2 TB of memory. As such, approximate wall runtimes can be obtained by dividing all numbers in Table <ref> by 128.

Despite the relatively small data volume spanned by , e.g., 86 GB for 30 GHz vs. 13 GB for -band, the CPU time is comparable to each of the LFI channels.  The single largest reason for this is the mapmaking step, Eq. (<ref>), which requires looping over the entire dataset for each matrix multiplication, a process which must be repeated $\sim20$ times. As discussed in Sec. <ref>, this is vastly sped up by the use of a low resolution preconditioner, reducing the number of iterations by an order of magnitude.

Additionally, operations that require the creation of  timestreams for each detector, i.e., TOD projection, sidelobe evaluation, and orbital dipole projection, take much longer than expected from a pure data volume scaling. Part of this is due to  evaluating the sky in two pixels simultaneously, doubling the expected workload, but the other issue is that we are unable to benefit from the ring-clustering based TOD distribution scheme used for LFI. Due to 's more complex scan strategy and detector geometry, it is impossible to cluster scans with similar pixel coverage onto a single core, which makes pixel-space lookup operations less efficient in this case.

Gain sampling and correlated noise sampling include multiple FFTs. Typical LFI TODs are of length $\sim200\,000$, an order of magnitude smaller than the  TODs of length $\sim5\,000\,000$. Despite the TOD lengths being pre-determined to be $2^N$, this extra length still results in longer run times for equivalent data volumes, but does yield noise information on much longer time scales than we have for LFI. Note that  was had typical $f_\mathrm{knee}$'s over 100 times smaller than LFI's, so TODs that were over 100 times longer are necessary for characterizing its noise PSD properties.

Trace plots of the 113 gain and noise parameters for a single scan starting on MJD 52285.2. The two colors correspond to the two independent Markov chains produced in this analysis.

Noise parameter correlation matrix. We average over all Gibbs samples of the noise parameters $\boldsymbol\xi_n=\{\alpha,f_\mathrm{knee},\sigma_0\}$ for each PID. We then find the correlation in time between these averages for the different bands and detector. The results here are for the calibrated white noise level, $\sigma_0[\mathrm{mK}]$. The values for each detector are ordered 13, 14, 23, and 24.

For the current analysis, which aims primarily to derive posterior-based  frequency maps, we produce a total of 500 main Gibbs samples, divided into two chains. Noting that the computational cost of the -channel carries almost half of the total expense of the  TOD processing, while being of less scientific importance than, say, the -band, we choose to only reprocess this channel every fourth main sample. Likewise, we only reprocess the -band every other main sample, and the LFI 70 GHz sample every fourth sample. The total cost for producing 500  , , ,  30, and 44 GHz samples, 250 -band samples, and 125 -band and 70 GHz samples is 210k CPU-hrs, and the total wall time is 33 days. Noting that the  analysis required 4000 samples to reach full convergence in terms of the optical depth of reionization <cit.>, a corresponding complete LFI+ analysis will cost about 1.7M CPU-hrs, and take about 9 months of continuous runtime on two cluster nodes. While entirely feasible, this is sufficiently expensive that we choose to perform the analysis in two stages; first we present preliminary frequency maps in the current paper, and use these to identify potential outstanding issues, either in terms of data model or Markov chain stability. An important goal of this phase is also to invite the larger community to study these preliminary maps, and thereby identify additional problems that we may have missed. Then, when all issues appear to have been resolved, we will restart the process, and generate sufficient samples to achieve full convergence.

§ INSTRUMENTAL PARAMETERS

In this section and Sect. <ref> we present the main results from the  DR1 analysis, which may be summarized in terms of the joint posterior distribution. For organizational purposes, we will discuss instrumental parameters, frequency maps, and astrophysical results separately in this and the following two sections. It is important to remember that these results are all derived from one single highly multivariate posterior distribution, and every parameter is in principle correlated with all others. In this section, we focus on instrumental parameters, starting with visual inspection of the basic Markov chains and posterior means, before considering each instrumental parameter in turn.

§.§ Markov chains, correlations and posterior mean statistics

To build intuition regarding the general Markov chain properties, we show in Fig. <ref> the Markov chains for the gain and noise parameters for one arbitrary diode (113) and scan. Each panel corresponds to one single parameter, and the observed variation quantify the uncertainty in that single parameter due to the combination of white noise and correlations with other parameters. Here we see that the different parameters have quite different correlation lengths; the gain (in the top panel) has a very short autocorrelation length, as in just a few samples, while the noise parameters have typical correlation lengths of a few tens of samples. Even for these parameters, however, the full set of 500 samples provides a fairly robust estimate of the full marginal mean and uncertainty.

The bottom panel shows the reduced normalized $\chi^2$ for the same scan in units of $\sigma\equiv\sqrt{2n_\mathrm{TOD}}$, and we see that this also shows similar correlation lengths as the noise parameters. This makes sense since the TOD residual at the level of a single sample is strongly noise dominated. In contrast, small variations in either the sky signal or gain have relatively small impacts on this particular $\chi^2$; the goodness of fit of such global parameters is better measured through map-level residuals and $\chi^2$'s. In this respect, we also note that the absolute value of the TOD-level $\chi^2$ is for this particular scan about $-7.5\,\sigma$, which at first sight appears as a major goodness of fit failure. However, it is important to recall that a typical scan contains about five million data points, and this statistic is therefore extremely sensitive to any deviation in the noise model. Specifically, the reduced $\chi^2$ for this particular scan is $\chi_{\mathrm{raw}}^2/n_{\mathrm{TOD}}=0.993$, which corresponds to an overestimation of the white noise level of only 0.3 %. Furthermore, as discussed in Sect. <ref>, we currently assume a strict $1/f$ noise model for the  noise, while the true  noise is known to exhibit a very slight non-white noise excess at high frequencies <cit.>. Properly modeling such non-white high-frequency noise is therefore an important goal for the next  data release.  Such work is also a vital step in preparing for integration of other types of experiments with non-white noise into the framework, such as  HFI. However, in absolute terms, the impact of this model failure is very limited, and not likely to significantly affect any astrophysical results; it is primarily a limitation for TOD-level goodness of fit testing.

For a survey of the entire experiment's noise properties, Fig. <ref> shows pairwise correlations between the various noise parameters for all DAs, averaged over all Gibbs samples and scans. It is important to note that a non-zero correlation in this plot does not indicate that the specific noise realization is correlated between DAs, but only that the noise PSD parameters are correlated. This is expected due to the  satellite motion around the Sun, which induces an annual variation in the system temperature. This correlation plot therefore primarily quantifies the sensitivity to this common-mode signal for each diode. Most notably, we see that the 2 DA exhibits particularly strong correlations, and we also note that the calibrated white noise $\sigma_0[\mathrm{mK}]=\sigma_0[\mathrm{du}]/g$ is generally more susceptible to these variations than $f_{\mathrm{knee}}$ and $\alpha$.

Overview of 113. The red solid lines in the first and second panel  are the regressed gain and baseline from , while the black lines in all panels are posterior means from the  Gibbs chain. The red dashed and yellow dashed lines are reported $\sigma_0$ and $f_\mathrm{knee}$ values from the first-year  data analysis and GSFC measurements, respectively.

Next, in Fig. <ref> we show posterior mean values for each instrumental parameter for the same 113 diode, in this case plotted as a function of time throughout the entire mission. The panels show, from top to bottom, 1) gain; 2) the difference between the baseline mean and its full-mission average; 3) the baseline slope; 4) the white noise level; 5) the correlated noise knee frequency; 6) the correlated noise slope; and 7) the TOD-level $\chi^2$. The  results are shown as black curves, while the  results are (for the gain and baseline) shown as red curves; dotted red and orange line corresponds to the first-year  and Goddard Space Flight Center (GSFC) laboratory measurements, respectively. Note that the gain and baseline are nearly indistinguishable – their differences will be discussed in more detail in Sect. <ref>. For brevity, we have only shown the results for one single diode here. However, a complete survey of all instrumental parameter posterior means for all 40 diodes is provided in Appendix <ref>, and all individual samples are also available in a digital format as part of the  DR1.

§.§ Gain and baselines

We now consider the gain and baseline parameters in greater detail, and aim to compare our estimates with the  products. The  gain and baseline estimates are not directly available in terms of public data products, but only in terms of the general parametric models. For instance, the  gain model reads <cit.>

\begin{equation}
\label{eq:wmap_gain}
g=\alpha\frac{\overline V-V_\circ-\beta(T_\mathrm{RXB}-290\,\mathrm K)}
{T_\mathrm{FPA}-T_\circ}+(m\Delta t+c),
\end{equation}

where $\overline V$ represents the radio frequency bias powers per detector; $T_\mathrm{RXB}$ and $T_\mathrm{FPA}$ are the receiver box and focal plane assembly temperatures, which are recorded every 23.04 s; $\alpha$, $\V_\circ$, $\beta$, $T_\circ$, $m$, and $c$ are all free parameters that are fit to a constant value across the mission for each diode. Evaluating this model as a function of $T_\mathrm{RXB}$ and $T_\mathrm{FPA}$ requires the housekeeping data for the thermistor that was physically closest to the relevant radiometer's focal plane on the satellite. The free parameters are fully tabulated in the  Explanatory Supplement <cit.>, but the physical layout of the thermistors in the focal plane and receiver box is not readily available. We therefore do not attempt to reproduce the gain model given in Eq. (<ref>).

Difference between the  $\boldsymbol d_\mathrm{cal}=\boldsymbol d/g-\boldsymbol b - \boldsymbol s_\mathrm{sl}$ and the delivered calibrated TOD from . Note that this is for a small segment of the 113 TOD displayed in Fig. <ref>.

Rather, we estimate the gains and baselines by comparing the uncalibrated  data with the calibrated  data, after subtracting a far sidelobe contribution convolved with the delivered  DA maps plus the Solar dipole.  We find that the calibrated and uncalibrated data can be related by

\begin{equation}
\end{equation}

where the second term is a cubic polynomial with coefficients $c_i$ referenced to the time at the beginning of the scan $t_0$. The red curves in the top two panels of Fig. <ref> correspond to these estimates. The agreement between the  and  gain and baseline models appear reasonable at this level and for this diode.

A complete comparison between the   and  gain and baseline models for all diodes is provided in Appendix <ref>. In particular, Fig. <ref> shows the baseline differences as a function of time, and here we see that most diode differences scatter around a constant value that is close to zero; the precise constant value is of limited importance, since that only corresponds to a difference in the overall monopole of the maps, which for  is determined through post-processing. However, there are a few notable features. First, we see that the two 11 diodes exhibit large variations at the very beginning of the mission, with typical values of a few du's. Individual scans show notable spikes for many diodes. These are all relatively isolated in time, and will therefore have relatively minor impact on the final maps. Far more significant are the -band differences, for which one sees both slow drifts and abrupt changes. Furthermore, in many cases they vary notably between diodes within the same DA, which translates into differences in the large-scale polarization maps derived from the two pipelines.

Similarly, Fig. <ref> compares the gain solutions directly, while Fig. <ref> shows the fractional differences in units of percent. Overall, we see that the two gain models agree to typically about 0.5 % in an absolute sense, and better than typically 0.1 % in terms of relative agreement between neighboring scans. The most striking feature in this plot is an annual variation that traces the  satellite's motion around the Sun. In general, such an oscillatory gain behaviour is entirely expected, because of known temperature variations in the satellite. However, the difficulty lies in estimating the magnitude of the oscillations, as different radiometers can respond differently to these temperature variations. In this respect, it is useful to recall that the  and  gain estimation algorithms differ at a fundamental level. The  analysis considers each DA in isolation, fitting seven instrumental parameters, defined by Eq. (<ref>), to the orbital dipole seen by each DA. The  analysis considers the problem globally, and attempts to fit all gain parameters to the full sky signal (including both the Solar and orbital CMB dipole) simultaneously, without the use of a strong instrumental model prior. Returning to the absolute gains shown in Fig. <ref>, it is difficult to determine visually which approach is better at this level alone, as the two models are quite similar; in some cases, such as 124 and 214, the  model oscillates more strongly than the  model, while in others, such as 113 and 114, the opposite is true. We also see the impact of the strong instrumental priors in the  solution particularly well in -band, where the  gains are far more noisy than the  gains.

The impact of these differences at the TOD level is illustrated in Fig. <ref>, which shows the calibrated  timestream $\boldsymbol d/g-\boldsymbol s_\mathrm{sl}-\boldsymbol b$ minus the  calibrated signal in units of microkelvin. The most prominent feature is a $\sim25\,\mathrm{\mu K}$ offset, which is unsurprising, given the different treatment of baselines in our two pipelines. The second obvious difference is a series of spikes associated with Galactic plane crossings. Differences of order $50\,\mathrm{\mu K}$ are seen where the absolute sky brightness is about $10\,\mathrm{mK}$, equivalent to $\sim0.5\,\%$ deviations in the gain solution. This is twice as large as the 0.2 % uncertainty estimated in <cit.> based on end-to-end simulations.

Another interesting feature in Fig. <ref> is slow correlated variations at a timescale of $\sim$ 20 sec timescale. There is nothing in the  instrument model that varies on such short timescales, and this must therefore come from . The most likely explanation is the fact that the  gain model depends directly on housekeeping data that are recorded with a 23.04 sec sample rate, and these values appear to have been applied without any smoothing, resulting in sharp jumps in the final  gain model, as well as increases in the observed level of scatter. At the same time, it is also important to note that the  gain model does not include any time-varying structure within a single scan, and if any artifacts resulting from this are identified in the current products, it may be worth incorporating housekeeping data in a future   data release.

§.§ Transmission imbalance

Transmission imbalance parameters for each  radiometer as estimated in the current analysis (second column) and in the official 9-year  analysis (third column). Our uncertainties indicate $1\,\sigma$ marginal posterior standard deviations.

to 1.8cm#2em
Radiometer     x_im^CG    x_im^𝑊𝑀𝐴𝑃3pt5pt
11       0.00028 ±0.00019       -0.00067 ±0.00017- 12       0.00442 ±0.00054        0.00536 ±0.00014 11       0.00338 ±0.00027        0.00353 ±0.00017 12       0.00138 ±0.00022        0.00154 ±0.00008 11       0.00076 ±0.00043       -0.00013 ±0.00046- 12       0.00530 ±0.00041        0.00414 ±0.00025 21       0.00993 ±0.00057        0.00756 ±0.00052 22       0.01314 ±0.00071        0.00986 ±0.00115 11      -0.00016 ±0.00066-        0.00053 ±0.00020 12       0.00221 ±0.00063        0.00250 ±0.00057 21       0.00281 ±0.00059        0.00352 ±0.00033 22       0.00411 ±0.00087        0.00245 ±0.00098 11       0.01145 ±0.00146        0.01134 ±0.00199 12       0.00338 ±0.00131        0.00173 ±0.00036 21       0.01534 ±0.00180        0.01017 ±0.00216 22       0.01618 ±0.00155        0.01142 ±0.00121 31      -0.00105 ±0.00107-       -0.00122 ±0.00062- 32       0.00375 ±0.00133        0.00463 ±0.00041 41       0.02596 ±0.00252        0.02311 ±0.00380 42       0.01886 ±0.00203        0.02054 ±0.00202 5pt5pt

Closely related to the gain model is the transmission imbalance factor,
$x_{\mathrm{im}}$, quantifying the difference between responsivity in the two
horns, as described in Sect. <ref>. These are listed for
each radiometer in Table <ref> for both  and ; for
the reported values correspond to marginal posterior means and
standard deviations. The same information is plotted in Fig. <ref>.
Overall, there is a reasonable agreement between the  and
estimates, with the same overall magnitude for each individual
radimoeter.  At the same time, we do observe several notable differences.
There are nearly $4\,\sigma$ differences between the derived 11 radiometer,
while for -band and 2, the $x_\mathrm{im}$ factors are consistently larger
by about $1\,\sigma$ for all radiometers. The comparison between -band and
-band is especially noteworthy, as these are the bands with the most
prominent transmission imbalance modes in the  maps, especially 4.
The uncertainties in  are also larger for the lower frequency
channels than in , which can be explained by the dependence of
$x_\mathrm{im}$ on the varying sky model within . At the same
time, the ratio of uncertainties in -band varies across radiometers, somewhat
depending on the knee frequency of the radiometer in question. Overall, the
variation in the noise levels and values is expected because of the larger
number of degrees of freedom in the  model, while the amplitude of
the relative agreement shows the robustness of the data to the specific
pipeline choices.

Comparison of the transmission imbalance factors, $x_{\mathrm{im}}$, estimated by  (dark colors) and  (light colors) for each radiometer.

While these differences may appear to be subtle, they could still be important both in terms of final sky map biases and error propagation. The reason for this is that $x_{\mathrm{im}}$ couples directly to the astrophysical sky signal, and in particular to the bright 3 mK Solar CMB dipole. Even an inaccuracy at the $\mathcal{O}(10^{-3})$ level can therefore excite correlated large scale artifacts at the microkelvin level, which are comparable to (or larger than) the expected cosmological reionization of about 0.5.

It is also important to note that even though the error bars reported in Table <ref> are of similar order of magnitude, the detailed manners in which these uncertainties are propagated into higher-level astrophysical results are very different in the two pipelines. Specifically, while the  sampling approach accounts for all couplings between the specific value of $x_{\mathrm{im}}$ and all other parameters (gain, baselines, correlated noise, CMB dipole, large scale polarization, etc.) at every single step of the Markov chain, the  approach only marginalizes over two linear templates in the low-resolution covariance matrix. These two templates are derived by changing $x_{\mathrm{im}}$ for one diode pair in a DA by +10 % and the other diode pair by $\pm10$ % with respect to their mean values. This linear low resolution approach can obviously only capture a limited subvolume of the full nonlinear effect of transmission imbalance uncertainties. Even cases for which the mean estimates formally agree within $1\,\sigma$ may therefore in practice result in significantly different sky maps and error estimates. We will return to the impact of transmission imbalance in Sect. <ref>.

Summary of noise properties.

to 2.cm#2em
3Sensitivity, $\sigma_0$ (mK $\sqrt{\mathrm{s}}$)
3Knee frequency, $f_{\mathrm{knee}}$ (mHz)     -3pt
3   3   3pt
Radiometer     Diode     GSFC     WMAP     CG/√(2)     GSFC     WMAP     CG/√(2)     CG Slope, α3pt5pt
11       1       0.72       0.66        0.704 ±0.002       6.13       0.4         0.82 ±0.20        -1.01 ±0.10      2                  0.708 ±0.003                   0.63 ±0.14        -0.95 ±0.10 12       1       0.87       0.75        0.796 ±0.004       5.37       0.51         0.42 ±0.19        -0.93 ±0.12      2                  0.780 ±0.005                   0.71 ±0.15        -1.02 ±0.10 11       1       0.75       0.71        0.788 ±0.001       1.66       0.71         1.20 ±0.22        -1.02 ±0.09      2                  0.777 ±0.001                   1.19 ±0.22        -1.02 ±0.09 12       1       0.77       0.72        0.788 ±0.003       1.29       0.32         0.62 ±0.16        -0.99 ±0.11      2                  0.784 ±0.001                   0.63 ±0.13        -1.01 ±0.11 11       1       0.99       0.92        0.998 ±0.002       3.21       1.09         1.06 ±0.16        -1.09 ±0.09      2                  0.992 ±0.002                   1.06 ±0.16        -1.10 ±0.09 12       1       0.95       1.02        1.159 ±0.007       3.13       0.35         0.45 ±0.47        -0.98 ±0.11      2                  1.146 ±0.007                   0.83 ±0.14        -1.00 ±0.09 21       1       0.89       0.85        0.908 ±0.002       1.92       5.76         2.88 ±0.37        -1.10 ±0.07      2                  0.906 ±0.002                   3.22 ±0.56        -1.10 ±0.06 22       1       1.04       0.99        1.074 ±0.004       4.61       8.62         3.95 ±0.54        -1.11 ±0.06      2                  1.064 ±0.003                   4.05 ±0.64        -1.11 ±0.06 11       1       1.25       1.22        1.551 ±0.003       2.56       0.09         1.27 ±0.15        -0.90 ±0.06      2                  1.539 ±0.003                   1.19 ±0.14        -0.89 ±0.06 12       1       1.07       1.11        1.398 ±0.002       4.49       1.41         2.11 ±0.20        -0.97 ±0.05      2                  1.432 ±0.002                   1.88 ±0.17        -0.96 ±0.05 21       1       1.01       0.97        1.241 ±0.298       2.43       0.88         1.50 ±0.24        -0.95 ±0.07      2                  1.217 ±0.294                   1.60 ±0.26        -0.97 ±0.06 22       1       1.13       1.1        1.443 ±0.300       3.06       8.35         4.01 ±0.85        -1.00 ±0.08      2                  1.415 ±0.316                   3.08 ±0.65        -1.01 ±0.08 11       1       1.18       1.35        1.938 ±0.005       16.2       7.88         5.59 ±0.53        -0.94 ±0.05      2                  1.895 ±0.005                   8.99 ±0.85        -0.95 ±0.04 12       1       1.41       1.61        2.301 ±0.005       15.1       0.66         3.91 ±0.42        -0.89 ±0.05      2                  2.345 ±0.006                   4.81 ±0.53        -0.89 ±0.05 21       1       1.38       1.61        2.225 ±0.007       1.76       9.02        13.57 ±1.47        -0.89 ±0.03      2                  2.292 ±0.006                   5.06 ±0.95        -0.93 ±0.05 22       1       1.44       1.72        2.291 ±0.006       0.77       7.47         3.02 ±0.53        -0.98 ±0.05      2                  2.232 ±0.007                   7.26 ±1.05        -0.95 ±0.04 31       1       1.47       1.65        2.328 ±0.005       1.84       0.93         1.30 ±0.46        -0.99 ±0.07      2                  2.322 ±0.006                   1.97 ±0.28        -0.98 ±0.06 32       1       1.69       1.86        2.707 ±0.015       2.39       0.28         1.59 ±0.29        -0.98 ±0.07      2                  2.579 ±0.015                   1.40 ±0.39        -1.00 ±0.07 41       1       1.6       1.71        2.519 ±0.010       8.46       46.5        26.81 ±1.83        -0.92 ±0.04      2                  2.479 ±0.009                  24.75 ±1.63        -0.92 ±0.04 42       1       1.43       1.65        2.221 ±0.017       5.31       26.0        16.10 ±1.09        -0.94 ±0.04      2                  2.202 ±0.015                  17.11 ±1.19        -0.94 ±0.04 5pt5pt

§.§ Instrumental noise

Next, we consider the instrumental noise parameters, $\boldsymbol\xi_n=\{\sigma_0, \fknee, \alpha\}$. In this case, we recall three major differences between the  and  analysis. First, while we model the noise explicitly with a $1/f$ noise profile in Fourier domain, the  analysis adopts a model independent approach by simply measuring the autocorrelation function directly. A notable advantage of the latter approach is that it naturally accounts for the non-white noise at high frequency without algorithmic modifications, while this has to be added manually in the parametric  approach. A second difference is the fact that while  uses 1- or 24-hour segments to estimate the noise model, we use 3–5 days, and are therefore able to trace noise correlations to longer timescales. Thirdly, while  assumed the noise filter to be constant within each year of operation, we allow it to vary between scans, that is, on a timescale of days.

With these differences in mind, Figs. <ref>–<ref> provides a complete overview of the noise parameters for all 40  diodes. As in Fig. <ref>, the solid black lines show  results, while the dotted red and orange lines show the corresponding 1-year and GSFC measurements, where available. Starting with the white noise level, we see that these are overall relatively constant in time, although with slight traces for annual variations in some channels (e.g., 113); slight instabilities near the beginning and/or end of the mission in other channels (e.g., ); and slight drifts in yet others (e.g., 12 and 32).

When comparing the  values with the  values, it is worth noting that  only published results for each diode-pair, not for individual diodes. All  values are therefore the same for each diode pair. Still, from the  results, which are reported individually for each diode, we see that diode pairs generally have quite similar white noise levels and vary at most by a percent.

To facilitate a more quantitative comparison, Table <ref> compares the  posterior mean results with the reported  results. Note that for $\sigma_0$, the  values have been scaled by a factor of $\sqrt{2}$, in order to account for the fact that these apply to individual diodes, as opposed to diode-pairs. Both in Table <ref> and Fig. <ref>, we see that about half of the  values lie between the two  results, while the other half are higher. In particular the -band noise levels are much higher in the  solution, sometimes by as much as 50 %.

In this respect, it is worth recalling from Sect. <ref> that the white noise level in raw du is in  not strictly sampled from the full posterior distribution, but rather estimated deterministically from the highest frequencies. This makes our estimate more sensitive to possible colored noise at high frequencies <cit.>. At the same time, the calibrated white noise level $\sigma_0[\mathrm K]=\sigma_0[\mathrm{du}]/g$ depends on the gain, and this  allows us to test the effects of the calibration on the instrument sensitivity itself. The calibrated white noise level follows a biannual trend indicative of a system temperature variation, which is to be expected given the radiometer equation

\begin{equation}
\sigma_0[\mathrm V]\propto gT_\mathrm{sys}.
\end{equation}

Aside from an overall amplitude shift due to the absolute calibration variation, the shape of the white noise level is stable throughout the Gibbs chain.

Another issue worth pointing out is the fact that we are not yet accounting for correlations between the white noise in diode pairs. However, the correlation coefficient between residuals is relatively small, with e.g., values of roughly 5 % for -band. Diode pairs in  and  have higher correlation of $\sim25\,\%$, but have similar orders of correlation with diodes from the other pair, indicating that the correlation is driven by unsubtracted sky signal.

In summary, we have not yet been able to identify the cause of the major difference in reported white noise levels in the -band; while we do detect goodness of fit failures of as much as 5–10 $\sigma$ for many of these diodes at the TOD level (see Sect. <ref>), such significances correspond to sub-percent errors in the white noise level. For comparison, a white noise misestimation of 50 % would translate into an 800 $\sigma$ $\chi^2$ failure. This is left to be understood through future work, but we do not expect it to indicate a real failure in either analysis, but it is more likely just a matter of different conventions.

Turning our attention to the low frequency parameters, we see in Table <ref> and Fig. <ref> that our knee frequencies lie between the  ground and laboratory measurements, almost without exception, which on the one hand indicates generally good agreement between the two analyses. However, our values are in fact closer to the  laboratory measurements than the  flight measurements. This may be due to the longer time-scales used in the  analysis.

PSD for diode 413 that spans MJDs 52252.3–52254.8. The power spectrum of the blue line corresponds to the residual, while the gray line is the residual with a correlated noise realization removed.

PSD for diode 413 that spans MJDs 52252.3–52254.8.
The black dashed line is a sample of the theoretical PSD, while the blue line is the smoothed residual power spectrum.

Most diodes have constant $f_\mathrm{knee}$ throughout the mission, with a few notable exceptions. First, all -band channels display some amount of temporal variation that does not seem to be associated with any sinusoidal features. Second, all 2 channels, 223, and 224 all display a similar asymptotic drift in time. We have not found any instrumental effects that share this feature.
The PSD slope $\alpha$ is around $-1$ for each diode, albeit with high scatter for the lower frequencies. As expected, the uncertainty in $\alpha$ decreases as $f_\mathrm{knee}$ increases, since there are more datapoints to fit below $f_\mathrm{knee}$ where the constraining power on $\alpha$ is the strongest.

(top): Difference between the  $\boldsymbol d_\mathrm{cal}=\boldsymbol d/g-\boldsymbol b - \boldsymbol s_\mathrm{sl}$ and the delivered calibrated TOD from . (Bottom): Raw correlated noise (gray) and smoothed data with Gaussian kernel (black). This shows the hourly baseline subtraction from the  treatment.

For completeness, Fig. <ref> shows a summary of the reduced
normalized $\chi^2$ for all diodes. The most striking features in
these figures are the amplitude and semiannual periodicity.  Given
the noise model and data residual, we can evaluate the goodness of fit
in the form of the relative $\chi^2$. Here, we find that approximately
half of the diodes have a $\chi^2$ value at least $6~\sigma$ above
or below the expected value.  Given perfect Gaussian residuals, we
would expect these values to be within $\pm1~\sigma$ 68 % of the time. For
a typical -band scan of length $n_\mathrm{TOD}=2^{22}$, a $10~\sigma$
model failure corresponds to $\chi^2/n_\mathrm{TOD}=1.003$. It is
therefore exceedingly difficult to look at any given  scan in
the time domain and identify a model failure. To illustrate this,
Fig. <ref> compares the observed noise PSD with the
best-fit model for the 413 diode. This is a $7~\sigma$ outlier;
despite this, the $1/f$ model appears to perform exceedingly well over
seven decades in frequency.

TOD corrections for -band for a single Gibbs sample, projected into maps. Columns show Stokes $T$, $Q$, and $U$ parameters. Rows show, from top to bottom, 1) correlated noise; 2) the orbital dipole; 3); bandpass mismatch leakage; and 4) sidelobe corrections. The bottom row shows the residual obtained when binning the sky and systematics-subtracted TOD into a sky map. Note that the correlated noise and residual have been smoothed by a $2^\circ$ Gaussian beam.

Only with aggressive smoothing does the model failure become apparent
at frequencies 1–10 Hz. This is illustrated in
Fig. <ref>, which shows exactly the same underlying
data as in Fig. <ref>, but heavily smoothed. Here, it is
clear that despite fitting the data well at the highest and lowest
frequencies, it is in the intermediate range of 1–5 Hz where the
$1/f$ model is a less accurate fit to the residual power
spectrum. Part of the cause of this failure is that the white noise
level is fixed by the value of the power spectrum at the
Nyquist frequency, as it was computed by differencing adjacent
samples. The power spectrum has a downward trend above 1 Hz,
indicating that the data would be better fit by a polynomial in powers of
$f^\alpha$. This is phenomenologically similar to the
collaboration's approach of describing the time-space
autocorrelation as a cubic polynomial in $\log\Delta t$

In practice, the failure of the $1/f$ model has a small effect on the final data products, and
was not visible in noise models when we modeled the data in one day scans
rather than the longer 3–7 day scans due to the lower $n_\mathrm{TOD}$ giving
a higher uncertainty on the relative $\chi^2$.  Therefore, although this
strictly constitutes a deficiency in the model, it is in practice too small to
affect the results of the rest of the chain. The downturn of the noise PSD at
high frequencies is also present in, e.g., the  HFI data
<cit.>, so improved modeling of this form will be a
necessity in future  endeavors, and will be used to improve the
data processing.

Pseudo-spectrum standard deviation for each instrumental
systematic correction shown in
Fig. <ref> (thin
colored lines). For comparison, thick black lines show spectra
for the mean of the full frequency map; thick red lines show their standard deviations (i.e., the full systematic uncertainty); and gray lines show white noise realizations.
Columns show results for ,
, 1, and 4, respectively, while rows show results for each of
the six polarization states ($TT$, $EE$, $BB$, $TE$, $TB$, and
$EB$). All spectra have been derived outside the CMB confidence
mask presented by <cit.> using the HEALPix
utility, correcting only for sky fraction and not for mask mode

Before concluding this section, we recall the close relationship between the correlated noise component and the baseline model. This is illustrated in Fig. <ref>, which shows the difference between the calibrated  and  TOD data for 113, i.e., the same as Fig. <ref>, but plotted for 24 hours instead of 10 minutes. The bottom panel shows the  correlated noise realization for the same period, both raw and smoothed. The most prominent feature in this figure is a varying signal of amplitude $0.2\,\mathrm{mK}$. This is due to the hourly baseline subtraction mentioned above, which contrasts with the  approach of assigning a linear baseline solution for the entire scan, and then accounting for the nonlinearity through $\n^{\mathrm{corr}}$. The variations are commensurate with the correlated noise correlation length, which for 113 has ${f_\mathrm{knee}\sim0.5\,\mathrm{mHz}}$, corresponding to a little over half an hour. Therefore, the hour-long baseline subtraction acts as a destriper, removing an estimate of the correlated noise.

§.§ Instrumental corrections in map domain

Returning to the global parametric data model in Eq. (<ref>), it is useful for intuition purposes to project each of the  TOD-level correction terms into sky maps. This is done for -band in Fig. <ref>, in the same form as was done for LFI by <cit.>. Columns correspond to Stokes $T$, $Q$, and $U$ parameters, while rows show, from top to bottom, 1) correlated noise; 2) the orbital CMB dipole; 3) bandpass leakage; and 4) sidelobe corrections. The bottom row shows the residual obtained after subtracting all modeled terms from the raw TOD. All maps correspond to one single and randomly selected Gibbs sample.

Starting with the correlated noise in the top row and the residual in the bottom row, we note that these are the only terms that are fundamentally stochastic in nature; all the instrumental terms are either primarily deterministic in nature, as they rely only on the sky model coupled to a small number of instrumental parameters, such as the gain, bandpass, or beam, whereas the random realizations of the sky model are highly constrained by the multifrequency component separation. As such, the correlated noise and residual maps act as the “trash cans” of Bayesian CMB analysis; together they show anything in the TOD that is not explicitly modeled. For the -band channel, we see that the correlated noise is limited to less than 1  over almost the full sky, while the residual appears like random noise over most of the sky. The main exceptions to this are strong residuals near the Galactic plane in temperature, which indicates the presence of unmodeled foreground features; typical candidates to explain this would be curvature in the synchrotron spectral index or a more complicated AME SED than that assumed here. Secondly, at high latitudes we see traces of a small dipole with an amplitude of 1–2 aligned with the CMB dipole. This indicates the presence of an absolute calibration deviation of about 0.03–0.06 %; this is within the uncertainty of the absolute -band calibration prior of 0.1 % discussed in Sect. <ref>, and when inspecting different Markov chain samples of the same type, one can see that the amplitude and sign of this residual scatter is around zero as expected.

Next, the amplitude of the orbital dipole is about 250 in temperature and 2.5 in polarization. This component by itself is exceedingly well known, as it depends only on the satellite velocity and the CMB monopole temperature. However, when actually fitting this to the real data, it obviously also depends on both the gain, and sample-to-sample variations in this map, and  is therefore a good tracer of gain uncertainties. In addition, its physical appearance also in principle depends on the sidelobe model, but we do not yet propagate sidelobe uncertainties anywhere in the analysis.

The third row shows the corrections for bandpass leakage.[In principle, this also accounts for leakage due to beam differences, as discussed by <cit.>, but for  this effect is much smaller than the bandpass effect.] As for  LFI, this term is by far the largest correction in polarization, with an amplitude that is almost an order of magnitude larger than any of the others. It is highly nontrivial to compare this component to the  analysis, since this effect was accounted for by solving for a spurious $S$ map as part of their mapmaking. However, it is important to note that whether one models this effect explicitly, as we do, or project it out during mapmaking, as the  pipeline did, the accuracy of the bandpass corrections depends directly on the accuracy of the gain and transmission imbalance calibration.

The fourth row shows the sidelobe contribution. Here we see that the temperature amplitude reaches 50 , which corresponds to 1.5 % of the CMB Solar dipole amplitude. If it should turn out that the sidelobe model were incorrect by, say, 30 %, this translates directly into an error in the absolute -band calibration of about 0.5 %, which is significantly greater than the statistical error shown above. Given the degeneracies discussed in Sect. <ref>, there is thus also a strong degeneracy between the AME dipole and the -band sidelobe. For polarization, the amplitude is mostly smaller than 1, and therefore of minor importance for this channel. We note that only the intensity component of the  sidelobe model has been published, and polarized sidelobes are not accounted for in the current processing. However, as shown by <cit.>'s Table 2, the amplitude of this contribution is relatively small, with a high-latitude mean of $0.8\,\mathrm{\mu K}$ for -band and $<0.1\,\mathrm{\mu K}$ for all other bands.

§.§ Instrumental uncertainties in power spectrum domain

We conclude this section by estimating the uncertainty of each instrumental effect in terms of angular power spectra, using the same methodology as <cit.> for  LFI. That is, we first compute the power spectra for each individual instrumental correction Markov chain sample, as illustrated in Fig. <ref>, and then compute the standard deviation of all these samples. The results from this computation are summarized in Fig. <ref> for four DAs, namely , , 1, and 4. In each panel the black line shows the power spectrum of the full coadded frequency sky map (including both signal and noise), while the gray line shows the white noise level. The thick dark red line shows the sum of all variations, while the thin colored lines show the contribution from individual correction terms.

On large angular scales in the $TT$ spectrum, we see that the dominant uncertainty comes from the orbital dipole (blue lines), which trace gain uncertainties. This is reasonable, since these data are strongly signal-dominated; indeed, for -band even the sidelobe contribution (green lines) is higher than the correlated noise.

For large-scale $E$-mode polarization, the dominant term varies from channel to channel. Specifically, for -band the bandpass leakage (thin red lines) and gain fluctuations are significantly larger than the correlated noise, while for 1 and -band the correlated noise dominates for most multipoles, although the gain fluctuations are comparable for some $\ell$'s.

An important conclusion to be drawn from these measurements is that a simple uncertainty model that primarily accounts for correlated noise is likely to be suboptimal for detailed cosmological analysis of large-scale polarization. Both gain and bandpass uncertainties are at least as important for the lowest multipoles, and simultaneously accounting for all of these contributions is important in order to derive robust cosmological results.

§ FREQUENCY MAPS

In this section, we discuss the reprocessed  frequency maps and their properties. In Sect. <ref> we present the reprocessed  maps themselves, commenting on notable features in the maps themselves, as well as differences with the  results. In Sect. <ref> we focus on the consistency of our new maps, both internally among the  channels, and with respect to . Section <ref> describes the efficiency of template-based transmission imbalance uncertainty propagation.

§.§ Map survey

Posterior mean -band map produced with the  pipeline. Rows show Stokes $T$, $Q$, and $U$, respectively. The temperature map is shown at full resolution, while the polarization maps are smoothed with a $2^{\circ}$ FWHM Gaussian beam to reduce small-scale noise.

Same as Fig. <ref>, but for -band.

Same as Fig. <ref>, but for -band.

Same as Fig. <ref>, but for -band.

Same as Fig. <ref>, but for -band.

Difference maps between the  and 9-year  frequency maps. Columns show Stokes $T$, $Q$, and $U$ parameter maps, while rows show -, -, -, -, and -band maps. The maps are all smoothed with a $2^\circ$ FWHM Gaussian beam.

Posterior variation maps for -band. Columns show the Stokes parameters and the correlation coefficient between $Q$ and $U$, while the rows show (top) the white noise rms per pixel and (bottom) the posterior standard deviation. The rms maps are unsmoothed, while the standard deviations have been smoothed to $7^\circ$.

Difference between two -band Gibbs samples, smoothed to $7^\circ$.

TOD-level residual map for -band, smoothed with a $5^\circ$ FWHM Gaussian beam.

We start by showing the co-added frequency -band, -band, -band, -band, and -band posterior mean maps in Figs. <ref>–<ref>, all defined in thermodynamic  $\mu\mathrm{K_{CMB}}$ units. All maps are produced at the DA level, and in these figures the , , and -band maps are generated by inverse-variance weighting the individual DAs. The temperature maps are presented at full angular resolution, while the polarization maps have been smoothed with a $2^\circ$ Gaussian beam to reduce small scale structure and make the larger scale effects more apparent. Overall, the temperature maps behave as expected from the official  analysis, with falling foreground amplitudes with frequency. Furthermore, it is very difficult indeed, if not impossible, to see visual differences between the  and  maps by eye when switching rapidly between them. However, it is important to note that the  frequency maps retain the Solar CMB dipole, following <cit.> and <cit.>, while it is removed in the  official maps. Similarly, we see that the amplitudes of the polarized maps decrease as expected from –-band following the expected synchrotron behavior, with a slight increase at -band due to the contribution of thermal dust.

Next, in Fig. <ref>, we show corresponding difference maps between the official 9-year  maps and the maps produced in this work. The color scale in this plot is linear with range $\pm5\muK$, and we see that the differences are thus quite small, typically smaller than 1–2 for most channels. The main exception to this is -band polarization, for which the differences are generally larger than $5\muK$.

Going into greater detail and starting with total intensity, we see first that the -band difference is dominated by a dipole with a $\sim$ 2.5 amplitude that is anti-aligned with the CMB Solar dipole. In addition, the Galactic plane is slightly brighter in  than . Both of these suggest that our total absolute -band calibration is lower than the  value by about 0.1 %; given the major differences in methodology described in Sect. <ref>, this degree of agreement is a major validation of the data and the two  pipelines. A similar small dipole difference is also seen in the -band.

For the remaining channels, and in particular for the - and -bands, the main intensity difference takes the form of a quadrupole with an amplitude of 2–3 aligned with the CMB dipole. Naively, one could suspect this to be due to different treatments of the relativistic quadrupole. However, as noted by <cit.>, the  maps retain the kinematic quadrupole, as does ; in our framework, this signal term is accounted for through the signal model defined in Eq. (<ref>). This is notably different from the  maps, which do remove the relativistic quadrupole from the frequency maps <cit.>. Additionally, even though the observed quadrupole has the expected shape, the frequency dependence is not consistent with the expected functional form $x\coth x$ where $x=h\nu/(2kT_\mathrm{CMB})$ <cit.>. For now, we speculate that these differences are rather due to second-order gain or baseline differences, possibly associated with the annual oscillatory structures seen in Fig. <ref>.

In polarization, we note large scale differences in both Stokes $Q$ and $U$. These differences do not match known Galactic component morphologies, but are more reminiscent of the poorly measured transmission imbalance modes discussed by <cit.>, although the map-space morphologies are not identical. In general, such large mode differences are due to at least three main effects: 1) incomplete polarization angle coverage for a few large-scale modes; 2) errors in transmission imbalance coupled with the Solar dipole; and 3) interplay between the transmission imbalance, the far sidelobe, and the Solar dipole, as briefly described in Sect. <ref>. The scale of these effects is most pronounced in the -band polarization results, where we see the largest differences between the two processing pipelines.

From these differences alone, it is not possible to determine whether the excess structures are present in the  or  maps, or both. However, Appendix <ref> provides a complete survey of the  frequency maps, and in particular Fig. <ref> compares these with the  maps. In this case, one clearly sees that the large-scale modes are predominantly present in the  maps, rather than .

Returning to the internal properties of the  posterior distribution, we show in the top panel of Fig. <ref> the -band white noise standard deviation per pixel in Stokes $T$, $Q$ and $U$; the fourth column shows the correlation coefficient between the $Q$ and $U$ coefficients. The bottom panel shows the corresponding posterior standard deviation per pixel. It is important to note that the white noise is not a free parameter in the data model, and there is no white noise component in the Gibbs sampler described by Eqs. (<ref>)–(<ref>). That also implies that there is no marginalization over white noise in the resulting frequency map ensemble. Rather, the full marginal uncertainty per frequency map pixel must be obtained by adding the two rows in Fig. <ref> in quadrature. However, a preferable approach to performing error propagation for higher level scientific analyses is to analyze each sample separately, taking into account only white noise for each sample, and then use the full sample ensemble as the final result. This is the only robust way of fully accounting for all correlations between the various effects.

The white noise pattern for $T$ follows the usual pattern with highest sensitivity at the North and South ecliptic poles, as well as circles around the poles corresponding to times when the partner horn is observing the opposite ecliptic pole. There are also regions of higher noise level corresponding to planets crossing the ecliptic, and regions of higher emission $\simeq140^\circ$ away from the Galactic center, which correspond to the times when the partner horn lies within the processing mask.

Comparison of the $C_\ell^\mathrm{TT}$, $C_\ell^\mathrm{EE}$, and $C_\ell^\mathrm{BB}$ from  and . Each row corresponds to a different DA, with frequency increasing from top to bottom. (left): ratio of $C_\ell^\mathrm{TT}$ from  compared to . (middle/right): $C_\ell^\mathrm{EE/BB}$ power spectra with  in blue and  in orange.

The polarized RMS maps share all of these characteristics, but with an overall amplitude shift due to polarization measurements having half the effective number of observations per pixel. In addition, the poles have a characteristic “X”-like structure that is rotated $45^\circ$ degrees between $Q$ and $U$, corresponding to different polarization orientations. There are also characteristic large scale structures visible in Galactic coordinates, corresponding to polarization modes poorly constrained by the  scan strategy.

While the maps in the top row of Fig. <ref> are directly comparable to the corresponding  products, the posterior standard deviation shown in the bottom row has no equivalent in the official  release. These maps can be considered the “systematic” error contributions, as their variation depends on the sampled instrumental parameters, i.e., gain, imbalance parameters, correlated noise, and sidelobe correction. The temperature map contains a clear quadrupole signature. This is due to the variation in the absolute calibration $g_0$, which changes the Solar dipole in the final map. In addition to the quadrupole, the Galactic plane also varies due to the gain solution's fluctuations. As expected, the white noise patterns associated with the scan strategy also appear in the polarization maps, which have much lower signal-to-noise ratio than the temperature map.

Another useful quantity is the difference between two arbitrary samples, which we show in Fig. <ref>. In temperature, the most striking term is a dipole, corresponding to the absolute gain difference, and the Galactic plane. There are also additional weaker lines associated with the scanning strategy, corresponding to different correlated noise and time variable gain realizations. In polarization, gain variations,  bandpass uncertainties, and correlated noise dominate the differences between two samples, as quantified in Fig. <ref>. The polarization differences are aligned with 's scans, modulated by the polarization angle.

Finally, the quality of the model in map space can be evaluated quite well by looking at the calibrated residual map, i.e., mapping the time-ordered residual $\boldsymbol r\equiv(\data-\ncorr)/\boldsymbol g-\boldsymbol s^\mathrm{tot}$. We display this TOD residual for the -band in the bottom panel of Fig. <ref> ($2^\circ$ FWHM) and Fig. <ref> ($5^\circ$ FWHM). In temperature there is a large residual along the Galactic plane, which is to be expected for both temperature and polarization, due to the complexity of the Galactic center. The temperature residual also contains a $\sim3\,\mathrm{\mu K}$ dipole due to the prior sampling of $g_0$ – drawing a sample based on the sky model would track the dipole much more closely, whereas the prior sampling by definition does not use the sky model to draw $g_0$. Other than the Solar dipole, Galactic plane, and point sources, both the temperature and polarization maps are visually consistent with white noise across the entire sky.

Again, we have for brevity primarily focused on the -band in this discussion. For completeness, however, similar plots for all DAs are shown in Appendix <ref>. In particular, Fig. <ref> compares the  and  DA polarization maps, Figs. <ref> and <ref> shows the white noise and posterior rms's, Fig. <ref> shows sample differences, and Fig. <ref> shows TOD residual maps.

In Fig. <ref> we compare angular power spectra computed from both   and  frequency maps. These spectra are derived using the  <cit.>[<https://github.com/LSSTDESC/NaMaster>]  routine, while applying the extended  temperature analysis mask which allows a sky fraction of 68.8 %. As the $TT$ power spectrum is strongly signal-dominated for $\ell\lesssim200$ for all DAs, it is particularly informative to consider ratios, which we show in the left column of Fig. <ref>. Here we see that the $TT$ spectra derived from the two pipelines are consistent to sub-percent level at all but the very largest and smallest scales for all DAs. We speculate that the large scale differences are due to different CMB Solar dipoles – as noted above, the  maps retain the Solar CMB dipole, and an estimate of this must be subtracted before evaluating these spectra. In contrast, the  maps have this contribution removed at the TOD level; small differences due to these different treatments are not unexpected. The small scale differences above $\ell\sim200$ can be attributed to different data selections and low-level processing. For instance, the  maps exploit about 1 % less data than ;  fits one $\sigma_0$ parameter per scan, while  assumes it to be constant for each year; the  gain model varies every 23 sec, while the  model assumes constant gain per scan etc.

Internal  difference maps, smoothed by $10^\circ$. The two left columns are Stokes $Q$, and the two right columns are Stokes $U$, with the  and  maps alternating between columns. The top to bottom rows are difference maps in increasing frequency.

Difference maps between similar  and  frequency maps. The comparison plots go, by column: Stokes $Q$ for the -produced  and  sky maps, Stokes $Q$ for official  and  data products, Stokes $U$ for the  sky maps, and Stokes $U$ for the official data products. (Top row):  LFI 30 GHz minus -band, scaled by the synchrotron power-law. (Top middle row):  -band minus LFI 30 GHz, also scaled by the synchrotron power-law. (Middle row):  -band compared to the LFI 44 GHz sky maps, scaled by the synchrotron power-law. (Bottom middle row):  -band minus LFI 70 GHz, with unit scalings for each band. (Bottom row): The  DR4 100 GHz map minus the  -band  also with unit scalings for each band.

The $E$-mode power spectra, displayed in the second column of Fig. <ref>, are mainly dominated by noise and polarized synchrotron emission. As expected, the large scale foreground-dominated multipoles decrease in amplitude according to the relative amplitude of the synchrotron spectrum. Overall, the  and  power spectra appear fairly consistent for the – channels, while at - and -band there is noticeably more scatter at low multipoles in the  spectra than in the  spectra.

The $B$-mode power spectra, displayed in the third column of Fig. <ref>, are expected to follow a similar pattern, but since foregrounds are generally reduced by a factor of $\simeq2$–4 <cit.>, this spectrum is less signal-dominated, and therefore more susceptible to instrumental systematics. For instance, the $C_{\ell=3}^\mathrm{BB}$ mode has been identified as being particularly poorly constrained due to its symmetry aligning with $\gtrsim10\,\mathrm{min}$ signals in the TOD induced by the  scan strategy <cit.>. In this figure, it appears that these low-$\ell$ modes appear significantly better constrained in the  maps than in  for - and -bands, and the overall large scale noise level is lower by one or two orders of magnitude.

Full sky half-difference spectra. The red lines are the power spectra of the  difference maps, while the black lines are the same for the reprocessed  maps.

§.§ Consistency tests through inter-channel difference maps

As described in Sect. <ref>, the - and -bands each had two DAs, while the -band had four DAs, and computing differences between the corresponding DA maps can highlight mismodeled systematics. While the -band and -band have different central frequencies, they are close enough that we can compare them by scaling -band assuming a polarized synchrotron power law SED of $\beta_\mathrm s=-3.1$. Similarly, internal differences between scaled , , and LFI 30 GHz maps provide an important null test. In particular, the –30 difference has received significant attention ever since the  2015 data release <cit.>, showing clear signatures of instrumental systematics. These were gradually reduced through improved  processing in the  2018 <cit.>, PR4 <cit.>, and  <cit.> data releases. Still, even after all these developments, large-scale residuals remained that were difficult to resolve through further LFI improvements <cit.>. In this section, we revisit this question for the  products.

We start by inspecting internal  half-difference maps of the form $(Q1-Q2)/2$ etc. These are plotted in Fig. <ref>. Here we see that the -band and -band half-difference maps from  have virtually no trace of poorly measured modes, and the differences appear to be well-traced by the rms maps. In contrast, the  half-difference maps show clear evidence of large-scale residuals. The largest visual improvement is in the -band, where the  case is almost entirely consistent with instrumental noise, as opposed to the  difference that is dominated by large-scale residuals.

Next, Fig. <ref> shows comparisons between the  $K$-
and $Ka$-bands and the LFI 30 GHz channel, between the  $Q$-band and LFI
44 GHz, and finally between  $V$-band and LFI 70 GHz. When comparing
maps with  LFI, we use  products, which represent the
cleanest version of  LFI published to date. For the  map
comparison, both  and  maps were produced by this joint analysis.
Additionally, we compare the mean -band maps with the  HFI DR4
100 GHz channel. It is worth noting that this 100 GHz map has had no input
from  so this difference map is an independent comparison
between two datasets and processing methods.

Starting with the  maps, we see in the first and third columns of
Fig. <ref> that the magnitude of the differences are small
in both Stokes $Q$ and $U$. Overall, across all five frequency map comparisons
we see small levels of variation, with structure contained to the Galactic
plane. Notably, however, there is a larger sky signal within the
$\mathit{Ka}-30$ Stokes $Q$ comparison.
This large-scale difference also exists in the $\mathit Q-44$ Stokes
$Q$, but it did not appear in the internal  half-difference map.

Columns two and four of Fig. <ref> show corresponding differences
between the official  and  LFI frequency maps. Similar to the
sky map comparisons, we see differences in the Galactic center,
and to a lesser degree along the Galactic plane due to the slight differences
in the frequency coverage. When comparing the official  maps,
particularly for -band, we see structures sweeping across large angular
scales across the sky, likely due to the poorly measured modes in -band.

Of particular note is the $100-\mathit W$ difference map. The difference maps here have a similar level of white noise and Galactic
contamination as the $\mathit V-70$ maps,  whereas the  differences
are driven by obvious transmission imbalance modes, each with an opposite sign and
magnitude. The difference between 100 GHz and  demonstrates that the good
agreement between the  and  LFI is not simply due to fitting
low-level parameters in a joint analysis framework – by obtaining -band maps
that are consistent with an independent 100 GHz polarization map, we have
shown that the -LFI agreement is the result of a genuine improvement in data processing.

Transmission imbalance template amplitudes for each  radiometer as estimated by fitting the official templates to low-resolution difference maps between  and . The templates are provided in mK, and the template amplitudes are therefore dimensionless. The fourth column lists the fractional decrease in standard deviation, $\Delta\sigma/\sigma\equiv\sqrt{\sigma_{\mathrm{raw}}^2 - \sigma_{\mathrm{corr}}^2}/\sigma_{\mathrm{raw}}$,  after subtracting the best-fit templates.

to 1.8cm#2em
DA     a_1     a_2     Δσ/σ3pt5pt