Issue 
A&A
Volume 675, July 2023
BeyondPlanck: endtoend Bayesian analysis of Planck LFI



Article Number  A13  
Number of page(s)  27  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202243186  
Published online  28 June 2023 
BEYONDPLANCK
XIII. Intensity foreground sampling, degeneracies, and priors
^{1}
Institute of Theoretical Astrophysics, University of Oslo, Sem Sælands vei 13, 0371 Oslo, Norway
email: k.j.andersen@astro.uio.no; andersen.kristian.joten@gmail.com
^{2}
Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, Milano, Italy
^{3}
INAFIASF Milano, Via E. Bassini 15, Milano, Italy
^{4}
INFN, Sezione di Milano, Via Celoria 16, Milano, Italy
^{5}
INAFOsservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, Trieste, Italy
^{6}
Planetek Hellas, Leoforos Kifisias 44, Marousi, 151 25
Greece
^{7}
Department of Astrophysical Sciences, Princeton University, Princeton, NJ, 08544
USA
^{8}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA
^{9}
Department of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{10}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{11}
Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
^{12}
Haverford College Astronomy Department, 370 Lancaster Avenue, Haverford, PA, USA
^{13}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{14}
Dipartimento di Fisica, Università degli Studi di Trieste, Via A. Valerio 2, Trieste, Italy
Received:
25
January
2022
Accepted:
29
August
2022
We present the intensity foreground algorithms and model employed within the BEYONDPLANCK analysis framework. The BEYONDPLANCK analysis is aimed at integrating component separation and instrumental parameter sampling within a global framework, leading to complete endtoend error propagation in the Planck Low Frequency Instrument (LFI) data analysis. Given the scope of the BEYONDPLANCK analysis, a limited set of data is included in the component separation process, leading to foreground parameter degeneracies. In order to properly constrain the Galactic foreground parameters, we improve upon the previous Commander component separation implementation by adding a suite of algorithmic techniques. These algorithms are designed to improve the stability and computational efficiency for weakly constrained posterior distributions. These are: (1) joint foreground spectral parameter and amplitude sampling, building on ideas from MIRAMARE; (2) componentbased monopole determination; (3) joint spectral parameter and monopole sampling; and (4) application of informative spatial priors for component amplitude maps. We find that the only spectral parameter with a significant signaltonoise ratio using the current BEYONDPLANCK data set is the peak frequency of the anomalous microwave emission component, for which we find ν_{p} = 25.3 ± 0.5 GHz; all others must be constrained through external priors. Future works will be aimed at integrating many more data sets into this analysis, both map and timeordered based, thereby gradually eliminating the currently observed degeneracies in a controlled manner with respect to both instrumental systematic effects and astrophysical degeneracies. When this happens, the simple LFIoriented data model employed in the current work will need to be generalized to account for both a richer astrophysical model and additional instrumental effects. This work will be organized within the Open Sciencebased COSMOGLOBE community effort.
Key words: cosmic background radiation / cosmology: observations / cosmology: miscellaneous
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The cosmic microwave background (CMB) represents one of our best sources for knowledge of the early universe. The intensity of the CMB peaks within the microwave frequency range at about 161 GHz (Mather et al. 1994) and a long line of experiments have targeted this frequency range since the CMB was first discovered by Penzias & Wilson (1965). However, there are many sources of radiation in the microwave sky that obscure our view of the CMB, both from within the Milky Way Galaxy and from distant sources (see, e.g., Delabrouille et al. 2013; Planck Collaboration IV 2020, and references therein). Each of these sources must be modeled to high accuracy in order to establish a clean estimate of the CMB sky.
A modeling of the Galactic foreground emission in the microwave frequency range has been a vital venture in the CMB field as modern observations require higher accuracy to properly characterize the fluctuations in the CMB (Planck Collaboration X 2016). During the Planck mission, several different component separation software were created and compared to clean the microwave sky (Planck Collaboration IV 2020). Since the official Planck analysis ended, considerable development has occurred within the Bayesian Commander component separation software, implementing a comprehensive scheme which feeds the results from component separation back into the timedomain analysis, resulting in well defined posteriors (BeyondPlanck Collaboration 2023).
The current paper outlines the additions to the component separation portion of Commander, highlighting the new functionalities and their feasibility through a set of simulations. The test case for these features is within the context of the BEYONDPLANCK framework, which aims to reprocess the Planck Low Frequency Instrument (LFI) data in an endtoend Bayesian approach. Here, we present the lowfrequency foreground sky model and component separation algorithms used for intensity analysis within the BEYONDPLANCK framework, as applied to the Planck (Planck Collaboration I 2020) LFI (Planck Collaboration II 2020). As the focus of the BEYONDPLANCK analysis is limited to LFI, we limit the amount of ancillary data used to assist in the component separation effort, aiming for LFI driven posteriors.
BEYONDPLANCK is a novel endtoend Bayesian CMB analysis framework that builds on decades of experience gained within the Planck collaboration and its main defining feature is that instrument characterization and calibration is performed jointly with mapmaking and component separation, resulting in a single statistically consistent model for the full data set. As such, the foreground sky model plays a key role in the process, feeding directly into many aspects of the analysis, from gain and correlated noise estimation via leakage corrections and mapmaking to final component maps, CMB estimates, and cosmological parameters. An overview of the full process is provided in BeyondPlanck Collaboration (2023) and its companion papers.
As the goal of the BEYONDPLANCK project is to focus on endtoend error propagation with respect to the Planck LFI data set, a limited set of external data is concerned in order to ensure that the statistical weight of LFI drives the results. The ancillary data included here are described in Sect. 3, and come in two flavors, namely as foreground amplitude priors and pixelized frequency maps to add constraining power to spectral index parameters. In order to properly characterize the foregrounds within this limited data set, a suite of algorithms are introduced to the Commander Gibbs sampling framework. These algorithmic implementations are the focus of the current paper. These algorithms are not exclusively useful to the minimal BEYONDPLANCK data set and they are widely applicable to future analyses.
The algorithm used in this paper derives most closely from a similar Commanderbased (Eriksen et al. 2004, 2008) Bayesian analysis of the Planck, WMAP, and Haslam et al. data presented by Planck Collaboration IX (2016). The main differences between the two analyses are as follows. First, in the previous analysis there was limited feedback between the gain estimation and component separation, as only a single overall absolute calibration factor was fitted during the component separation process. In the current analysis, a full timedependent timeordered data model is propagated throughout the analysis. Second, the previous analysis provided only a single maximumlikelihood solution for each component, together with marginal perpixel uncertainties. In the present analysis, we provide a full Monte Carlo ensemble of samples drawn from the full posterior distribution, which represents full endtoend propagation of all uncertainties. Third, the computational framework used in the Planck 2015 analysis required uniform resolution across all frequency bands, thus limiting the resolution to that of the instrument with the poorest resolution. This was improved upon in the Planck 2018 analysis (Planck Collaboration IV 2020), where a new computational framework that allows for different resolution and smoothing scales was developed. However, the corresponding analysis only included a single joint lowfrequency component in intensity and did not provide new estimates of the individual lowfrequency components. In the current analysis, we provide fullresolution parameter maps for all components listed above, limited only by the signaltonoise ratio (S/N) of the data in question.
In this paper, we also describe the implementation of four important new algorithmic features into the Commander framework that are all designed to improve sampling efficiency and stability for weakly constrained posterior distributions. The first of these is a joint spectral parameter and amplitude sampler that employs the marginal spectral parameter posterior distribution to move quickly through the multidimensional parameter space. This algorithm was first introduced in the CMB literature by Stompor et al. (2009) and later implemented in the MIRAMARE component separation code by Stivoli et al. (2010). The main advantage of this approach is a significantly reduced Markov chain Monte Carlo (MCMC) correlation length and overall lower computational costs as a result. In this paper, we demonstrate this sampler on the peak frequency of the anomalous microwave emission (AME) spectral energy density (SED). However, the importance of this new step will increase significantly when additional spectral parameters are explored in the future and will be vital in probing, for instance, the thermal dust spectral index and temperature efficiently with the Planck High Frequency Instrument (HFI).
The second algorithmic improvement is componentbased monopole determination. As discussed at length by, for instance, Planck Collaboration X (2016) and Wehus et al. (2017), one of the most important challenges regarding intensitybased spectral parameter estimation is an accurate determination of the zerolevels at each frequency band; any error in this will translate directly into a bias in spectral parameters, which typically manifests itself as a spatial correlation between the spectral parameter map and the corresponding component amplitude map. At the same time, it is important to note that few (if any) CMB experiments actually have sensitivity to the true sky monopole, and they have therefore typically instead resorted to morphologybased algorithms to determine meaningful zerolevels. In this paper, we point out that this problem may be entirely circumvented by instead focusing on determining the zerolevels of the astrophysical component maps, rather than individual frequency maps. The resulting global sky model may then be used to determine the frequency map offsets. This approach is significantly more transparent from a physical point of view (e.g., the CMB temperature perturbation component may be assumed to have an identically vanishing monopole), it automatically guarantees consistency between the zero levels at different frequency channels and it gives zero statistical weight to the frequency monopoles in the actual fitting procedure.
Nevertheless, there is a formal degeneracy between the spectral parameters and the frequency monopoles at each step in the algorithm and to eliminate the Markov chain correlation length increase from this degeneracy, we additionally implement a new joint spectral parameter and monopole sampler as our third algorithmic improvement.
Finally, we generalize the concept of informative spatial component map priors that was introduced by Planck Collaboration IV (2020) and Planck Collaboration Int. LVII (2020) and use the results to apply informative physical Gaussian spatial priors. This can be leveraged to significantly reduce correlations between various sky components on small angular scales and, in particular, degeneracies between AME, freefree and CMB (Colombo et al. 2023) may be alleviated in this manner. Of course, the ideal approach to resolve such degeneracies is not through the use of informative priors, but rather by integrating additional data. The current algorithms allow for a gradual and controlled introduction of such data sets, without introducing pathological artifacts along the way.
The rest of the paper is organized as follows. Section 2 gives a short overview of the BEYONDPLANCK framework, and sky signal model. Section 3 describes the data set used within this analysis, the motivation for the data selection, and the simulation data used for algorithm validation. Section 4 describes the main sampling algorithms for intensity foreground parameters. Validations, through simulations, of the new sampling algorithms are presented in Sect. 5. The main results are summarized in Sect. 6 and we present our conclusions in Sect. 7. We note that polarizationbased component separation results are discussed separately by Svalheim et al. (2023b).
2. Overview of the BEYONDPLANCK sampling framework
The BEYONDPLANCK project has implemented an integrated endtoend data analysis pipeline for CMB experiments (BeyondPlanck Collaboration 2023), connecting all steps going from raw timeordered data to cosmological analysis in a selfconsistent Bayesian framework, finally realizing ideas originally proposed almost 20 years ago by Jewell et al. (2004) and Wandelt et al. (2004). This methodology allows us to characterize degeneracies between instrumental and astrophysical parameters in a statistically welldefined framework, with uncertainties propagating consistently through all stages of the pipeline. It also seamlessly connects lowlevel instrumental quantities like gain (Gjerløw et al. 2023) and correlated noise (Ihle et al. 2023), bandpasses (Svalheim et al. 2023a), and far sidelobes (Galloway et al. 2023b) via Galactic parameters such as the synchrotron amplitude and spectral index (current paper and Svalheim et al. 2023b), to the angular CMB power spectrum and cosmological parameters (Colombo et al. 2023; Paradiso et al. 2023). In this section, we provide a brief review of the BEYONDPLANCK sky model, data selection, and sampling scheme; we refer to the various companion papers for more details on the model.
2.1. Data model and Gibbs chain
In BEYONDPLANCK, the most basic data sets are raw uncalibrated timeordered data (TOD), which are modeled as follows:
Here j represents a radiometer label; t indicates a single time sample; p denotes a single pixel on the sky; and c represents one single astrophysical signal component. Furthermore, g denotes the instrumental gain; P denotes the pointing matrix; B^{symm} and B^{asymm} denote the symmetric and asymmetric beam matrix, respectively; a represents the astrophysical signal amplitudes; β shows the corresponding spectral parameters; Δ_{bp} are the bandpass corrections; M_{cj} denotes the bandpassdependent component mixing matrix; s^{orb} is the orbital dipole; s^{fsl} are the far sidelobe corrections; s^{1 hz} represents electronic 1 Hz spike corrections; n^{corr} is the correlated noise; and n^{w} is the white noise. This simple model has already been demonstrated to be an excellent fit to the Planck LFI measurements by Planck Collaboration II (2014), Planck Collaboration II (2016), Planck Collaboration X (2016), Planck Collaboration II (2020), although the current work stands as the first time it has been formulated in terms of one single equation. It is also important to note that the reason this model actually does work for Planck LFI is to a large extent, thanks to the fact that the LFI radiometers have very well behaved systematics; for other detector types, more complex models are very likely needed. For a more detailed explanation of the current model, we refer to BeyondPlanck Collaboration (2023) and companion papers.
Data sets may also be included in the form of preprocessed pixelized sky maps, m_{ν}, in which case the above data model is simplified to:
When producing m_{ν}, all timedependent quantities (i.e., the far sidelobe, orbital dipole, 1 Hz spike, and correlated noise contributions) in Eq. (1) are explicitly subtracted from the TOD prior to mapmaking, leaving only sky stationary contributions in the final pixelized map. However, at this level only a very limited set of instrumental parameters may be accounted for per frequency, namely an overall absolute calibration factor, g_{ν}, an azimuthally symmetric beam, B^{symm}, and white noise, . This latter expression may be written on the following compact form:
where A_{ν}(β) is an effective mixing matrix that takes into account both the frequency scaling of each component and beam convolution.
The goal of the Bayesian approach is to sample from the joint posterior distribution,
This is a large and complicated distribution, with many degeneracies. However, exploiting the Gibbs sampling algorithm (Geman & Geman 1984) we may factorize the sampling process into a finite set of simpler sampling steps. In this algorithm, samples from a multidimensional distribution are generated by sampling from all corresponding conditional distributions. The BEYONDPLANCK Gibbs chain may be written schematically as follows (BeyondPlanck Collaboration 2023),
where ← indicates sampling from the distribution on the righthand side.
We note that not all of these steps follow the strict Gibbs approach of conditioning on all other parameters. Most notably for us, this is the case in Eq. (9) for the spectral parameters sampler, P(β ∣ d, …), which places conditions on all parameters, except for a. Instead, we effectively sample a and β jointly by exploiting the definition of a conditional distribution, as detailed in Sect. 4.1. The advantage of a joint sampling step is a significantly shorter Markov correlation length as compared to standard Gibbs sampling.
A very convenient property of Gibbs sampling is its modular nature, as the various parameters are sampled within each conditional distribution, but joint dependencies are explored through the iterative scheme. In this paper, we are therefore only concerned with the sampling of two of the above steps, namely Eqs. (9) and (10). For all other sampling steps, we refer to BeyondPlanck Collaboration (2023) and references therein.
2.2. Astrophysical sky model
The dominant astrophysical foreground components in the Planck LFI frequencies are synchrotron, AME, freefree, thermal dust emission, and compact radio sources. The modeling of each of these components is detailed in BeyondPlanck Collaboration (2022), so we only review the relevant details in this paper. In Table 1, we summarize the models for each component in terms of free parameters, priors, and SEDs. In addition, each diffuse component is modeled in terms of an amplitude sky map, a, at a given reference frequency ν_{0} in brightness temperature units. Scaling to arbitrary frequencies is performed through the SEDs, such that the actual observed signal at a given frequency, ν, may generally be written as:
Summary of main parametric signal models for the temperature analysis.
where i denotes the specific component, ν_{0, i} is the reference frequency of the given component, β_{i} is a set of componentspecific spectral parameters, and f_{i} is the SED. For diffuse components, a_{i} is defined in terms of spherical harmonic space with a maximum multipole, ℓ_{max} defined for each component, depending on the signaltonoise ratio (S/N) and angular resolution of the data sets supporting that component. For instance, synchrotron and AME have lower values of ℓ_{max} than the thermal dust and CMB. In addition, as discussed in Sect. 4.4, we regularize the highℓ multipoles of each component with some smoothing prior, either derived from the known physical behaviour of the respective component (e.g., Planck Collaboration X 2016) or by a Gaussian smoothing operator. For compact sources, a_{i} represents simply the flux density in mJy, with the spectral index α also defined in mJy, and an explicit unit conversion factor, U_{mJy}, converts from flux density to brightness temperature units.
With this notation, the astrophysical sky model used for the current BEYONDPLANCK analysis may be written as follows:
where a_{CMB} and a_{quad} are given in thermodynamic temperature units (K_{CMB}), a_{j, src} in flux density units (mJy), and all other amplitudes a_{i} are given in terms of brightness temperature (K_{RJ}). The amplitude of component i is equal to that observed at a monochromatic frequency, ν_{0, i}. The sum in Eq. (18) runs over all compact sources brighter than some flux threshold as defined by an external source catalogue. In particular, we adopt the same catalogue as Planck Collaboration IV (2020), which is a hybrid of the AT20G (Murphy et al. 2010), GB6 (Gregory et al. 1996), NVSS (Condon et al. 1998) and PCCS2 (Planck Collaboration XXVI 2016) catalogs, comprising a total of 12 192 individual sources.
When comparing the results from the above model with previous work, it is important to note that we fit a straight powerlaw for the synchrotron SED. This means that the effect of any potential negative curvature between 408 MHz and 30 GHz, as, for instance, assumed by Planck Collaboration X (2016), will instead by interpreted as a slightly steeper spectral index in the current analysis. This is more explicitly demonstrated in Sect. 6.
For further information regarding this model and a brief discussion of each individual component, we refer to BeyondPlanck Collaboration (2023) and references therein. The main goal of the present paper is to establish efficient sampling algorithms for the amplitudes and spectral parameters in Eq. (13)–(18).
3. Data selection
3.1. BEYONDPLANCK data selection
As discussed in Sect. 2 of BeyondPlanck Collaboration (2023), the only data set which is considered at the timeordered level is the Planck LFI data. With the minimal sky model, discussed in Sect. 2.2, a total of four unpolarized astrophysical sky components are considered. Seeing as LFI only contains three frequency channels, it is clear that the LFI data itself is unable to properly constrain this sky model. A set of selected external data is therefore included in the foreground analysis in order to constrain the sky model presented within BEYONDPLANCK.
Table 2 provides an overview of all frequency maps included in the intensity component separation procedure. We note again that the main motivation underlying the BEYONDPLANCK analysis is not to derive a novel stateoftheart intensity sky model, but rather to develop and demonstrate the Bayesian endtoend analysis framework using Planck LFI as a worked case. Accordingly, to ensure that the main results are dominated by Planck LFI, all CMBdominated Planck HFI bands and the WMAP Kband channel are excluded from the analysis; the data summarized in Table 2 represent a minimum set that is able to algebraically resolve all main foreground components relevant for Planck LFI. All sky maps (including LFI and others) are discretized using the HEALPix^{1} (Górski et al. 2005) pixelization.
Frequency band summary for the BEYONDPLANCK intensity analysis.
For all nonLFI bands, we adopt nominal bandpass profiles as recommended by the respective references. However, as an exception, we adopt the simplified (and commonly used) delta function approximation for the Haslam 408 MHz channel. For LFI 30 GHz, we allow for both an absolute bandpass shift for the full frequency band and relative differences between individual detectors; whereas for the 44 and 70 GHz channels, we only allow relative detector shifts, but with no overall absolute shifts; see Svalheim et al. (2023a) for further details.
The noise is assumed to be uncorrelated and Gaussian for all channels except Planck LFI, with a spatially varying root mean square (RMS) as defined by the number of hits per pixel. Again, the only exception is Haslam 408 MHz, which is nominally strongly signaldominated per pixel and dominated by systematic uncertainties, not statistical. In this case, we instead adopted a noise rms model that is the sum of an isotropic 0.8 K term (representing statistical uncertainties) and 1% of the actual map itself, representing multiplicative uncertainties; this is the same approach as taken by Planck Collaboration X (2016). For LFI, correlated noise is accounted for on all angular scales through explicit timedomain sampling, as discussed by Ihle et al. (2023).
In the data model given in Eq. (1), only the orbital CMB dipole and the far sidelobes are modeled with the full asymmetric beams. For astrophysical component modeling, all beams are assumed to be azimuthally symmetric, with window functions, b_{ℓ}, provided individually by each experiment. No uncertainties on these are propagated in the current analysis, but support for this will be added in future work. The Haslam 408 MHz and Planck DR4 (NPIPE) 857 GHz maps are smoothed from their native resolutions to 60′ and 10′, respectively. The latter is additionally repixelized to a HEALPix resolution of N_{side} = 1024 to reduce CPU and memory requirements; note that this channel still has a higher angular resolution than the 70 GHz LFI channel, which is the highest resolution channel of the main BEYONDPLANCK analysis.
An important and novel aspect of the current analysis is componentbased monopole (or “zerolevels” or “offsets”) determination, as discussed in Sect. 4.2. Rather than attempting to set the monopoles for each channel before component separation, we impose physical priors on the monopole for each astrophysical component. This astrophysical model is then used to determine deterministically the zerolevel for each frequency map. Reasonable priors may be defined for all components except synchrotron emission, and for this component we instead adopt explicit literature values. Specifically, we adopted a monopole value of 8.9 ± 1.3 K_{CMB} for synchrotron emission at 408 MHz, as estimated by Wehus et al. (2017) and we thereby neglect possible contributions from freefree emission to Haslam 408 MHz outside the very conservative Galactic mask employed by Wehus et al. (2017). We also applied the dipole corrections to the Haslam map derived by the same analysis.
Any additional preprocessing applied to the various maps was kept at a minimal level. Specifically, for WMAP we added the WMAP solar CMB dipole of (d, l, b) = (3355 μK, 263.99° , 48.26° ) to each map (Hinshaw et al. 2009); while for Planck 857 GHz, we apply a zodiacal light correction, following Planck Collaboration Int. LVII (2020).
No calibration corrections are applied to any nonLFI data sets, and we thus rely on the calibration of the original analyses for these channels. This is particularly important with respect to the WMAP channels, which have a nonnegligible impact on the solar CMB dipole; consequently, the final BEYONDPLANCK solar dipole estimate represents a noiseweighted average between the WMAP and BEYONDPLANCKbased LFI estimates.
3.2. Simulation data
Proper testing of the algorithms presented in the current paper constitutes a vital component in the verification of the feasibility of these algorithms within the Commander framework. As such, a suite of simulated data (with controlled noise) and total offsets was constructed. These simulated data are created to represent the full BEYONDPLANCK sky model, with frequency channels equivalent to the conditions placed on the data selection (discussed in Sect. 3.1 and presented in Table 2).
The simulations are entirely created within the mapspace domain. Using a sample from the BEYONDPLANCK sky model ensemble, mock frequency maps were created. Using Commander we output each of the foreground components, as determined by the BEYONDPLANCK sky model, at each of the input frequency channels listed in Table 2. As a result, we were able to create mock frequency maps by coadding each of these sky components at the corresponding frequencies. Thanks to this simple method, the full sky model, including the spectral parameters, are encapsulated in the set of simulated frequency maps.
Noise was then added to each of the simulated frequency maps by taking a random realization of the noise rms sky maps, assuming that the noise is white. For the noise realizations, we utilize the actual noise rms data as used and produced within the full BEYONDPLANCK results.
4. Commander extensions for efficient intensitybased component separation
The specific BEYONDPLANCK computer code implementation is called Commander3 (Galloway et al. 2023a), and this is a direct generalization of the code first introduced for CMB power spectrum estimation purposes by Eriksen et al. (2004) and later generalized to also account for astrophysical component separation by Eriksen et al. (2008), Seljebotn et al. (2014,2019). It was one of four main component separation algorithms adopted by the Planck collaboration (Planck Collaboration XII 2014; Planck Collaboration X 2016; Planck Collaboration IV 2020; Planck Collaboration Int. LVII 2020). Unless it is useful for context, we do not distinguish between the different code versions and we simply refer to all versions as Commander. In this section, we describe the four algorithmic improvements we have made to Commander, as introduced in Sect. 1.
4.1. Joint amplitude and spectral parameter sampling
As summarized in Eqs. (5)–(11), the BEYONDPLANCK pipeline implements a Gibbs sampling chain iterating over all free parameters in the data model. While Gibbs sampling in general is a very powerful method for exploring complicated distributions, its main weakness is the inability to probe degenerate distributions. This problem is illustrated for a toy example in the left panel of Fig. 1: the black contours represents the 68 and 95% confidence limits of a twodimensional Gaussian distribution with a Pearson’s correlation coefficient of ρ = 0.99. The purple point indicates the starting position of a Markov chain, while the dashed, grey horizontal line indicates the baseline of the corresponding conditional distribution P(β ∣ a_{init}), which itself is shown as an orange curve. Since the Gibbs sampling algorithm works by moving according to conditional distributions alone, the allowed step size in each iteration is very narrow compared to the full marginal distribution, shown as a blue distribution.
Fig. 1. Illustration of conditional (pure Gibbs; orange) and marginal sampling (blue) algorithms for a highly correlated (Pearson’s correlation coefficient of ρ = 0.99) twodimensional Gaussian distribution (black contours). The initial position (β_{init}, a_{init}) is indicated by a purple dot. Left: comparison of unnormalized conditional P(β ∣ a_{init}) distribution evaluated at the initial position and the corresponding marginal P(β) distribution; note that the latter is much wider than the former. Assuming β samples drawn at the 10th percentile, the graphs along the vertical lines represent unnormalized conditial distributions P(a ∣ β) evaluated at the β values drawn with the conditional (orange) and marginal (blue) distributions of β; note that the marginal sampling case results in a much longer step length between the initial and final sample values. Right: samples of a standard Gibbs sampling chain (orange) using conditional sampling for both a and β, and a sampling chain (blue) using marginal sampling for β and conditional sampling for a. Both cases show the first 100 samples initialized from the purple point. 
To illustrate the step size effect in the Gibbs sampler, we consider a Gibbs move in Fig. 1, starting with β from the purple point. The relevant conditional distribution for β is shown as an orange curve along the horizontal dashed gray line passing through the purple point and we can thus draw a random value from this distribution. This could for instance be the value indicated by the rightmost vertical dashed gray line. According to the Gibbs sampling algorithm, we are required to draw a sample from the corresponding conditional distribution for a, which is indicated by the orange distribution aligned with the vertical dashed line. One possible outcome of this, after completing one full Gibbs iteration, is the orange point. Now, because each conditional distribution is much narrower than the corresponding marginal distributions, the relative Gibbs step size is very short, and it takes a very long time to move from one side of the joint distribution to the opposite. The result is poor Markov chain mixing and a very long correlation length. As a realworld illustration of this, the orange points in the right panel of Fig. 1 show the 100 first steps of an actual Gibbs chain with this precise target distribution. We see that less than half of the distribution is actually explored, and many thousands of samples will be required in order to probe the full distribution with this algorithm.
This problem is directly relevant for modern Bayesian intensitybased CMB component separation. For experiments such as Planck and WMAP, characterized by very high S/N, there are strong degeneracies between the foreground amplitudes, the foreground spectral parameters, and maplevel monopoles. Explicitly, if one assumes that all spectral parameters and monopoles are known, then the conditional amplitude uncertainty is very small. Conversely, if we assume the amplitude and monopoles to be known, then the conditional spectral parameter uncertainties are small. However, when all parameters are unknown, the full uncertainties are significant.
For this reason, Gibbs sampling should usually be considered a last resort to handle an otherwise intractable distribution. If direct joint sampling methods are available, then those are usually more efficient. Fortunately, the Gibbs sampling method can be interleaved by any combination of conditional and joint steps while still maintaining the requirement of detailed balance (Geman & Geman 1984); also, the more steps that can be handled jointly, the more efficient the overall Gibbs chain will be. For the purposes of intensitybased CMB component separation, we therefore introduce a new specialpurpose joint amplitude–spectral parameter step by exploiting the definition of a conditional distribution as follows,
The first distribution on the right hand side is the marginal distribution of β with respect to the data m_{ν}, and the second distribution is the conditional distribution of a with respect to β. This equation therefore implies that we may generate a joint sample by first drawing β from its “marginal” distribution, and then sample a from the corresponding “conditional” distribution:
We note the absence of a in the first distribution. We then need to derive sampling procedures for each of these two distributions. According to Bayes’ theorem:
where P(m_{ν}) is just a normalization factor (often called the “evidence”), P(a, β) denotes optional priors, and the final factor is the likelihood function, P(m_{ν} ∣ a, β)≡ℒ(a, β). From the compact data model in Eq. (3), we note that
and since is assumed to be zeromean and Gaussian distributed with known variance, we can immediately use the following expression for the likelihood function:
where N_{ν} is the noise covariance matrix of band ν, which is diagonal in the case of pure white noise. The priors are less welldefined, and are left to the user to determine. In the following, we adopt Gaussian priors for spectral parameters and for notational convenience, we assume no spatial amplitude priors, S^{−1} = 0.
We first consider the marginal spectral parameter distribution, P(β ∣ m_{ν}). This is derived by integrating Eq. (22) with respect to a, and this was done by Stompor et al. (2009) and Stivoli et al. (2010) as part of developing the MIRAMARE component separation code. The result is expressed as:
where all terms should be interpreted as sums over frequencies. The same authors also introduced a socalled “ridge likelihood,” in which one does not marginalize over a, but rather sets a equal to its maximum likelihood value for a given value of β. This may also be analytically evaluated, and is in fact identical to the above expression, with the exception of the last determinant term being excluded. We implemented support for both options in our codes. We note that this expression requires all data to be defined at the same angular resolution and so, all data must be smoothed to a common resolution before evaluating Eq. (23).
The second required distribution is P(a ∣ m_{ν}, β). This is a simple multivariate Gaussian distribution in a, for which there are efficient samplers readily available (see, e.g., Appendix A of BeyondPlanck Collaboration 2023 for details). One particularly efficient sampling equation is as follows:
where η is a vector of random Gaussian N(0, 1) variates. This equation may be solved efficiently using preconditioned Conjugate Gradient methods (Shewchuk 1994), as discussed by Seljebotn et al. (2019).
These equations are integrated into the main BEYONDPLANCK Gibbs sampling loop according to the following steps. First, we run a short (typically a few hundred steps) standard Metropolis sampler (see Appendix A of BeyondPlanck Collaboration 2023) for each spectral parameter, using the product of Eq. (23) and any desired priors to define the accept rate, that is, the relative number of Metropolis proposals being accepted (which should preferably stay between 0.3 and 0.7 for an efficient sampler). All data are smoothed to a common angular and pixel resolution before evaluating the expression. Immediately following the last Metropolis step, we draw one sample from P(a ∣ m_{ν}, β) using Eq. (24); it is critically important that no other parameters are updated between β and a, as the previous value of a is completely inconsistent with the new β value, which is drawn marginally with respect to a.
Returning to Fig. 1, the improvement achieved by this joint twostep sampler is illustrated as blue distributions and points. Starting with the left panel, the fundamental difference between the joint and Gibbs samplers is that the first step in the βdirection is drawn from the full marginal distribution (horizontal blue distribution) instead of the conditional distribution (horizontal orange distribution). This is much wider, and covers by construction the full width of the underlying target distribution. One single proposal may therefore move from one side of the distribution to the other, and there is no memory of the previous parameter state. However, to obtain a valid sample, it is critically important to draw a corresponding sample from the appropriate conditional amplitude distribution (vertical blue distribution) immediately after the marginal move. Correspondingly, the blue points in the right panel shows 100 samples drawn with the joint sampler. In this case, they cover the full distribution very efficiently.
Figure 2 shows a similar comparison for AME ν_{p} for a test chain that only explores the AME parameters in the BEYONDPLANCK data model with both methods. Also in this case, we see that the marginal sampler explores the full range much more efficiently than the conditional sampler.
Fig. 2. Comparison of AME ν_{p} chains derived using the conditional (orange; Eq. (22)) and marginal (blue; Eq. (23)) samplers discussed in the text. For the purposes of this illustration, all other parameters than a_{AME} and ν_{p} are fixed. 
4.2. Componentbased monopole determination
Next, we considered the problem of monopole determination for CMB experiments, which has long been one of the main challenges for parametric component separation methods (see, e.g., Planck Collaboration X 2016; Wehus et al. 2017, and references therein). The problem stems from the following challenge: For traditional CMB experiments and maximum likelihood mapmaking methods, there are no datadriven constraints on the monopoles in the derived frequency sky maps. For example, WMAP is explicitly differential in nature, measuring only differences between pairs of points, and therefore cannot by construction constrain the zerolevel. For Planck, the high level of 1/f noise prohibits any useful constraints on the zerolevels. An important exception to this is COBE/FIRAS (Mather et al. 1994), which is absolutely calibrated; but its mKlevel uncertainties are still orders of magnitude too large to be useful for modern CMB component separation purposes.
For this reason, several indirect methods have been established to determine the frequency map monopoles based on the morphology of the maps themselves. Four examples are mean subtraction in a small region (Planck Collaboration V 2014); fitting a planeparallel cosecant model (Bennett et al. 2003; Planck Collaboration II 2016); imposing foreground SED consistency between neighboring frequencies (Wehus et al. 2017); and crosscorrelation with external data sets with known zerolevels (Planck Collaboration VIII 2016). However, all of these methods have in common the fact that they operate on the basis of frequency maps and are aimed to determine the zerolevel at a given frequency channel, before feeding these into traditional component separation algorithms. In this study, we have made the observation that it is, in fact, much simpler to determine the monopoles of the component amplitude maps and to then use these to deterministically set the frequency map monopoles through the resulting sky model. The frequency map zerolevels have thus no independent impact on any higher order analyses (most notably, the spectral parameters), but simply adjust to whatever the model dictates at any given moment.
As a result, the question that immediately rises considers how we may, in fact, determine the component monopoles. This must be done on a casebycase basis, applying the most natural prior for each component. We note that any true monopole signal in the components that do not agree with the chosen priors will end up in the frequency map monopoles. Starting with the CMB case, this can either be set to zero or 2.7255 ± 0.0006 K (Fixsen 2009), depending on whether we want sky maps without or with the CMB monopole. In practice, we additionally account for suboptimal foreground modeling by applying a mask. For the current analysis, we derived the CMB monopole mask from a set of smoothed component amplitude maps, namely, by thresholding the sum of synchrotron, AME, freefree and thermal dust emission, all smoothed to 10° FWHM. In addition, we masked out radio sources and any pixel with a reduced normalized χ^{2} higher than 5 σ. The resulting mask is shown in Fig. 3, and has an accepted sky fraction of f_{sky} = 0.64. The monopole of the CMB component map is set to zero (or 2.7255 K) outside this mask, while simultaneously fitting for (but not modifying) the dipole component.
Fig. 3. Masks used for signal amplitude zerolevel and band monopole sampling: Frequencyband monopole and CMB amplitude zerolevel (top) and freefree amplitude zerolevel (bottom), with a freefree amplitude sample at an angular resolution of 30′ FWHM plotted underneath. 
Regarding the freefree component, Planck Collaboration X (2016) found that the measured emission is strongly noisedominated over large areas of the sky, with no detectable amplitude. Also, in this case, we therefore set the monopole to have zero mean outside a conservative mask. In this case, the mask is derived from the freefree amplitude map itself, evaluated at 30 GHz and smoothed to 10° FWHM, and truncated at 5 μK_{RJ}. In addition, we excluded areas where the other foreground signals are high to account for signalleakage into the freefree component, similarly to the CMB mask, but while thresholding the sum of all other components evaluated at 44 GHz. The resulting mask is shown in the bottom panel of Fig. 3, accepting f_{sky} = 0.50. We do note that any true unmasked freefree signal is by definition positive and this can bias the monopole also in the noisedominated regime. Future works should aim to correct for this bias by directly estimating the residual freefree monopole in the unmasked region. We do note, however, that this bias will decrease as more high sensitivity data become available, as more and more of the freefree signal may be masked directly.
In contrast, synchrotron emission as observed by the Haslam 408 MHz map is highly diffuse on the sky and there are no regions on the sky that can be assumed to be approximately clean of synchrotron emission, namely, exhibiting no synchrotron signal. In this case, the best estimates of the synchrotron amplitude are those already estimated for the Haslam 408 MHz map itself. In this paper, we adopt the zerolevel correction of 8.9 ± 1.3 K_{CMB} derived by Wehus et al. (2017). We marginalized over the uncertainty by drawing a random offset correction in every Gibbs iteration, as defined by a Gaussian distribution with the quoted mean and standard deviation.
For the thermal dust emission, we adopted essentially the same approach as the Planck HFI DPC (e.g., see Planck Collaboration III 2020), setting the zerolevel through crosscorrelation with H I column density observations (e.g., see Lenz et al. 2017), although with a few minor variations. First, we applied this method to the thermal dust component map, as opposed to individual frequency maps. Second, for Planck DR4 (Planck Collaboration Int. LVII 2020), the HFI 545 GHz zerolevel was set through a linear fit for pixels with N_{H I} < 4 × 10^{20} cm^{−2}. However, as the 545–N_{H I} scatter plot appeared to be nonlinear around values of N_{H I} = 1.5 × 10^{20} cm^{−2}, they also performed a second degree fit. We show in Fig. 4 a similar scatter plot between N_{H I} and the thermal dust amplitude from one of our Gibbs samples, where we performed both a first and second degree polynomial fit to the plot at N_{H I} < 4 × 10^{20} cm^{−2}. Furthermore, we also performed a linear fit at N_{H I} < 2 × 10^{20} cm^{−2} and we see that the intersection of the linear fit with a lower threshold is close to the intersection of the second degree fit. This raises the question of the uncertainty of the threshold value for the linear fit, as Planck Collaboration XXIV (2011) found a good correlation up to at least N_{H I} = 2 × 10^{20} cm^{−2}. We therefore implement this crosscorrelation method as a prior on the thermal dust zerolevel using a range of thresholds, and for each Gibbs iteration we perform linear fits of N_{H I} and the thermal dust amplitude with N_{H I} threshold values ranging from 1.5 to 4 [10^{20} cm^{−2}] with increments of 0.5. Then we draw the intersection value from a Gaussian distribution given the mean and variance of the linear fits, which is subtracted from the dust amplitude. This way, the uncertainty of the thermal dust amplitude zerolevel also propagates through the pipeline.
Fig. 4. Crosscorrelations between the H I column density N_{H I} and (top) the thermal dust amplitude and (bottom) the AME amplitude. The thermal dust crosscorrelation is evaluated at HEALPix resolution N_{side} = 64 and a common angular resolution of 60′ FWHM, while the AME crosscorrelation is evaluated at HEALPix resolution N_{side} = 16 and a common angular resolution of 10° FWHM. The lines represent bestfit lines for pixels with H I column densities less than 2 × 10^{20} cm^{−2} (blue) or 4 × 10^{20} cm^{−2} (orange). The green curve is the bestfit second degree polynomial to pixels with H I column densities less than 4 × 10^{20} cm^{−2}. The contour lines are plotted at 0.001, 0.01, 0.05, and 0.1 ; where only the lower three contour line values are plotted for AME. The contours have been smoothed for visualization. 
The zerolevel of the AME component is determined using the same procedure, noting that Planck Collaboration X (2016) demonstrated a very tight spatial correlation between AME and thermal dust emission on large angular scales. The only difference with respect to the thermal dust procedure is that we smooth all maps to a common angular resolution of 10° FWHM and a HEALPix resolution of N_{side} = 16, and adopt thresholds of 2 to 4 [10^{20} cm^{−2}], with increments of 0.5; both the smoothing and the lower resolution are imposed to reduce the impact of instrumental noise. A scatter plot between AME and N_{H I} is shown in Fig. 4 for one arbitrary Gibbs sample.
With the introduction of componentbased monopole priors, all frequencyband monopoles become free parameters and can be deterministically fitted. Explicitly, for each frequency channel we first subtract the predicted sky model as defined by Eq. (2) and then fit and subtract the residual monopole outside some mask. The mask should be defined such that it excludes areas on the sky prone to foreground mismodeling, hence we adopt the same mask as we do for the CMB monopole prior, shown in the top panel of Fig. 3. We note that this monopole adjustment needs to be done immediately after any change in any of the component maps, in order to not break the Gibbs chain, fully analogously to the immediate amplitude update that must follow any marginal spectral parameter move discussed in the previous section.
4.3. Joint spectral parameter and frequencyband monopole sampling
Returning to the AME–H I crosscorrelation plot in Fig. 4, we notice that the zerolevel is associated with a large statistical uncertainty. When sampling the AME peak frequency, ν_{p}, this uncertain monopole is also directly affected by the resulting SED changes, and corresponding monopole offsets are induced at all frequencies. If ν_{p} is sampled conditionally with respect to the band monopoles, these will therefore tend to pull ν_{p} towards the old value and thereby increasing the overall Markov chain correlation length.
This inefficiency may be alleviated by exploiting the new componentbased monopole sampler described in Sect. 4.2. Since all frequencyband monopoles are now deterministically defined by the sky model, these can be adjusted jointly whenever that is modified. We therefore implement internal estimation of frequency band monopoles during the spectral parameter sampling algorithm, such that for each proposed spectral parameter value, we estimate a new band monopole value conditioned on the input amplitude and the proposed parameter value. The frequency band monopoles are updated together with the final parameter at the end of the sampling.
To illustrate the usefulness of this combined sampling step, we generate an idealized simulation that includes only AME signal and noise at each frequency channel, with an angular resolution of 10° FWHM and a HEALPix resolution of N_{side} = 16. We choose an input peak frequency of ν_{p} = 26 GHz, and adopt a zerolevel from H I column density crosscorrelation to 2 μK_{RJ}. Figure 5 shows the resulting loglikelihood (or χ^{2}) distributions as evaluated from the marginal definition in Eq. (23), for cases both with (green curve) and without (orange curve) marginalizing over the frequency band monopoles. We see that by marginalizing over the band monopoles the loglikelihood function widens by a factor of 1.5–2. This translates into correspondingly longer Metropolis step sizes in the spectral index sampling steps, and thereby faster exploration of the full posterior distribution. The higher conditional S/N a given component has, the more important this effect will be.
Fig. 5. Comparison of AME ν_{p}χ^{2} distributions with (green) and without (orange) monopole marginalization. These distributions are evaluated using the marginal spectral parameter likelihood given in Eq. (23), but the same qualitative behaviour holds irrespective of which spectral parameter distribution is used (conditional, ridge, or marginal): the distribution becomes significantly wider when marginalizing over band monopoles. 
4.4. Breaking smallscale degeneracies through spatial priors
The final algorithmic improvement presented in this paper is the introduction of informative spatial priors for foreground components, either in the form of purely algorithmic smoothing power spectrum priors or as actual informative Gaussian priors with a nonzero mean. The first of these has already been used in the latest Planck analyses (Planck Collaboration IV 2020; Planck Collaboration Int. LVII 2020), but in the following, we generalize the approach to nonzero cases and we systematically show how different choices affect the final results.
Mathematically speaking, the only difference between an informative prior and a smoothing prior is whether a preexisting mean map is assumed for the astrophysical component in question (in which case the prior is called “informative”) or whether the prior mean is assumed to be zero. Practically speaking, however, there is also an important difference between the prior variances in the two cases, since for informative priors the variance quantifies the allowed level of fluctuations around the mean map; while for smoothing priors, it quantifies the allowed level of fluctuations around zero. Thus, for informative priors a prior variance of zero is fully acceptable, in which case the output component map will be identical to the prior mean; while for a smoothing prior the variance should be larger than the actual component fluctuations in order to avoid oversmoothing.
We start by revisiting the sampling equation for the component amplitude maps, as defined by Eq. (24). This equation provides a sample from a posterior defined only by the likelihood itself. If we additionally want to impose a Gaussian prior on the amplitudes, as defined by a multivariate Gaussian distribution, N(μ, S), then this is generalized to:
We refer to Appendix A in BeyondPlanck Collaboration (2023) for an explicit derivation. In this expression, μ has the same dimension as a, and represents the prior mean for a, while S is an associated prior covariance matrix that defines the “strength” of the prior, fully analogous to the usual standard deviation, σ, of a Gaussian univariate prior. Thus, if S = 0, the final solution for a will be identical to μ, while if S → ∞ (or, equivalently, S^{−1} = 0), then the prior term vanishes, and one is left with the original likelihood expression in Eq. (24). For reference, we note that the previous Planck analyses (Planck Collaboration IV 2020; Planck Collaboration Int. LVII 2020) set μ = 0 in this equation and only used S to impose smoothness on a.
Computationally speaking, introducing informative priors with nonzero means in Eq. (25) represents no additional algorithmic complications compared to the priorfree case: the equation is in both cases solved using the same preconditioned conjugate gradient implementation. If anything, the equations are actually a bit easier to solve with informative priors, as they reduce degeneracies between different parameters, and thereby reduce the condition number of the coefficient matrix on the lefthand side of the equation. As pointed out by Seljebotn et al. (2019), from an algorithmic point of view an informative prior defined in terms of a mean map with a specified covariance may simply be considered to be a new independent data set with sensitivity only to the component in question and it therefore provides orthogonal information with respect to the likelihood term contributions. Rather, the main challenge regarding priors is how to define them in a useful and controlled manner that does not significantly bias or contaminate the final posterior distribution and this must be assessed on a casebycase basis.
In the current BEYONDPLANCK analysis, we followed Planck Collaboration IV (2020) and define S in harmonic space, giving different prior weights to different angular scales. Explicitly, each component map is defined in terms spherical harmonics:
and we define the prior covariance matrix as:
where P_{ℓ} is an angular prior power spectrum for a, which, again, is fully analogous to the standard deviation of a Gaussian prior, but now defined per angular multipole.
4.4.1. Algorithmic smoothing priors
Starting with the algorithmic smoothing priors adopted by Planck 2018, we note that if we define to be an estimate of the true angular power spectrum of a, then the following power spectrum prior,
simply represents a Gaussian smoothing prior of a with a smoothing kernel width equal to σ, which often is defined in terms of . We call this an algorithmic smoothing prior, as it explicitly pushes the solution to be smooth on small angular scales. We also note that does not have to be an accurate estimate of the true component power spectrum, but it should in general be greater than the true spectrum, in order to prevent the prior from being overly constraining. This is again fully analogous to choosing a prior width that is wider than the expected target distribution in standard univariate analysis problems. Thus, simply setting to a constant that is a few times larger than the expected spectrum is usually a perfectly good choice.
It is important to note that Eq. (28) is indeed just a prior, and not a deterministic postprocessing smoothing operator. This has both advantages and drawbacks that are important to be aware of when using the products from the analysis. The main advantage is that signaldominated localized objects (for instance point sources) are not excessively smoothed when applying the smoothing as a prior. The main disadvantage is that the effective angular resolution of the amplitude map becomes spatially varying and depends on the local S/N in a given pixel; if the S/N is high, the angular resolution will be determined by the resolution of the data, while if the S/N is low, it is given by θ_{FWHM}. This is similar to the GNILC method (Remazeilles et al. 2011), which also implements S/N dependent angular resolution. The main difference between the two methods is that while GNILC requires regions of different resolutions to be predefined, the current approach automatically and dynamically adopts the resolution while solving Eq. (25).
In the current BEYONDPLANCK analysis, where the main scientific target is CMB power spectrum and cosmological parameter estimation, we do not impose any priors on the CMB component amplitudes, but we do impose a Gaussian smoothing prior for thermal dust emission in intensity with θ = 5′ FWHM and at 545 GHz. For studies that are primarily interested in the angular power spectrum of thermal dust emission, it would be more useful to instead impose a LambdaColdDarkMatter (ΛCDM) spectrum on the CMB amplitude map, and no priors on the thermal dust spectrum. Then the resulting thermal dust power spectrum would be an unbiased estimator; with the current analysis, the thermal dust spectrum will be biased low on small angular scales due to the smoothing prior. We also impose a Gaussian smoothing prior on synchrotron emission in intensity, with θ = 60′ FWHM and at 408 MHz.
4.4.2. Informative Gaussian spatial priors
For AME and freefree emission, we adopt informative priors with μ ≠ 0. The reason for this is simply that the limited data combination considered in the current BEYONDPLANCK analysis (see Table 2) is inadequate for constraining all of AME, freefree, and CMB separately without additional information; when future observations from, for instance, CBASS (King et al. 2010) and QUIJOTE (GénovaSantos et al. 2015) become publicly available and integrated in the analysis, this will hopefully no longer be necessary. As an illustration of the problem, Fig. 6 compares the AME and freefree amplitude maps derived without any priors near the Galactic South Pole; even at a visual level, these two maps are nearly perfectly anticorrelated, with no true constraining power on their own. This also makes other components (most importantly the CMB) susceptible to small systematic residual mismatches between the AME and freefree components and it significantly increases the CMB noise.
Fig. 6. Comparison of AME (left) and freefree (right) amplitude maps derived without any spatial priors in a 20 × 20° field centered on the Galactic South Pole, (l, b) = (0° , − 90° ) at an angular resolution of 60′ FWHM. We note the striking anticorrelation between the two maps. 
Starting with the AME case, we first note that the priorfree Planck 2015 analysis (Planck Collaboration X 2016) found a very strong spatial correlation between their AME and thermal dust component maps at an angular resolution of 1° FWHM. In general, a high degree of correlation between these components is expected from current theoretical AME models (e.g., Erickson 1957; Draine & Lazarian 1998; AliHaïmoud 2010; Silsbee et al. 2011; Hensley & Draine 2020), although the specific correlation coefficient depends on model details. Based on these observations, we adopted the Planck DR4 857 GHz map^{2} (Planck Collaboration Int. LVII 2020) as a spatial mean template for AME, that is, μ in Eq. (25). However, before it can be inserted into Eq. (28), it must be adjusted in amplitude to account for the mean SED difference between thermal dust emission at 857 GHz and AME at 22 GHz. To do this, we solved Eq. (28) using the 857 GHz map scaled by a range of values between α = 2 × 10^{−5} and 4 × 10^{−5} as the AME prior. We then took the difference between the derived amplitude map and the input prior and adopted the scaling factor for which the difference is smallest as our default prior. Example difference maps are shown in Fig. 7 and we adopted α = 3 × 10^{−5} as our default prior.
Fig. 7. Difference maps between the derived amplitude and prior maps for AME for three different 857 GHz scaling factors. From top to bottom, the three panels: scaling factors of α = 2 × 10^{−5}, α = 3 × 10^{−5}, and α = 4 × 10^{−5}. 
The final step is to define the strength of this prior, as given by P_{ℓ}. Ideally, we want the prior to be stronger (i.e., S to be smaller) for the noisedominated small angular scales and looser for the signaldominated large angular scales. To quantify these considerations, we must define the following a powerlaw prior power spectrum for AME:
where q is an overall amplitude at a pivot multipole, ℓ_{0}, and β is a tilt parameter. A negative (positive) β results in a stronger (weaker) prior at high multipoles and a smaller (higher) amplitude, q, gives a stronger (weaker) prior on all angular scales.
Figure 8 compares the resulting AME amplitude maps for various choices of q and β (bottom panels) with the prior mean map (top panel). Here we see that a high value of leads to very similar solutions for β = −1 and β = 1, indicating that the prior is largely irrelevant, and the solution is data dominated. We also note that there is substantial instrumental noise at high latitudes. For , the maps are notably smoother at high latitudes, but we also see clear smoothing artefacts near the Galactic plane, corresponding to harmonic space ringing from the Galactic plane. A value of β = −1 leads to sharper edges than β = 1.
Fig. 8. Effects of the spatial prior on the sampled AME amplitude. AME amplitude prior map, derived by scaling the Planck DR4 857 GHz by α = 3 × 10^{−5} and smoothing to 10′ FWHM (top). Derived AME amplitude maps for four different spatial prior combinations (bottom), . Rows: results for and , respectively, while columns show results for β = −1 and 1. 
To actually determine the parameters used for the final prior, we solve Eq. (28) over a grid in q and β and evaluate, as follows,
for each configuration, where s is the derived sky model in each case. The results from this evaluation are shown in the top panel of Fig. 9 in terms of the normalized reduced where n_{d.o.f.} is the number of degrees of freedom; for a perfect model fit and n_{d.o.f.} ≫ 1, this quantity should be distributed approximately as a Gaussian distribution with zero mean and unit standard deviation. For and β = 0, we see that , which indicates a clear excess residual. However, for and β = 0, this excess is greatly diminished, while there still is some effect from the prior. In the following, we adopt this latter combination as our default AME prior.
Fig. 9. Average per N_{side} = 16 pixel as a function of AME (top) and freefree (bottom panel) prior amplitude, q, where and n_{dof} = 15 400. Colored solid lines show results for different tilt parameters, β, and the intersection between the thick green curve and dashed line indicates the prior combination adopted for the main BEYONDPLANCK analysis. 
For freefree emission, there are no corresponding fullsky independent spatial templates available in the literature. Observations of Hα (Finkbeiner 2003) or radio recombination lines (RRL; Alves et al. 2015) might serve useful roles, but both are associated with significant shortcomings for the purposes of the current analysis: the Hα observations lack most of the Galactic plane signal due to dust absorption, while the RRL observations only cover a part of the Galactic plane. For now, the best available fullsky freefree tracer is in fact the Planck 2015 freefree map (Planck Collaboration X 2016), which is based on the same data set as studied in the current paper, but additionally (and critically) the Planck HFI observations as well. We therefore adopted this map as a spatial prior in the current analysis, while recognizing that this is strictly speaking not admissible in the Bayesian framework; some data (i.e., LFI, WMAP, and Haslam) are used twice to constrain freefree emission and the resulting uncertainties will therefore be underestimated. In practice, this solution is a way of integrating HFI observations into the analysis without directly affecting the CMB component. A critical goal for nearfuture work is to integrate HFI observations directly into the analysis in the form of frequency maps and at that point, this informative freefree prior will be removed.
We adopted the same parametric function for freefree emission as for AME, defined by Eq. (29) and adjusted the free parameters in the same way. The results from this optimization are shown in the bottom panel of Fig. 9 and we adopted and β = 2 in this case. A comparison of different prior choices with the actual input prior map is shown in Fig. 10.
Fig. 10. Freefree amplitude prior map, adopted from the Planck 2015 analysis which includes HFI observations (top). Derived freefree amplitude map for four different spatial prior combinations (bottom), . Rows show results for and , respectively, while columns show results for β = 1 and 3. 
The only algorithmic difference with respect to AME is that we additionally imposed a Gaussian smoothing for freefree emission, as per Eq. (28) with θ = 30′ FWHM. This is done to account for the fact that the distribution of freefree emission is highly localized on the sky and therefore requires a high maximum multipole moment to capture all significant structures; the additional Gaussian smoothing ensures that no ringing emerges from the highℓ truncation.
5. Validation by simulations
In order to validate the component separation implementations, we ran Commander on the simulated data as described in Sect. 3.2. In this section, we describe the efficacy of the algorithmic developments by comparing the input foreground sky maps with the resulting mean of an ensemble of 300 component separation samples as produced by Commander. As described in Sect. 6.1, the only spectral parameter which is sampled and constrained by data is the AME peak frequency, ν_{p}, which is determined as a full sky value.
For the amplitude maps, we defined a new variable to represent the relative difference between the input amplitude map and the results of running the simulation through the Commander component separation procedure. We define ϵ as the relative difference given by:
where, again, c is an astrophysical sky component, a^{in} is the input simulation amplitude map, and a^{out} is the mean amplitude map of the component separation samples.
The results of the amplitude maps, represented in terms of the relative difference maps defined in Eq. (31), are summarized in Fig. 11. Inspecting the difference maps shows that there is good agreement with the input amplitude map, with departures from small relative differences in parts of the sky which are easy to explain. Notably, both the freefree and AME components show high levels of noisy departures off of the Galactic plane. Compared to the relative difference between two samples within component separation, the relative differences seen here are small. We see that in the high S/N observations along the Galactic plane, the agreement is excellent, with the outline of each components amplitudes clearly defined.
Fig. 11. Relative difference maps, ϵ_{i}, for each of the sky components within the simulated BEYONDPLANCK dataset. Top row: AME, CMB, and thermal dust maps. Bottom row: freefree and synchrotron maps. 
For the other three sky components, we see even better agreement. Unsurprisingly, the thermal dust amplitude map shows excellent agreement over the majority of the sky, though with deviations at the 10–20% level in the low signaltonoise regions of the sky. Finally, the synchrotron amplitude map shows very small differences in ϵ_{s}, though the imprint of the component is significantly less notable here than in the other components.
The results of the 300 samples for the full sky spectral parameters and band monopoles can be seen in Figs. 12 and 13 respectively. Much like the sky component maps, we see excellent agreement with the input values for both spectral parameters and the band monopoles. We note that both the spectral indices and the band monopoles show a few samples which are significant deviations from the input value. This is to be expected, and is in fact an intention of the algorithms implemented in this work. As described in Sect. 4.3, the marginalization over the band monopoles allows for a broader loglikelihood, corresponding to a more complete exploration of the underlying distribution.
Fig. 12. Full sky spectral parameters as a function of sample (blue) for the controlled simulation. The nominal input value of the simulation is overlayed as the black dashed line. 
Fig. 13. Band monopoles for each of the simulated frequency bands as a function of sample. The nominal input value is given by the overlaid black dashed line. Right panel: the distribution of the samples. 
6. BEYONDPLANCK analysis and posterior distributions
We now turn our attention to the actual BEYONDPLANCK analysis and intensity component posterior distributions derived from the data combination discussed in Sect. 3. The results described in this section represent the intensity foreground results of the BEYONDPLANCK project and are a practical demonstration of the algorithms described and tested in Sects. 4 and 5 respectively. We once again note that the goal of the BEYONDPLANCK project is not to derive a new stateoftheart astrophysical component model within the relevant CMB frequency range (given that critical data sets such as Planck HFI and WMAP Kband are not included), but rather to lay the algorithmic groundwork for a statistically robust communitywide sky model, as will be implemented through the Open Science COSMOGLOBE^{3} community effort. As far as Bayesian intensity sky modeling is concerned, Planck Collaboration X (2016) still represents a state of the art approach.
6.1. Spectral parameter prior tuning
Before presenting the BEYONDPLANCK Markov chains and posterior distributions, there is still one task that must be completed before the algorithm described in Sect. 2 is carried out to completion, namely, finalizing the informative spectral parameter priors. With the reduced number of data sets included in this work, we have reduced constraining power when sampling spectral parameters, and strong priors are required for most spectral SED parameters. With this in mind, we can assume that all free parameters can only be fitted with a single constant value over the full sky, at least for now. Already at this stage, we fix the thermal dust temperature, T_{d}, at the sky map derived by Planck DR4 (Planck Collaboration Int. LVII 2020), noting that LFI has no constraining power for this particular parameter.
Even though LFI should have some constraining power of the thermal dust spectral index, β_{d}, the thermal dust and AME components are found to be highly degenerate with the limited data set used in this work. A joint fit of the AME ν_{p} and β_{d} would therefore lead to unphysical results, with preliminary analyses showing that the β_{d} diverged to values β_{d} > 2, raising ν_{p} to much higher frequencies. The uncertainty in the β_{d} value is important for error propagation and, thus, we must simply marginalize over the adopted prior, instead of trying to constrain β_{d} with the current data set.
For each free spectral parameter, we created a dedicated sampling mask, where we exclude regions on the sky where the other components are strong in order to reduce potential modeling mismatch errors to propagate between the various components. These masks are created from the amplitude maps of the modeled components, evaluated at 44 GHz and smoothed to 10° FWHM. In addition, we masked radio sources by thresholding the Planck 30 GHz compact source map at three different angular resolutions, namely, at native resolution and at 1° and 10° FWHM. For β_{s} and ν_{p}, we excluded regions of the sky where any other component signals is greater than 40 μK_{RJ}; while for T_{e}, we excluded the areas where the other components are greater than 50 μK_{RJ}. For all parameters, we exclude any pixel for which the smoothed radio sources are stronger than 30 μK_{RJ}. Finally, we masked out regions of the sky contributing to the largest 15% of a 10° FWHM χ^{2} map to exclude regions with large known modeling errors. The accepted sky fractions of the resulting masks are f_{synch} = 0.66, f_{ff} = 0.74, and f_{AME} = 0.66, respectively.
For each free parameter, we adopted a Gaussian informative prior with some mean and standard deviation, . The prior parameters are informed by literature results and listed in Table 1. For synchrotron emission, we note that few intensitybased constraints are available for frequencies higher than 30 GHz in the literature and we therefore adopted β_{s} = −3.3 ± 0.1, as derived from polarization measurements in Planck Collaboration V (2020). We also note that we have attempted to use flatter mean values of β_{s} = −3.1 and −3.0, as suggested from lowfrequency surveys, but these result in obvious artefacts in the current analysis in the form of a significantly overestimated synchrotron amplitude at 30 GHz. As an example, Fig. 14 compares the AME and freefree amplitude maps derived for two different values of β_{s}, namely, β_{s} = −2.8 and −3.3. While neither of these solutions produce an excess χ^{2} and, therefore, a free likelihooddriven fit is unable to distinguish between them, it is obvious from visual inspection that the former spectral index leads to clear synchrotron leakage into both the AME and freefree components. With a prior of β_{s} = −3.3 ± 0.1, the nonphysical flatindex solutions are largely excluded, while some parameter space is still allowed toward the steeper end to explore degeneracies. Another approach of reducing the predicted amplitude at CMB frequencies is by introducing a negative curvature in the synchrotron SED (as discussed in Sect. 2.2) and when comparing the results from this paper with previous results, it is important to ensure that the models are compatible.
Fig. 14. AME (left) and freefree (right) amplitude maps derived for two different synchrotron spectral indices. Top row: results for β_{s} = −2.8, while the bottom row shows results for β_{s} = −3.3. In the top row, we see a negative synchrotronlike (see Fig. 19) imprint, which is not present in the bottom row. 
Given the listed priors, we performed a coarse χ^{2} grid evaluation for each free parameter, conditioning on all other spectral parameter means, but allowing for amplitudes and band monopoles to adjust to the given spectral parameter. The resulting χ^{2}’s (evaluated at a HEALPix resolution of N_{side} = 16) are shown in Fig. 15.
Fig. 15. χ^{2} distributions (blue curves) from coarse grid evaluation of each free spectral parameters. From top to bottom, the panels show (1) AME peak frequency, ν_{p}; (2) synchrotron spectral index, β_{s}; and (3) freefree electron temperature, T_{e}. The minimum χ^{2} value for each parameter has been subtracted in each case. Orange curves show the priors adopted for the given component; see Table 1. The prior values for β_{s} and ν_{p} have been scaled by a factor of 100 to fit in the plot with the derived χ^{2} values. 
Starting with AME ν_{p}, shown in the top panel, we see that the χ^{2} is welldefined with a typical bestfit value around 22 GHz. The actual χ^{2} values show rapid increases at both lower and higher values with variations ranging in the thousands, and correspondingly, the prior (which is on the order of unity) is therefore largely irrelevant. It is clear that the current data combination has significant constraining power for ν_{p}.
The second panel shows similar results for β_{s}. In this case, we see a rapid χ^{2} increase for β ≲ −2.7, but it is otherwise slowly increasing for smaller values of β_{s} in the region β_{s} < −3.2, becoming almost flat. Additionally, we already know from Fig. 14 that spectral indices flatter than β_{s} ≲ −2.8 lead to clearly contaminated AME component maps, even if the χ^{2} is not able to identify this. At the same time, the actual χ^{2} variations are indeed larger than the prior, and this typically indicates that the algorithm prefers to use this unconstrained degree of freedom to fit other degrees of freedom, for instance, modeling errors in the thermal dust model. To avoid pathological solutions, we instead chose to marginalize explicitly over the prior and disable the likelihood term entirely when sampling this parameter. In other words, we simply marginalized over the adopted prior, but we did not attempt to constrain β_{s} with the current data set.
The same considerations hold to an even greater extent for the last parameter. For the electron temperature, T_{e}, the χ^{2} variations are entirely spurious and we therefore disable the likelihood term for T_{e}.
In summary, the only spectral parameter that the current data set is able to robustly constrain in intensity is the AME peak frequency, ν_{p}. All others are either drawn from their corresponding priors in the current analysis or frozen. We note that introducing additional data sets to constrain these parameters is a critically important next step for future works.
With β_{s} and β_{d} drawn from priors and T_{e} frozen, we find that the mask used in the coarse grid sampling of ν_{p} is too conservative when sampling ν_{p}, masking out too much of the galactic plane and leading to large dustlike residual signals in the LFI 30 GHz and the WMAP Ka channels. Additionally, a more dustlike signal was found to be leaking into the freefree component amplitude. In order to limit these effects, a less conservative mask had to be used. The mask used to sample ν_{p} in the final BEYONDPLANCK production is generated by excluding all pixels with values above 200 μK_{RJ} of the freefree amplitude at both 30′ and 2° FWHM angular resolution evaluated at 40 GHz; all pixels with values above 100 μK_{RJ} of the point source amplitudes smoothed with a 1° FWHM beam evaluated at the LFI 30 GHz band frequency; and the regions of the sky contributing to the largest 2.5% of the χ^{2}. The accepted sky fraction of the resulting mask is f_{AME} = 0.91.
6.2. Markov chain trace plots and correlations
With the data, sampling algorithms, and priors in place, we are ready to consider the actual Markov chain results. As described by BeyondPlanck Collaboration (2023), for the final analysis, we produced two independent chains, each with 750 samples.
Figure 16 shows the first 250 samples for a selection of relevant parameters. Several points are worth noting in this figure. First, we note that the burnin period for most of the foregroundinduced parameters is very short, while it is slightly longer for some of the global gain and bandpass instrumental parameters. We removed the first ten samples for burnin for the following analysis. This short burnin is a combination of the novel sampling algorithms introduced in Sect. 4 and the prior choices discussed above. At the same time, we do observe a weak shift in the average of ν_{p}, where the chains split away from each other around sample 120, which appears to trace some of the more slowly varying instrumental parameters, primarily the total bandpass correction of LFI 30 GHz; this is quite typical, as many global instrumental parameters tend to have long Markov correlation lengths and these directly affect foreground residuals.
Fig. 16. Trace plots of a set of sampled component separation parameters and selected instrument parameters, including absolute calibration (g), bandpass shift (Δ_{bp}), and the Solar CMB dipole () from a naive mono and dipole estimate using the band monopole mask, as described in Fig. 3. The parameter a_{ame} is the AME amplitude evaluated at 10° FWHM and N_{side} = 16 of the pixel at Galactic coordinates (l, b) = (31° ,5° ). Parameters from chain 1 are plotted in blue and chain 2 in orange. 
In Fig. 17, we plot Pearson’s correlation coefficients between the various parameters. For this particular plot, we have subtracted a running average with a length of ten samples (five samples prior and succeeding) from each function before evaluating the correlation coefficient, in order to highlight samplebysample correlations; two parameters may appear to be spuriously correlated if there are longterm gradients, irrespective of their origins.
Fig. 17. Correlation coefficient plot of local deviation, from a running mean of five prior and five succeeding samples of the sampled spectral parameters, component crosscorrelation intersections, frequency band monopoles, absolute gain calibration, LFI 30 GHz bandpass shift, and CMB dipole amplitude, all as described in Fig. 16. 
Several interesting observations may be made from this plot. First, we note that there is a very high correlation between β_{d} and both the AME peak frequency ν_{p} and the CMB dipole amplitude . This is not unexpected, given the critically important role of thermal dust emission at all CMB frequencies. At the same time, this also serves as a useful reminder that several important BEYONDPLANCK results depend directly on official Planck results through the use of the HFI 857 GHz frequency band and the assumed thermal dust spectral index prior of β_{d} = 1.56 ± 0.03 (Planck Collaboration X 2016; Planck Collaboration IV 2020; Planck Collaboration Int. LVII 2020), and the systematic errors in these are not propagated properly in the current analysis. Integrating HFI observations into the BEYONDPLANCK framework at the TOD level is clearly an important goal of nearterm work.
Next, we note that all microwave frequency monopoles are internally strongly correlated and, notably, anticorrelated with the component monopoles. Both of these observations make intuitive sense, as if one increases a given component monopole, the frequency monopoles have to decrease to result in a net zero change to the overall model. In addition, all frequency monopoles have to change in the same direction to a given change in the component monopoles. Accurately accounting for all of these correlations in terms of MCMC samples is perhaps the most important advantage of the Bayesian endtoend BEYONDPLANCK framework.
6.3. Goodnessoffit statistics
Before presenting the actual posterior distributions, we consider the goodnessoffit of the derived model. First, Fig. 18 shows posterior mean residual maps on the form r_{ν} = d_{ν} − s_{ν} for each of the ten frequency bands included in the analysis, where the average is evaluated over all samples in the Markov chain.
Fig. 18. Mean residual maps for the different frequency channels included in the BEYONDPLANCK analysis. All maps have been smoothed to a common angular resolution of 2° FWHM. The Planck LFI residuals are plotted with a range of 3 μK_{CMB}, and the WMAP residuals are plotted with a range of 5 μK_{CMB}. Bottom right panel: shows the mean reduced chisquared per N_{side} = 16 pixel of the BEYONDPLANCK Gibbs chain. 
Starting with the LFI channels shown in the top section of the figure, we see that the residuals at high Galactic latitudes are largely consistent with instrumental noise, except for scattered point source residuals at 30 GHz, while at low Galactic latitudes there is an obvious dust morphology residual coming from either AME or thermal dust. At 30 GHz it is also possible to see some bright negative freefree residuals. Overall, though, the fits are performing quite well, with typical residuals smaller than 3 μK over most of the sky, and the remaining artefacts are relatively easy to mask through Galactic and point source thresholding.
The middle section shows the WMAP channels, plotted on a color range of ±5 μK. The most striking residual feature in this case is a strong dust residual in the Ka (33 GHz) channel. This, combined with the weaker dust residuals observed in the LFI channels, strongly suggests that the current singlecomponent AME model adopted for the BEYONDPLANCK analysis is not a statistically adequate model for the actual AME sky. In fact, this shortcoming was already pointed out by Planck Collaboration X (2016), who introduced a second AME component to fit the full contribution. Doing the same with the current data selection would lead to a massively increased noise level for all derived components, and we instead accept the foreground mismodeling here and we instead simply make sure to mask out the contaminated regions of the sky in higherorder analyses.
Other notable features in the WMAP residuals are large regions of lowlevel residuals at high Galactic latitudes that do not obviously trace known Galactic components. As discussed by Barnes et al. (2003), an important challenge regarding this data set on large angular scales is sidelobe modeling and this may also be relevant for the residuals we see in Fig. 18. A reanalysis of the timeordered WMAP data within the BEYONDPLANCK framework is already ongoing (Watts et al. 2023).
The last two frequency channels, Haslam 408 MHz and HFI 857 GHz, show very uniform residuals. This is simply due to the fact that their unique S/N values massively dominate the synchrotron and thermal dust components, respectively, and any potential mismodeling would therefore leak directly into the component amplitude maps. A flat residual should thus not be interpreted as the absence of systematic errors, but rather simply as an indication that these two channels have no significant crosscheck in terms of their effect on the signal model by any other channel. We note that this is not entirely true for the Haslam 408 MHz map, which does have competitors in terms of freefree emission near the Galactic center in Planck and WMAP, and corresponding residuals may be seen here.
The bottom right panel of Fig. 18 shows the reduced and normalized χ^{2} per N_{side} = 16 pixel, as defined by:
where the sum runs over all pixels within a given lowresolution pixel, and n_{d.o.f.} = 15 400 is an estimate of the total number of degrees of freedom for each lowresolution pixel. Since n_{d.o.f.} is large, this quantity is expected to be N(0, 1) distributed in the ideal case, and Fig. 18 thus quantifies the agreement between the data and the model in units of standard deviations per pixel. Overall, we see that the distribution agrees with the expectation to about 1σ at high Galactic latitudes, except for some compact sources, while at low Galactic latitudes there is a strong residual with a clear dustlike morphology. This χ^{2} map serves as an important input for producing masks for higher order analyses.
6.4. Signal component posterior distributions
Next, we consider the signal component posterior distributions and we start with the frequency monopoles, as summarized in terms of posterior mean and standard deviations in the second column of Table 3. For comparison purposes, the third column shows corresponding results derived by the Planck team; all results except 857 GHz are reproduced from the Planck 2015 analysis (Planck Collaboration X 2016), while the 857 GHz result is taken from Planck DR4 (Planck Collaboration Int. LVII 2020).
Comparison of frequency channel monopole constraints from BEYONDPLANCK and Planck.
Several important differences between the two sets of results can be identified. First, we note that the BEYONDPLANCK LFI mean monopoles are all zero; this happens by construction during the mapmaking process, as the frequency band monopoles are determined directly from the sky model and any deviation from this is assigned to the correlated noise component (Ihle et al. 2023). On the other hand, we do see that the uncertainties of the LFI zerolevels are larger at both 30 and 44 GHz, reflecting the difficulty of uniquely determining the AME offset, as discussed in Sect. 4.2. In addition, the uncertainty of the Haslam, and thereby the synchrotron monopole propagates down to the lower frequency Planck and WMAP channels and contributes to increased uncertainties of these monopoles. We argue that the uncertainties determined through componentbased monopole determination are more realistic than those obtained through morphologically based frequency monopole determination.
For WMAP, we note that the our monopole corrections are considerably larger compared to those determined in the Planck 2015 estimates. This is caused by two main differences. First, the 2015 analysis adopted the WMAP Kaband explicitly as a fixed “anchor channel” (Planck Collaboration X 2016) that was not allowed to vary in the analysis. In the current analysis, for which we instead impose constraints directly on the component monopoles, this channel is instead associated with a 16 μK_{CMB} correction, which is then accounted for by an offset of −17 μK_{CMB} at LFI 30 GHz in the 2015 analysis. Second, we have changed the average synchrotron index from β_{s} = −3.1 in the Planck 2015 analysis (Planck Collaboration X 2016) to β_{s} = −3.3 in the current analysis. Since the synchrotron monopole is set to 8.9 K_{RJ} at 408 MHz (see Sect. 3), the predicted monopole at WMAP Kaband is 11 μK_{CMB} with β_{s} = −3.1 and 4 μK_{CMB} with β_{s} = −3.3, a difference of 7 μK_{CMB}. The frequency band monopole has to adjust for this difference. Similarly, differences in the zerolevel of the other components, especially the AME, will contribute to frequencyband monopole differences in equal fashion. Regarding the offsets for the 857 GHz channel, we see that these agree within 2σ as estimated by BEYONDPLANCK and the uncertainties are significantly larger (and, we believe, more realistic) in the new approach, which is a reflection of the fact that we are now propagating uncertainties in the thermal dust zerolevel as discussed in Sect. 4.2.
Regarding spectral parameters, the AME peak frequency is the only spectral parameter that is fully sampled from the data in this component separation work. For this, we find a posterior mean and standard deviation of ν_{p} = 25.3 ± 0.5 GHz. From the high correlation coefficient between ν_{p} and the thermal dust β_{d} seen in Fig. 17, it is clear that this result is highly dependent on the assumptions made regarding thermal dust modeling in this paper, and future work that includes Planck HFI observations will be critically important to make the AME model more robust.
Posterior mean and standard deviation maps for each of the four component maps (synchrotron, freefree, AME, and thermal dust emission) are shown in Figs. 19 and 20, each at their own angular resolution (120, 30, 120, and 10′) and reference frequency (408 MHz, 40 GHz, 22 GHz, and 857 GHz), and all shown in terms of brightness temperature. In Fig. 20, we notice the dark stripes following the Planck scanning pattern at high galactic latitudes in the standard deviation map of the freefree component. These stripes are barely visible in the AME standard deviation map while being completely absent from the synchrotron and thermal dust maps. This can be explained by the sampling of the amplitude zerolevels, where the freefree zerolevel is the only one not marginalized over in the BEYONDPLANCK Gibbs chain. The base value at higher latitudes of the amplitude RMS is equivalent to the RMS of the sampled zerolevels.
Fig. 19. Posterior mean amplitude maps for each of the four fitted foreground component; synchrotron (topleft), freefree (topright), AME (bottomleft), and thermal dust emission (bottomright). The angular resolutions of the four maps are 120, 30, 120, and 10′ FWHM, respectively. 
Fig. 20. Posterior standard deviations for the same maps as shown in Fig. 19. We note that synchrotron and dust emission are potted with linear scaling. 
Figure 21 shows mean and standard distribution plots for the compact source component for a 20° ×20° field of the sky in the Northern Galactic hemisphere for each of the three LFI frequency channels. The fourcomponent posterior mean maps may be compared to similar products from the Planck 2015 analysis (Planck Collaboration X 2016) and corresponding difference maps are shown in Fig. 22. All maps are smoothed to a common angular resolution of 2° FWHM before differencing, and a free offset has been fitted and subtracted using the frequencyband monopole mask discussed in Sect. 4.3. In addition, for the components that have different reference frequencies in the two analyses (i.e.., freefree, AME, and thermal dust emission), a single multiplicative factor has been fitted to take into account SED scaling differences.
Fig. 21. Partial sky plots of the mean (top) and standard deviation (bottom) amplitude of the fitted compact sources (point sources) as seen by the three Planck LFI detector bands; 30 GHz (left), 44 GHz (middle), and 70 GHz (right). The plots show gnomonic projections of a 20 × 20 degree patch of the sky centered on 90° longitude and 70° latitude, with the north galactic pole located towards the topcenter of the plots. All plots are at native angular resolution and pixelization, see Table 2 for details. 
Fig. 22. Difference maps between the component amplitude maps derived in BEYONDPLANCK and the Planck 2015 analysis, for synchrotron (topleft), freefree (topright), AME (bottomleft), and thermal dust emission (bottomright), respectively. All maps have been smoothed to a common angular resolution of 2° FWHM, a relative offset has been fitted and subtracted using the frequencyband monopole mask discussed in Sect. 4.3 and differences in reference frequencies (where relevant) have been accounted for by a single multiplicative scaling factor. 
Starting with the synchrotron case, we note that typical highlatitude differences are small compared to the overall amplitude of the Haslam 408 MHz map, typically less than 1 K_{RJ}. This is of course entirely expected, since the two analyses are both dominated by the same map. However, we do see differences at both low and high Galactic latitudes; for high latitudes we note that the current analysis explicitly models individual point sources, while in the Planck 2015 analysis there was no such component and compact sources were therefore part of the diffuse components. At low latitudes, the differences are dominated by differences in the freefree, AME, and thermal dust emission models.
For freefree emission, we see clear negative imprints of the freefree amplitude map itself at low Galactic latitudes. In this case, we note that while the Planck 2015 analysis fit the electron temperature pixelbypixel, we adopt a single constant value of T_{e} = 7000 K in the current analysis. Second, we see clear positive imprints of AME or dust around the negative freefree imprint at low Galactic latitudes, resulting from degeneracies between the component signals. Lastly, we also note that the current analysis suffers significantly with respect to freefree emission due to the absence of the Planck HFI channels, which provide both angular resolution and sensitivity to this component.
The AME component residual map is largely dominated by a dipole component aligned with the direction of the Galactic center. In this case, we note that the BEYONDPLANCK analysis estimates the absolute calibration of each Planck LFI frequency channel through a joint fit with all frequencies and this provides more robust estimates of the CMB dipole. We also see a negative freefree imprint at low Galactic latitudes similar to what we see in the freefree difference map. One key feature is the almost circular blob just above and slightly to the left of the Galactic center, which is also clearly seen as a dark blob in the mean map in Fig. 19. This coincides with a similarly bright circular blob in the freefree amplitude. A similar signal can be seen in the Planck 2015 freefree signal (Planck Collaboration X 2016), although it did not affect the AME amplitude in the same way. If we look at the residuals in Fig. 18, we notice the same blob in the LFI 30 GHz and the WMAP Ka maps, indicating a degeneracy between the modeled AME and freefree components. To break this degeneracy, either more data sets needs to be introduced or the freefree T_{e} and the AME ν_{p} needs to be sampled with spatial variance, or both. This is, however, not possible with the limited data set in this analysis and therefore needs to be left for future analyses.
Finally, for the thermal dust emission amplitude map, which is essentially defined by the Planck 857 GHz frequency map, the differences are explained very closely by the morphology of the thermal dust spectral index map β_{d}(p) presented by Planck Collaboration X (2016). This is also entirely expected, given the fact that we only model β_{d} in terms of a single spatial constant over the full sky.
7. Summary and outlook
The main goal of the current paper is to establish an efficient Monte Carlo sampling scheme for intensity foregrounds within an endtoend Bayesian CMB analysis pipeline, as implemented in the BEYONDPLANCK framework. This sampling scheme must be able to operate in both low and high signaltonoise regimes, and it must be able to incorporate both algorithmic and informative priors in a controlled and transparent manner. Degeneracies between different parameters must be explored properly and it must be possible to propagate corresponding uncertainties to higher level products.
Most of the algorithmic elements in the machinery used in this paper were developed and applied within the context of the official Planck project, as described in, for instance, Planck Collaboration X (2016), Planck Collaboration IV (2020), and Planck Collaboration Int. LVII (2020). In this paper, we have added four new algorithmic components to this machinery, namely: (1) joint amplitude and spectral index sampling, borrowing heavily from ideas already introduced by Stompor et al. (2009) and Stivoli et al. (2010); (2) componentbased monopole determination; (3) joint spectral index and monopole sampling; and (4) the application of informative spatial priors. Each of these steps significantly improve the computational efficiency and robustness of the Gibbs samplingbased Commander approach.
We stress that the current BEYONDPLANCK results are not intended to define a new stateoftheart model of the astrophysical sky. Rather, the current framework and analysis constitute a “skeleton” to which additional data sets, both legacy and future, may be added in a controlled fashion, keeping track of both instrumental and astrophysical modeling errors. As more and more data sets are added, the dependency on external priors may be gradually lifted until, hopefully, all key parameters of the model become data driven. Obviously, the single most important step towards realizing this goal is the introduction of Planck HFI TOD observations, which still define the stateoftheart in terms of fullsky CMB sensitivity to date. On a longer term, we argue that all key data sets in the community should be integrated into the model, allowing one set to break the degeneracies of others. This is the goal of the Open Source COSMOGLOBE project, and the transition from BEYONDPLANCK to COSMOGLOBE may be illustrated in Fig. 23: the top panel of this figure shows the standard deviation brightness temperature as a function of frequency for each primary intensity CMB foreground as colored bands. The BEYONDPLANCK frequency bands are indicated by vertical bars. Looking at this figure, noting the large unexplored frequency ranges, it is strikingly obvious why the current data model is significantly prior dominated; thus, more data are desperately needed. The bottom panel of Fig. 23 shows an alternative scenario, in which almost the entire frequency range is covered by including data from several past and planned sky surveys. Providing a computationally and organizationally efficient platform to make this happen is the goal of COSMOGLOBE and the BEYONDPLANCK project represents an important step towards realizing this promise.
Fig. 23. Brightness temperature RMS as a function of frequency and astrophysical component for temperature. Each component is smoothed to an angular resolution of 1° FWHM, and the lower and upper edges of each band are defined by masks covering 27 and 88% of the sky, respectively. We note that foreground RMS values decrease nearly monotonically with sky fraction, whereas the CMB RMS is independent of sky fraction, up to random variations. The vertical bands represent the frequency range of detector data, where the top panel shows the data employed in this paper and the bottom panel shows some data available for future analysis. 
Acknowledgments
We thank Prof. Pedro Ferreira and Dr. Charles Lawrence for useful suggestions, comments and discussions. We also thank the entire Planck and WMAP teams for invaluable support and discussions, and for their dedicated efforts through several decades without which this work would not be possible. The current work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement numbers 776282 (COMPET4; BEYONDPLANCK), 772253 (ERC; BITS2COSMOLOGY), and 819478 (ERC; COSMOGLOBE). In addition, the collaboration acknowledges support from ESA; ASI and INAF (Italy); NASA and DoE (USA); Tekes, Academy of Finland (grant no. 295113), CSC, and Magnus Ehrnrooth foundation (Finland); RCN (Norway; grant nos. 263011, 274990); and PRACE (EU).
References
 AliHaïmoud, Y. 2010, Astrophysics Source Code Library [record ascl:1010.016] [Google Scholar]
 Alves, M. I. R., Calabretta, M., Davies, R. D., et al. 2015, MNRAS, 450, 2025 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, C., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 97 [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 BeyondPlanck Collaboration (Andersen, K. J., et al.) 2023, A&A, 675, A1 (BeyondPlanck SI) [Google Scholar]
 Colombo, L. P. L., Eskilt, J. R., & Paradiso, S. 2023, A&A, 675, A11 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693 [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Draine, B. T., & Lazarian, A. 1998, ApJ, 494, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Erickson, W. C. 1957, ApJ, 126, 480 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., O’Dwyer, I. J., Jewell, J. B., et al. 2004, ApJS, 155, 227 [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [Google Scholar]
 Finkbeiner, D. P. 2003, ApJS, 146, 407 [Google Scholar]
 Fixsen, D. J. 2009, ApJ, 707, 916 [Google Scholar]
 Galloway, M., Andersen, K. J., Aurlien, R., et al. 2023a, A&A, 675, A3 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galloway, M., Reinecke, M., Andersen, K. J., et al. 2023b, A&A, 675, A8 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Geman, S., & Geman, D. 1984, IEEE Trans. Pattern Anal. Mach. Intell., 6, 721 [CrossRef] [Google Scholar]
 GénovaSantos, R., RubiñoMartín, J. A., Rebolo, R., et al. 2015, MNRAS, 452, 4169 [CrossRef] [Google Scholar]
 Gjerløw, E., Ihle, H. T., Galeotta, S., et al. 2023, A&A, 675, A7 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
 Gregory, P. C., Scott, W. K., Douglas, K., & Condon, J. J. 1996, ApJS, 103, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1 [NASA ADS] [Google Scholar]
 Hensley, B. S., & Draine, B. T. 2020, ApJ, 895, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Weiland, J. L., Hill, R. S., et al. 2009, ApJS, 180, 225 [Google Scholar]
 Ihle, H. T., Bersanelli, M., Franceschet, C., et al. 2023, A&A, 675, A6 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jewell, J., Levin, S., & Anderson, C. H. 2004, ApJ, 609, 1 [Google Scholar]
 King, O. G., Copley, C., Davies, R., et al. 2010, SPIE Conf. Ser., 7741, 1 [Google Scholar]
 Lenz, D., Hensley, B. S., & Doré, O. 2017, ApJ, 846, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Mather, J. C., Cheng, E. S., Cottingham, D. A., et al. 1994, ApJ, 420, 439 [Google Scholar]
 Murphy, T., Sadler, E. M., Ekers, R. D., et al. 2010, MNRAS, 402, 2403 [Google Scholar]
 Paradiso, S., Colombo, L. P. L., Andersen, K. J., et al. 2023, A&A, 675, A12 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Penzias, A. A., & Wilson, R. W. 1965, ApJ, 142, 419 [CrossRef] [Google Scholar]
 Planck Collaboration XXIV. 2011, A&A, 536, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2014, A&A, 571, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2014, A&A, 571, A5 [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2016, A&A, 594, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2016, A&A, 594, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2020, A&A, 641, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. LVII. 2020, A&A, 643, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J.F. 2011, MNRAS, 418, 467 [Google Scholar]
 Seljebotn, D. S., Bærland, T., Eriksen, H. K., Mardal, K. A., & Wehus, I. K. 2019, A&A, 627, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seljebotn, D. S., Mardal, K.A., Jewell, J. B., Eriksen, H. K., & Bull, P. 2014, ApJ, 210, 24 [NASA ADS] [Google Scholar]
 Shewchuk, J. R. 1994, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, Edition 1¼, http://www.cs.cmu.edu/~quakepapers/painlessconjugategradient.pdf [Google Scholar]
 Silsbee, K., AliHaïmoud, Y., & Hirata, C. M. 2011, MNRAS, 411, 2750 [NASA ADS] [CrossRef] [Google Scholar]
 Stivoli, F., Grain, J., Leach, S. M., et al. 2010, MNRAS, 408, 2319 [CrossRef] [Google Scholar]
 Stompor, R., Leach, S., Stivoli, F., & Baccigalupi, C. 2009, MNRAS, 392, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Svalheim, T. L., Zonca, A., Andersen, K. J., et al. 2023a, A&A, 675, A9 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Svalheim, T. L., Andersen, K. J., Aurlien, R., et al. 2023b, A&A, 675, A14 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wandelt, B. D., Larson, D. L., & Lakshminarayanan, A. 2004, Phys. Rev. D, 70, 083511 [Google Scholar]
 Watts, D. J., Galloway, M., Ihle, H. T., et al. 2023, A&A, 675, A16 (BeyondPlanck SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wehus, I. K., Fuskeland, U., Eriksen, H. K., et al. 2017, A&A, 597, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Comparison of frequency channel monopole constraints from BEYONDPLANCK and Planck.
All Figures
Fig. 1. Illustration of conditional (pure Gibbs; orange) and marginal sampling (blue) algorithms for a highly correlated (Pearson’s correlation coefficient of ρ = 0.99) twodimensional Gaussian distribution (black contours). The initial position (β_{init}, a_{init}) is indicated by a purple dot. Left: comparison of unnormalized conditional P(β ∣ a_{init}) distribution evaluated at the initial position and the corresponding marginal P(β) distribution; note that the latter is much wider than the former. Assuming β samples drawn at the 10th percentile, the graphs along the vertical lines represent unnormalized conditial distributions P(a ∣ β) evaluated at the β values drawn with the conditional (orange) and marginal (blue) distributions of β; note that the marginal sampling case results in a much longer step length between the initial and final sample values. Right: samples of a standard Gibbs sampling chain (orange) using conditional sampling for both a and β, and a sampling chain (blue) using marginal sampling for β and conditional sampling for a. Both cases show the first 100 samples initialized from the purple point. 

In the text 
Fig. 2. Comparison of AME ν_{p} chains derived using the conditional (orange; Eq. (22)) and marginal (blue; Eq. (23)) samplers discussed in the text. For the purposes of this illustration, all other parameters than a_{AME} and ν_{p} are fixed. 

In the text 
Fig. 3. Masks used for signal amplitude zerolevel and band monopole sampling: Frequencyband monopole and CMB amplitude zerolevel (top) and freefree amplitude zerolevel (bottom), with a freefree amplitude sample at an angular resolution of 30′ FWHM plotted underneath. 

In the text 
Fig. 4. Crosscorrelations between the H I column density N_{H I} and (top) the thermal dust amplitude and (bottom) the AME amplitude. The thermal dust crosscorrelation is evaluated at HEALPix resolution N_{side} = 64 and a common angular resolution of 60′ FWHM, while the AME crosscorrelation is evaluated at HEALPix resolution N_{side} = 16 and a common angular resolution of 10° FWHM. The lines represent bestfit lines for pixels with H I column densities less than 2 × 10^{20} cm^{−2} (blue) or 4 × 10^{20} cm^{−2} (orange). The green curve is the bestfit second degree polynomial to pixels with H I column densities less than 4 × 10^{20} cm^{−2}. The contour lines are plotted at 0.001, 0.01, 0.05, and 0.1 ; where only the lower three contour line values are plotted for AME. The contours have been smoothed for visualization. 

In the text 
Fig. 5. Comparison of AME ν_{p}χ^{2} distributions with (green) and without (orange) monopole marginalization. These distributions are evaluated using the marginal spectral parameter likelihood given in Eq. (23), but the same qualitative behaviour holds irrespective of which spectral parameter distribution is used (conditional, ridge, or marginal): the distribution becomes significantly wider when marginalizing over band monopoles. 

In the text 
Fig. 6. Comparison of AME (left) and freefree (right) amplitude maps derived without any spatial priors in a 20 × 20° field centered on the Galactic South Pole, (l, b) = (0° , − 90° ) at an angular resolution of 60′ FWHM. We note the striking anticorrelation between the two maps. 

In the text 
Fig. 7. Difference maps between the derived amplitude and prior maps for AME for three different 857 GHz scaling factors. From top to bottom, the three panels: scaling factors of α = 2 × 10^{−5}, α = 3 × 10^{−5}, and α = 4 × 10^{−5}. 

In the text 
Fig. 8. Effects of the spatial prior on the sampled AME amplitude. AME amplitude prior map, derived by scaling the Planck DR4 857 GHz by α = 3 × 10^{−5} and smoothing to 10′ FWHM (top). Derived AME amplitude maps for four different spatial prior combinations (bottom), . Rows: results for and , respectively, while columns show results for β = −1 and 1. 

In the text 
Fig. 9. Average per N_{side} = 16 pixel as a function of AME (top) and freefree (bottom panel) prior amplitude, q, where and n_{dof} = 15 400. Colored solid lines show results for different tilt parameters, β, and the intersection between the thick green curve and dashed line indicates the prior combination adopted for the main BEYONDPLANCK analysis. 

In the text 
Fig. 10. Freefree amplitude prior map, adopted from the Planck 2015 analysis which includes HFI observations (top). Derived freefree amplitude map for four different spatial prior combinations (bottom), . Rows show results for and , respectively, while columns show results for β = 1 and 3. 

In the text 
Fig. 11. Relative difference maps, ϵ_{i}, for each of the sky components within the simulated BEYONDPLANCK dataset. Top row: AME, CMB, and thermal dust maps. Bottom row: freefree and synchrotron maps. 

In the text 
Fig. 12. Full sky spectral parameters as a function of sample (blue) for the controlled simulation. The nominal input value of the simulation is overlayed as the black dashed line. 

In the text 
Fig. 13. Band monopoles for each of the simulated frequency bands as a function of sample. The nominal input value is given by the overlaid black dashed line. Right panel: the distribution of the samples. 

In the text 
Fig. 14. AME (left) and freefree (right) amplitude maps derived for two different synchrotron spectral indices. Top row: results for β_{s} = −2.8, while the bottom row shows results for β_{s} = −3.3. In the top row, we see a negative synchrotronlike (see Fig. 19) imprint, which is not present in the bottom row. 

In the text 
Fig. 15. χ^{2} distributions (blue curves) from coarse grid evaluation of each free spectral parameters. From top to bottom, the panels show (1) AME peak frequency, ν_{p}; (2) synchrotron spectral index, β_{s}; and (3) freefree electron temperature, T_{e}. The minimum χ^{2} value for each parameter has been subtracted in each case. Orange curves show the priors adopted for the given component; see Table 1. The prior values for β_{s} and ν_{p} have been scaled by a factor of 100 to fit in the plot with the derived χ^{2} values. 

In the text 
Fig. 16. Trace plots of a set of sampled component separation parameters and selected instrument parameters, including absolute calibration (g), bandpass shift (Δ_{bp}), and the Solar CMB dipole () from a naive mono and dipole estimate using the band monopole mask, as described in Fig. 3. The parameter a_{ame} is the AME amplitude evaluated at 10° FWHM and N_{side} = 16 of the pixel at Galactic coordinates (l, b) = (31° ,5° ). Parameters from chain 1 are plotted in blue and chain 2 in orange. 

In the text 
Fig. 17. Correlation coefficient plot of local deviation, from a running mean of five prior and five succeeding samples of the sampled spectral parameters, component crosscorrelation intersections, frequency band monopoles, absolute gain calibration, LFI 30 GHz bandpass shift, and CMB dipole amplitude, all as described in Fig. 16. 

In the text 
Fig. 18. Mean residual maps for the different frequency channels included in the BEYONDPLANCK analysis. All maps have been smoothed to a common angular resolution of 2° FWHM. The Planck LFI residuals are plotted with a range of 3 μK_{CMB}, and the WMAP residuals are plotted with a range of 5 μK_{CMB}. Bottom right panel: shows the mean reduced chisquared per N_{side} = 16 pixel of the BEYONDPLANCK Gibbs chain. 

In the text 
Fig. 19. Posterior mean amplitude maps for each of the four fitted foreground component; synchrotron (topleft), freefree (topright), AME (bottomleft), and thermal dust emission (bottomright). The angular resolutions of the four maps are 120, 30, 120, and 10′ FWHM, respectively. 

In the text 
Fig. 20. Posterior standard deviations for the same maps as shown in Fig. 19. We note that synchrotron and dust emission are potted with linear scaling. 

In the text 
Fig. 21. Partial sky plots of the mean (top) and standard deviation (bottom) amplitude of the fitted compact sources (point sources) as seen by the three Planck LFI detector bands; 30 GHz (left), 44 GHz (middle), and 70 GHz (right). The plots show gnomonic projections of a 20 × 20 degree patch of the sky centered on 90° longitude and 70° latitude, with the north galactic pole located towards the topcenter of the plots. All plots are at native angular resolution and pixelization, see Table 2 for details. 

In the text 
Fig. 22. Difference maps between the component amplitude maps derived in BEYONDPLANCK and the Planck 2015 analysis, for synchrotron (topleft), freefree (topright), AME (bottomleft), and thermal dust emission (bottomright), respectively. All maps have been smoothed to a common angular resolution of 2° FWHM, a relative offset has been fitted and subtracted using the frequencyband monopole mask discussed in Sect. 4.3 and differences in reference frequencies (where relevant) have been accounted for by a single multiplicative scaling factor. 

In the text 
Fig. 23. Brightness temperature RMS as a function of frequency and astrophysical component for temperature. Each component is smoothed to an angular resolution of 1° FWHM, and the lower and upper edges of each band are defined by masks covering 27 and 88% of the sky, respectively. We note that foreground RMS values decrease nearly monotonically with sky fraction, whereas the CMB RMS is independent of sky fraction, up to random variations. The vertical bands represent the frequency range of detector data, where the top panel shows the data employed in this paper and the bottom panel shows some data available for future analysis. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.