© The Authors 2024

1 Introduction

The Galactic interstellar medium (ISM), which represents 1015% of the total mass of the Milky Way, is constituted of diffuse material in the Galactic plane and Galactic halo. Among its various components, cold molecular gas is concentrated in big molecular clouds. In such clouds, carbon monoxide (CO) is abundant, and rotational lines of CO emission, at multiples of v = 115.27 GHz, provide a convenient observable for tracing the distribution of cold molecular gas in the Milky Way. With the lines of emission being located in frequency bands commonly used for cosmic microwave background (CMB) observations, CO also is a foreground contaminant that must be cleaned-out for sensitive CMB observations. To that effect, reliable, high signal-to-noise-ratio, full-sky tracers of CO emission would be useful for the preparation and analysis of upcoming sensitive CMB survey data over large fractions of high-Galactic-latitude sky.

Numerous ground-based observations of CO line emission have been performed in the past decades. Notably, a large survey of the CO J = 1 → 0 line emission at low Galactic latitude was published by Dame et al. (2001), with a recent revision to extend the footprint of the survey (Dame & Thaddeus 2022). While many other observations exist on selected lines of sight or in limited regions of specific interest (Hartmann et al. 1998; Ikeda et al. 1999; Fixsen et al. 1999; Magnani et al. 2000; Sawada et al. 2001; Mizuno & Fukui 2004; Wilson et al. 2005; Jackson et al. 2006; Bieging et al. 2010; Yoda et al. 2010; Bieging & Peters 2011; Handa et al. 2012; Oka et al. 2012; Polychroni et al. 2012; Burton et al. 2013; Dempsey et al. 2013; Rigby et al. 2016; Colombo et al. 2019; Benedettini et al. 2020; Rennie et al. 2022; Dong et al. 2023; Park et al. 2023), no spectroscopic survey of the full sky Galactic CO emission in various rotational lines exists as of now.

The Planck space mission, launched by the European Space Agency (ESA) in 2009, was not designed to map CO emission. However, the excellent sensitivity of the Planck detectors allowed for the release of full-sky observations of CO emission in the J =1 → 0, J = 2 → 1, and J =3 → 2 rotational lines (Planck Collaboration XIII 2014). Two main data products have been released, dubbed TYPE 1 and TYPE 2 CO maps. The TYPE 1 maps make use of the different sensitivity to CO of detectors nominally in the same Planck frequency band. TYPE 2 maps assume a model of dust emission and fit, in each sky pixel, for a multicomponent model of Galactic emission that includes CO lines, the CMB, and thermal dust. An additional data product, the TYPE 3 CO map, assumes a constant line ratio between the various lines and solves for a single CO map in a multicomponent fit than maximizes the output CO map signal-to-noise ratio. CO products obtained with the Commander component separation method have also been obtained at HEALPix¹ N_side = 256 for the three lines, and with 7.5′ angular resolution and N_side = 2048 for the J = 2 → 1 line only, the latter assuming a fixed line ratio for the first three rotational lines (Planck Collaboration X 2016).

Published Planck CO maps suffer from significant limitations. The TYPE 1 maps, which use only single-frequency data to extract the CO, are in principle less prone to sky modeling errors, but the signal-to-noise ratio is low, in particular at high Galactic latitude where CO emission is faint. TYPE 2 maps, on the other hand, have reduced noise level, but are more prone to systematic errors due to variations of real sky emission scaling compared to the model that is assumed, and they are available only for the first two CO lines (J = 1 → 0 and J=2 → 1). TYPE 2 maps are also delivered at lower angular resolution (15′ for the two lines) instead of the native resolution of the Planck frequency channels, as is the case for TYPE 1 maps. The TYPE 3 map, which assumes a fixed CO line ratio, cannot be trusted in regions of strong CO emission, where variations of the line ratio are clearly detected. Finally, the Commander maps are either at low HEALPix N_side or not available for all lines, and they are also not obtained in a pipeline specifically optimized for mapping CO.

In this work, we aim to extract improved maps of CO line emission for the three lines detected by Planck (J = 1 → 0, J = 2 → 1, and J =3 → 2), with reduced noise contamination, and taking into account the variation of the line ratio across the sky. To that effect, we post-process some of the published Planck CO map products using a modified version of the generalized needlet Internal Linear Combination (GNILC) algorithm of Remazeilles et al. (2011), in which we assume increasing priors as the signal-to-noise ratio decreases in the needlet domains considered. This approach also produces maps at fixed high angular resolution, instead of the variable angular resolution that is obtained by the standard GNILC algorithm².

The rest of this paper is organized as follows. In Sect. 2, we discuss and compare available Planck CO products and their limitations. Section 3 presents the pipeline we used to produce reprocessed full-sky maps of the first three rotational lines of CO emission. Data products are presented in Sect. 4 and discussed in Sect. 5. We conclude in Sect. 6.

2 Existing Planck CO maps

In this work, we only considered the Planck legacy data products for the CO line emissions. Here, we discuss these legacy data products and their limitations.

2.1 Single-frequency TYPE 1 maps

The Planck TYPE 1 maps were generated for each CO line by forming a weighted linear combination of difference maps between the various pairs of bolometers associated with the single-frequency channel corresponding to the line’s frequency using the Modified Internal Linear Combination Algorithm (MILCA; Hurier et al. 2013; Planck Collaboration XIII 2014). Taking into account that CMB and dust emissions do not vary significantly across bolometers within a specific frequency channel, whereas the CO line response depends much on the spectral bandpass of each bolometer, the MILCA weights were selected to deproject CMB and dust emissions (assuming both exhibit a uniform spectrum across the bolometers in thermodynamic units), while preserving the CO signal.

The Planck TYPE 1 CO maps consist of three full-sky maps, corresponding to the rotational transition lines J = 1 → 0, J = 2 → 1 and J =3 → 2. The angular resolution and HEALPix N_side of the maps are those of the original Planck frequency maps, per frequency channel. The single-channel methodology employed for TYPE 1 CO maps effectively reduces the potential contamination from thermal dust, which could otherwise occur with multifrequency combinations. However, it does come with the cost of a lower signal-to-noise ratio in contrast.

2.2 Multifrequency TYPE 2 maps

The TYPE 2 maps were obtained by applying a multicompo-nent generalized least-squares (GLS) filter (Ruler algorithm; Planck Collaboration XII 2014; Planck Collaboration XIII 2014) to multifrequency Planck observations to disentangle CO emission from other astrophysical foregrounds. This requires the prior assumption of a parametric model for Galactic emission used in the GLS mixing matrix.

The Planck TYPE 2 CO maps consist of only two full-sky maps, corresponding to the lines J = 1 → 0 and J = 2 → 1. Both maps are produced at a common angular resolution of 15′, and HEALPix N_side of 2048. The multichannel approach enhances the signal-to-noise ratio for TYPE 2 CO maps in comparison to TYPE 1 CO maps. However, it is more prone to residual foreground contamination, stemming from thermal dust, free-free, and radio source emissions that are partially projected into the final products by reason of potential modeling errors across the frequency channels.

Specifically, TYPE 2 CO maps are obtained by assuming negligible CO emission (as compared to dust) in the 353 GHz Planck channel, fitting a linear mixture of CO, dust, CMB and free-free emission in each pixel, with fixed emission laws for each of the extra components besides CO (for details, see Sect. 2.2 of planck Collaboration XIII 2014). For free-free emission, the spectral index is assumed to be α_ff = −2.15 in brightness units, while the dust model assumes uniform dust temperature (T_dust = 17 K) and spectral index (β_dust = 1.6). Hence, the presence of a strong steep spectrum (α_source < −2.15) radio source would result in an over-estimation of the low-frequency foreground at the line frequency and, thus, an underestimation of the CO emission. Similarly, dust emission with spectral parameters very different from what is assumed will result in positive or negative leakage of dust emission in the TYPE 2 maps.

2.3 Commander maps

A J = 2 → 1 map is obtained with Commander at resolution 7.5′ and N_Side = 204 (Planck Collaboration X 2016). This map results from a global fit of all astrophysical components across Planck frequency channels. As such, it is more prone to modeling errors generating confusion between components than the TYPE 1 products. It is challenging to make a component separation pipeline that is optimal for all of the components present in the observations. The Commander method is quite flexible, but its implementation in planck Collaboration X (2016) was not specifically tuned for a reliable mapping of CO emission lines in the Planck data sets. Maps obtained with Commander at N_side = 256, which do not fully exploit the angular resolution of the Planck observations, are not used in the reprocessing performed here.

2.4 Limitations of existing maps

Figure 1 plots TYPE 2 versus TYPE 1 map values for all sky directions, showing a clear excess of TYPE 2 J =1 → 0 and deficit of TYPE 2 J=2 → 1 emission as compared to TYPE 1 products. This discrepancy is clearly visible for regions with large CO emission overall, but it is not obviously seen for lower emission pixels.

This calibration mismatch is also visible in a comparison of the power spectra of the TYPE 1 and TYPE 2 maps (Fig. 2). The discrepancy between the TYPE 1 and TYPE 2 maps is reduced when the power spectra are computed after masking the brightest sky regions using a small mask of the Galactic ridge and of the brightest CO regions, displayed in the top panel of Fig. 3. This is not the case, however, for the discrepancy between the Commander map and the other two, which is increased when the brightest regions are masked.

As discussed in Planck Collaboration XIII (2014), discrepancies between the TYPE 1 and TYPE 2 maps likely originate from a mixture of contamination by¹³ CO and other diffuse Galactic emission residuals, and local inadequacy of the assumed emission spectral parameters for the generation of the TYPE 2 data products. Residual contamination from thermal Sunyaev-Zeldovich (SZ) emission is also evident in the TYPE 2 maps, with massive galaxy clusters such as Coma clearly visible in these maps, but not in TYPE 1 maps. Since the TYPE 1 data products are obtained from linear combinations of observations in the same frequency band, for which the contribution of all Galactic emissions other than CO lines are canceled, TYPE 1 maps are more reliable tracers of the actual integrated CO emission (a linear combination of the main line and isotopo-logues), and less prone to contamination by other Galactic and extragalactic components.

The Commander map seems to lack power in low signal regions, and a visual inspection also shows significant contamination by SZ emission from galaxy clusters. As in addition the noise level is higher than that of the TYPE 2 map, the Commander map it is not used for the present work, as it brings little information that is not already present in TYPE 1 and TYPE 2 maps.

Only the TYPE 1 map is available for the J =3 → 2 emission line. The central part of the map is displayed in Fig. 4. A visual inspection shows large scale ringing around the Galactic plane, and residuals of scanning artifacts at higher Galactic latitude. Using comparisons with other maps, an attempt will be made to filter some of these systematic effects in our improved Planck CO data products.

Fig. 1 Scatter plot with pixel values of TYPE 1 CO map as abscissa and pixel values of TYPE 2 CO maps as ordinate. We smooth all maps to 30 arcmin resolution and increase the HEALPix pixel size to Nside = 256 for this figure. The top figure shows the scatter for CO J = 1 → 0 Planck maps, and the bottom figure shows it for CO J = 2 → 1 maps. For large signals, the signal in TYPE 2 maps is higher than in TYPE 1 maps for the J = 1 → 0 line, and lower for the J = 2 → 1 line. In the inset window, we show that the low-brightness pixels do not show obvious inconsistency. This is in agreement with the comparison made in Planck Collaboration XIII (2014).

Fig. 2 Power spectra comparison of TYPE 1 and TYPE 2 CO J =1 → 0 maps (top) and J = 2 → 1 maps (bottom), with and without an attenuation mask to reduce the amplitude of emission in the brightest 2% of pixels (shown in the top panel of Fig. 3). The TYPE 1 J =1 → 0 have lower power compared to the TYPE 2 map. For the CO J = 2 → 1 line emission, the TYPE 1 map has higher power compared to the TYPE 2 map. Both of them have higher power than the Commander CO J = 2 → 1 map. The discrepancy between TYPE 1 and TYPE 2 CO J = 2 → 1 power spectra largely goes away when masking the regions of strongest CO.

Fig. 3 Masks used in this work. Top: mask used to compute power spectrum after excluding the brightest CO emission regions in the Galactic ridge. Bottom: mask used for fitting global CO mixing coefficients at low Galactic latitude the Galactic plane, obtained from the sky fraction observed in the Dame et al. (2001) survey, excluding regions with CO emission larger than 3 KRJ km s−1. This results in the brightest 0.5% of pixels (in the Galactic ridge) being masked.

Fig. 4 TYPE 1 CO J = 3 → 2 data product from Planck, shown here after smoothing to 30′ resolution. We are showing a region ±90° in Galactic longitude and ±30° in Galactic latitude about the Galactic center. Ringing effects and systematic residuals are visible above and below the Galactic ridge.

3 Reprocessing pipeline

The different approaches that have been pursued to obtain TYPE 1, TYPE 2, TYPE 3, and Commander CO maps have advantages and drawbacks. As argued in the Planck CO paper (Planck Collaboration XIII 2014), TYPE 1 maps are more reliable (in terms of separating CO from other emissions), but more noisy than the other products. This motivates an attempt to denoise those TYpE 1 maps to improve their signal-to-noise ratio. To that effect, one can exploit the fact that the signal part in the TYpE 1 CO maps for the various CO emission lines should be strongly correlated between maps, while the noise should be essentially independent. Similarly, the higher signal-to-noise ratio TYpE 2 maps can be used to help identify CO emission clouds locally in the TYPE 1 maps.

3.1 Data model and methodology

For each of the TYPE 1 and TYPE 2 maps (a total of five maps indexed by α = 1,…, 5), we model the data d_α(p) as $$d_{\alpha }(p)=s_{\alpha }(p)+n_{\alpha }(p),$$(1)

where p denotes sky pixel, s_α(p) is the total sky signal in map α (mostly CO line emission, but possibly also some residuals of continuum Galactic emission, in particular for the TYPE 2 maps), and n_α(p) denotes the noise contamination. In vector format, this can be recast, for each pixel, as $$d(p)=s(p)+n(p),$$(2)

where all vectors are five-dimensional. In the discussion that follows, we do not explicitly write the pixel dependence, and unless explicitly stated all data vectors and matrices have a pixel dependence that is hidden.

The signal, s, which is mostly due to emission from the same CO molecular clouds, is strongly correlated between these maps. The key idea for denoising is that we can write the signal in each map as a linear combination of N_s templates (which may or may not be due solely to CO emission), t, that trace the sky signal: $$s_{\alpha }={\displaystyle \sum _{i=1}^{N_{\text{s}}}{A}_{\alpha i}{t}_i},$$(3)

or $$s=At.$$(4)

In general, for N_ch input maps and N_s signal templates, A denotes a N_ch × N_s mixing matrix. If N_s is strictly lower than N_ch (which is expected to be the case, considering that the signal is highly correlated among the various CO input maps), we can project the data on the subspace spanned by the N_s templates that are required to model the signal, and hence get rid of part of the noise contamination. Doing this with simultaneous localization in a harmonic space and pixel domain is the main idea behind GNILC (Remazeilles et al. 2011), which we used for the present work, but with a modification that we describe next³.

For some applications, one of the drawbacks of the standard GNILC procedure is that when the dimension N_s of the foreground subspace is zero in some needlet domains (i.e., in some needlet bands and in some regions of sky), the data for that needlet domain is completely discarded (projection onto a zero-dimensional subspace). This results in varying angular resolution over the sky, which complicates the map characterization and its scientific exploitation. This can be alleviated by imposing that the dimension of the foreground subspace be always at least equal to one, and more if the data requires so. If, however, one keeps the dimension corresponding to the largest eigenvalue of the covariance matrix, in regions of low signal one would keep the dimension with the highest noise excursion. Hence, in practice, when the GNILC procedure would suggest discarding the observations as too noisy (i.e., projecting on a zero-dimensional subspace), we instead used a one-dimensional prior for the signal subspace and combined all observations using a least-squares estimate of the signal assuming the corresponding signal color. This amounts to co-adding all observations with weights that maximize the signal-to-noise ratio for a pre-determined signal color, and then rescaling the resulting single CO map to inject its contribution in each individual line map.

3.2 Overview of the pipeline

The main steps of the processing pipeline can be summarized as follows.

1. Preprocessing: the input data to the pipeline are the set of five published Planck TYPE 1 and TYPE 2 CO maps, smoothed to a common 10′ angular resolution. We masked or saturated some of the brightest pixels to avoid ringing effects polluting the neighborhood when we perform harmonic transforms for the needlet decomposition; We also masked in TYPE 2 maps 186 bright galaxy clusters, the thermal SZ emission of which leaked in the Planck data products (this leakage is clearly visible for the Coma cluster, for instance). We injected noise from null maps (i.e., half-difference of half-ring Planck CO maps) in the masked regions to preserve the local statistics of the Planck maps noise.

2. Calibration: we computed a best-estimate ‘CO mixing vector’ for the set of five maps, by calibrating each map against the J = 1 — 0 map of Dame et al. (2001), reprojected into a Healpix map and available on the LAMBDA website. This provides us with a prior for the signal subspace (i.e. the five coefficients a_α = A_α1 for one single template in Eq. (3)) in needlet domains where no signal is detected above the noise in the auto- and cross-covariance of the input maps.

3. Extended GNILC pipeline (named xGNILC): we computed needlet coefficients for all five input CO maps and implemented a modified GNILC algorithm, in which a GLS solution assuming the CO scalings from the previous step is used in needlet domains where the Akaike information criterion (AIC) suggests a zero-dimensional GNILC signal subspace. We obtained output maps for each of the TYPE 1 and TYPE 2 line maps. Only outputs corresponding to the TYPE 1 inputs were kept in our final set of CO products. For all lines, the signal that had been subtracted from the input Planck CO maps in the masking or saturation process in step 1 was re-injected in the output maps.

4. Postprocessing: we postprocessed the output xGNILC J =3 → 2 map to subtract large-scale systematics originating from the input Planck TYPE 1 J =3 → 2 product.

All three maps were processed through the first three steps of the pipeline. Only the J = 3 → 2 map was processed with the additional postprocessing step. Details about each of these steps are provided in the next subsections.

3.3 Preprocessing

Although in principle the J = 2 → 1 and J =3 → 2 maps allow for an angular resolution on the order of 5^', the signal-to-noise ratio becomes very low on the smallest scales. The output data products will be maps at N_side of 1024, with ℓ_max of 2048. We choose a common resolution of 10′, but we ‘apodize’ the Gaussian beam to smoothly go to zero by ℓ_max of 2048, as shown in Fig. 5. In the text that follows, we denote this apodized beam simply as 10′ beam. Hence, as a first step, we smoothed all maps to a common 10′ angular resolution (very little information can be extracted from the observations on smaller scales in practice, and mostly in Galactic regions where other data from spectroscopic observations of the CO lines exist).

Next, we note that the Planck CO maps are very bright toward the Galactic center and toward a few compact regions of emission, mostly near the Galactic plane. These small regions of bright emission tend to produce ringing features when maps are decomposed into needlets. To avoid this unwanted feature, we identified these regions and masked and filled them. The exact process consisted of the following steps.

1. Identified the region of very bright CO emission and determine a brightness threshold in K_RJ km s⁻¹ that was required to create a clipping mask for the region.

2. Masked and set the masked pixels to the threshold value. Performed a median filter on this map, with the median computed over superpixels at N_side = 256, and smoothed the median-filtered map with a 15′ beam.

3. Used a 1° apodized mask to mask the original map and fill the gaps with the smoothed median-filtered map obtained above.

We saved the difference between the maps obtained in this way and the initial maps, for re-injection of these strong signals in the final data products.

The Planck TYPE 2 maps are contaminated by thermal SZ (tSZ) emission. To produce a mask of tSZ emission we used a simulation of known tSZ sources from a recent version (v2.3.2) of the Planck Sky Model (Delabrouille et al. 2013). We produced a mask where we set the 6.5 × 10⁻³ % brightest pixels of the simulated tSZ map to zero. We were careful not to mask anything with strong CO signal, and to that effect exclude from the masked pixels any pixel in the Galactic plane where either the Dame CO map has a value greater than 0.5 K_RJ km s⁻¹, or where in the faint clouds away from the Galactic plane the Planck TYPE 2 CO J = 2 → 1 map was larger than 0.25 K_RJ km s⁻¹. This excluded any masking that may be too close to any region containing strong CO signal⁴. We additionally masked circular regions with radii of 0.35° and 0.55° around the center of the Virgo and Coma clusters, respectively. This mask was finally apodized with a 0.6° cosine taper. We masked and filled the Planck TYPE 2 CO maps with the corresponding Planck TYPE 2 null maps using the apodized tSZ emission mask.

3.4 Calibration

Calibration against the Dame et al. (2001) map provides us with estimates of the ratio of CO emission in each Planck CO map to the CO emission in the Dame map. To perform this calibration, we first smoothed all maps to a common 30' angular resolution and downgraded them to HEALPix N_side = 512. We also generated, for all Planck products, corresponding statistical uncertainty maps downgraded to N_side = 512. We then considered two options for performing a global calibration.

1. In option 1, we selected a sky region in the Dame map using the mask displayed in the bottom panel of Fig. 3. For each Planck CO map, we then computed a global best fit for a multiplicative factor ${\overline{a}}_{\alpha }$ by sampling a Gaussian likelihood with an uniform prior of 0 ≤ a_α ≤ 2.5. For this fit, we assumed a model in which ${\text{CO}}_{\text{TYPE\hspace{0.17em}X\hspace{0.17em}}}^J=a_{\alpha }\times {\text{CO}}_{\text{DAME}}^{J=1\to 0}$, where X ∈ {1, 2} and J indexes the CO line.

2. For option 2, we repeated the same procedure, but using a fitting model with an offset, ${\text{CO}}_{\text{TYPE\hspace{0.17em}X}}^J=a_{\alpha }\times {\text{CO}}_{\text{DAME}}^{J=1\to 0}+K_{\text{TYPE\hspace{0.17em}X\hspace{0.17em}}}^J$. For the offset we assumed a flat prior, $K_{\text{TYPE\hspace{0.17em}X\hspace{0.17em}}}^J\in [-10,10]K_{\text{RJ}}\text{km\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}$.

We also investigate how the calibration varies across the sky. To that effect, we repeated the calibration after dividing the calibration region in subregions of equal sky area, obtained by sorting first the pixels used in the calibration in three equal-area Galactic latitude regions, and then each latitude region in ten equal-area longitude regions (for a total of 30 subregions). We performed the calibration independently in each of those regions. We then computed the mean, the median, and the standard deviation of the 30 calibration coefficients obtained in this way.

All calibration results are reported in Table 1, where they are also compared to expected values derived from the calibrations performed in the original Planck publication. We find reasonable consistency between the various results and choose the coefficients obtained as the mean of 30 calibrations in different regions to our fiducial a_α (from the fit with an offset). We use those values for the GLS solution in needlet domains of faint CO emission in the extended GNILC pipeline described later.

As a final investigation of the calibration of Planck CO products relative to the J =1 → 0 Dame map, we repeated the calibration locally over the entire sky, in local superpixels corresponding to HEALPix N_side = 16. For these local fits, performed for all sky superpixels where Dame data was present, we assumed, for a_α, a Gaussian prior with mean and standard deviation values from the 30-regions analysis. Resulting calibration maps are displayed in Fig. 6.

Fig. 5 Needlet bandpass windows, $h_{\ell }^j$, showing the selection weights in harmonic space (plotted in black). The “apodized” beam transfer function of the output map resolution of 10′ is also shown in red. The end of the beam transfer function is modified to smoothly go to zero at ℓmax = 2048.

3.5 An extended GNILC pipeline

The extended GNILC pipeline (dubbed xGNILC in the following) used in the present work has four main parts performed consecutively: (1) the needlet decomposition; (2) the determination of the CO subspace (in each sky pixel, for each needlet scale); (3) the projection of the input data into final needlet coefficient maps using either a multidimensional ILC or a GLS solution (the latter where the dimension of the CO subspace is estimated to be null, i.e., N_s = 0, in step 2); and (4) the recombination of the needlet maps into five processed maps. We describe each of these steps in more detail below.

3.5.1 Needlet decomposition

The so-called needlets are a kind of spherical wavelet, which allow us to compute statistics and perform map combinations locally in harmonic and in pixel spaces simultaneously (Narcowich et al. 2006). The originally introduced needlets have spectral band shapes specifically tailored for good localization. In practice, applications for CMB data analysis (e.g. Pietrobon et al. 2008; Faÿ et al. 2008; Delabrouille et al. 2009; Geller et al. 2008; Scodeller et al. 2011) have used various types of needlet bandpass windows, adapted to the angular resolution of the observations and the variability of the signals of interest and of the noise across the sky and in harmonic domain. For this work, we use needlets that are defined in harmonic space with eight Gaussian-shaped bandpasses $h_{\ell }^j$, such that ${\displaystyle {\sum }_j{\left(h_{\ell }^j\right)}^2}=1$. These bandpasses are written as $$h_{\ell }^j=\{\begin{array}{cc}\exp \left[-\frac{\ell (\ell +1){\theta }_1^2}{16\ln 2}\right]\hfill & j=1\hfill \\ \sqrt{\exp \left[-\frac{\ell (\ell +1){\theta }_j^2}{8\ln 2}\right]-\exp \left[-\frac{\ell (\ell +1){\theta }_{j-1}^2}{8\ln 2}\right]}\hfill & 1<j<N_{\text{bands\hspace{0.17em}}}\hfill \\ \sqrt{1-\exp \left[-\frac{\ell (\ell +1){\theta }_{N_{\text{bands\hspace{0.17em}}}-1}^2}{8\ln 2}\right]}\hfill & j=N_{\text{bands\hspace{0.17em}}}\hfill \end{array},$$(5)

where θ_j· is the full width at half maximum of a Gaussian smoothing function associated with each needlet band (in radians), and N_bands is the number of needlet scales. In this work, the choice of θ_j (in arcminutes) is {600′, 300′, 120′, 60′, 30′, 15′, 10′}. The harmonic windows obtained in this way are shown in Fig. 5. Using the definition for hį we can write the needlet basis function as $${\psi }_{jk}(p)=\sqrt{\frac{4\pi }{N_j}}\underset{\ell m}{{{\displaystyle \sum }}^{\text{}}}{h}_{\ell }^j{Y}_{\ell m}(p)Y_{\ell m}^{\ast }\left({\widehat{n}}_{jk}\right).$$(6)

In the HEALPix pixelization scheme, N_j is the total number of map pixels for jth needlet scale, and ${\widehat{n}}_{jk}$ gives the direction associated with the kth pixel of the jth needlet map. The data maps d, were transformed in terms of ψ_jk, giving the following coefficient maps: $$b_j\left({\widehat{n}}_{jk}\right)={{\displaystyle \int }}^{\text{}}d{\psi }_{jk}^{\ast }\text{d}{\Omega }_p.$$(7)

Finally, we computed the localized N_ch × N_ch covariance matrices from the needlet coefficient maps as $${\widehat{C}}_{jk}=\frac{1}{N_k}\underset{k^{\prime }}{{{\displaystyle \sum }}^{\text{}}}{w}_j\left(k,k^{\prime }\right){\mathbf{b}}_j\left({\widehat{n}}_{j{k}^{\prime }}\right){\mathbf{b}}_j\left({\widehat{n}}_{j{k}^{\prime }}\right),$$(8)

with the w_j(k, k′) function selecting the domain of N_k pixels around the kth pixel to compute the covariance at the kth pixel.

Fig. 6 Maps of calibration obtained by fitting Planck TYPE 1 (top row, from J = 1 → 0 to J = 3 → 2, left to right) and TYPE 2 (bottom row, from J =1 → 0 to J = 2 → 1, left to right) maps against Dame J = 1 → 0 map in Nside = 16 super pixels. We perform a fit in any super pixel if there are 50 or more original pixels at Nside = 512 with a high signal-to-noise ratio in the Dame map. For all super-pixels where this condition is not met, we assume the mean calibration given in Table 1. These maps illustrate the variation over the sky of the ratio between the CO detected in Planck TYPE 1 and TYPE 2 maps and the CO detected in the Dame J = 1 → 0 map.

3.5.2 Determination of the CO subspace

The signal model introduced in Eq. (3) assumes that we can represent the signal in the set of input maps with a mixing matrix A and templates t. The number of emission templates required to model the signal is determined by the dimension of the “signal subspace”. Multiple dimensions are needed to describe the signal when it is bright, complex, possibly multicomponent, and observed with high signal-to-noise ratio, as is the case in the Galactic plane, and in particular close to the Galactic center. One of the key features of GNILC is to determine the dimensions of the signal subspace for every needlet band as a function of sky position. This allows us to retain the differences between maps only where they are significant (i.e., detectable above the noise), but average down the noise contamination when the signal-to-noise ratio in regions of the pixel or harmonic space is low, and detectable astrophysical emissions nearly identical in the various observations (within uncertainties).

For the data model introduced for the Planck CO maps we can write the total covariance matrix $$C=C_s+C_n$$(9)

as the sum of the covariances of the signal and of the noise. We use the public Planck CO null maps from the Planck Legacy Archive⁵ to obtain an estimate for the noise covariance C_n. We can then “whiten” the total covariance as $$\tilde{C}=C_n^{-1/2}C{C}_n^{-1/2}=C_n^{-1/2}{C}_s{C}_n^{-1/2}+I.$$(10)

In general, we can then obtain a eigenvalue decomposition of the $\tilde{C}$ as $$\tilde{C}=UD{U}^{\text{T}}=U_s{D}_s{U}_s^{\text{T}}+U_n{U}_n^{\text{T}},$$(11)

where we use Eq. (10). Here U_s is a N_ch × N_s matrix with N_s columns representing eigenvectors for the eigenvalues associated with the signal subspace. As discussed in Planck Collaboration Int. XLVIII (2016), we can obtain the number of dimensions, N_s, of the signal subspace from the eigenvalues, D = diag$\left({\lambda }_1,\dots ,{\lambda }_{\alpha },\dots ,{\lambda }_{N_{\text{ch}}}\right)$, by minimizing the AIC: $$\underset{\alpha \in \left[1,N_{\text{ch}}\right]}{\min }\left[2\alpha +{\displaystyle \sum _{k=\alpha +1}^{N_{ch}}\left({\lambda }_k-\log {\lambda }_k-1\right)}\right].$$(12)

This gives an estimate of the CO subspace dimension, N_s , as a function of sky position for every needlet scale. In Fig. 7, we show the GNILC estimated CO subspace dimension for the 60′ and 10′ needlet scales. It is clear from the figure that for needlet scales covering smaller scales, the signal-to-noise ratio becomes poorer further away from the Galactic plane. Consequently, the traditional GNILC method sets the CO signal subspace dimension to be zero in these noise-dominated regions of needlet space.

The signal covariance matrix can be written as $$C_{\text{s}}=A\langle t{t}^{\text{T}}\rangle A^{\text{T}}.$$(13)

Using Eqs. (10) and (11), we can also write the signal covariance matrix as $$C_{\text{s}}=C_{\text{n}}^{1/2}{U}_{\text{s}}\left(D_{\text{s}}-I\right)U_{\text{s}}^{\text{T}}{C}_{\text{n}}^{1/2}.$$(14)

We then estimate the mixing matrix as $\widehat{A}=C_{\text{n}}^{1/2}{U}_{\text{s}}$.

Fig. 7 Maps of dimension of CO subspace as determined by GNILC for the 60′ needlet band (top) and the 10′ needlet band (bottom). At smaller angular scales, larger parts of the high-latitude sky becomes noise-dominated (dark blue in the figure). The standard GNILC algorithm assigns those regions a zero-dimensional signal subspace, which results in varying effective resolution in standard GNILC map. In our processing, in those needlet domains, we produce the CO maps using the GLS weights. With this modification to the GNILC algorithm, final xGNILC maps are produced with a uniform angular resolution.

3.5.3 Multidimensional ILC and GLS

We produced cleaned maps, ŝ, as a combination of N_ch input maps, d, with a set of optimization weights from our extended GNILC pipeline, W_xGNILC, given by $$\widehat{s}={\text{W}}_{\text{xGNILC}}d\text{.\hspace{0.17em}}$$(15)

Remazeilles et al. (2011) showed that the multidimensional ILC weights are $$W_{\text{ILC}}=\widehat{A}{\left({\widehat{A}}^{\text{T}}{C}^{-1}\widehat{A}\right)}^{-1}{\widehat{A}}^{\text{T}}{C}^{-1}.$$(16)

We use the mixing matrix estimated in the previous step to compute the multidimensional (N_s-dimensional) ILC weights.

In regions of the map that have poor signal-to-noise, the GNILC algorithm estimates the signal subspace dimension to be zero (Fig. 7). This typically happens for regions at high Galactic latitude, for needlet bandpasses that cover the smaller scales, where the CO line emission is very faint and the noise is comparatively high. As a consequence, the traditional GNILC algorithm would assign zero weights for all and throw away all information in that needlet domain. As mentioned earlier, this results in two issues: (1) variation of map resolution from location to location on the sky, and (2) loss of faint signal. To avoid this drawback, we instead performed an inverse noise weighted combination by inverting the system d = āt + n, where ā is the N_ch-dimensional “mixing vector determined in the calibration step (Sect. 3.4), and t is a single CO template. This inverse noise weighted solution is obtained through the GLS solution: $${\text{W}}_{\text{GLS}}=\overline{a}{\left({\overline{a}}^{\text{T}}{C}_{\text{n}}^{-1}\overline{a}\right)}^{-1}{\overline{a}}^{\text{T}}{C}_{\text{n}}^{-1}.$$(17)

In summary, we can write the xGNILC weights as $$W_{\text{xGNLLC}}=\{\begin{array}{c}{W}_{\text{ILC}}{N}_{\text{s}}>0\\ W_{\text{GLS}}{N}_{\text{s}}=0\end{array}.$$(18)

3.5.4 Recombination into processed maps

We obtained re-processed maps from the GNILC cleaned needlet coefficient maps by inverse needlet transform, at HEALPix N_side of 1024 and at 10′ resolution. We finally replaced the regions of very bright CO that were masked/saturated as described in Sect. 3.3 with data in the same patch from the corresponding original input Planck CO data product.

3.6 Postprocessing

The postprocessing step is only relevant for the CO J = 3 → 2 data product from our pipeline. As discussed in Sect. 2.4, the CO J =3 → 2 has clearly visible systematic artifacts around the Galactic plane. These artifacts are not removed by the xGNILC pipeline⁶. In the postprocessing stage, we reduced contamination by those residual systematics by filtering. In practice, we filtered a difference map: $$\text{Δ}=d_{\text{xGNILC}}^{J=3\to 2}-\left({\overline{a}}_{J=3\to 2}/{\overline{a}}_{J=2\to 1}\right)d_{\text{xGNILC}}^{J=2\to 1}$$(19)

(i.e., a J = 3 → 2 map in which the CO signal is reduced by subtracting the J = 2 → 1 signal rescaled with the ratio of the global calibration factors). As the scaling relation is only approximate, we find several CO-associated features in the residual map due to the CO signal not canceling out. We threshold that difference map and masked all pixels where the absolute value of the residual was above 0.8 K_RJ km s⁻¹. We then median-filtered the masked residual map in superpixels at N_side = 16 and smooth that median-filtered map with a 90′ Gaussian beam. The resulting map captures the large scale features of the residuals shown in Fig. 4. We subtracted this median-filtered masked residual map from the xGNILC CO J = 3 → 2 map to produce our released data product for J = 3 → 2.

3.7 Uncertainties

The error in the CO maps comprises a statistical uncertainty from the measurement noise and a systematic uncertainty due to the choice of the calibration factor that is used as a prior. The systematic uncertainty for CO J = 1 → 0 and J =2 → 1 maps are small compared to the statistical uncertainty.

3.7.1 Statistical uncertainties

The input Planck TYPE 1 and TYPE 2 maps come with their individual uncertainty maps, but TYPE 1 and TYPE 2 maps are not independent, and they have significant cross-covariances. In addition, there is noise correlation between map pixels, and the full statistical description of the final statistical errors thus requires a set of simulated noise maps. We modeled the noise statistics in the input map across pixels and scales using the individual per-pixel uncertainty maps to generate pixel-dependent noise, which we then post-processed (with a multivariate filter in harmonic space) for the simulated noise to match noise auto-and cross-spectra computed from the null maps (see Appendix A for details). We provide 100 such generated simulations.

In Table 2, we show the statistical uncertainty that we computed directly from the xGNILC null maps for high Galactic latitude. We checked that the errors computed in this way for xGNILC match the standard deviation that we can compute from the noise simulation suite (to within a percentage level of accuracy or better).

3.7.2 Systematic uncertainties

We also obtained maps of the systematic uncertainties, which are shown in Fig. 11. The systematic uncertainty captures the impact of calibration uncertainties that enter the pipeline through the choice of a prior on the mixing vector ā (but not any calibration uncertainties already impacting the original Planck data products). For the xGNILC CO J =3 → 2 map, the calibration uncertainty also impacts the postprocessing step as we use the global scaling ratio to compute the residual map. We used the individual fit over 30 subregions (with offset), as discussed in Sect. 3.4, to construct a covariance matrix for the scaling parameter a_α. This covariance matrix with the mean values for a_α, given in Table 1, is used to define a multivariate Gaussian distribution for the mixing vector.

Systematic uncertainties are then obtained by running the xGNILC pipeline with 200 different multivariate Gaussian realizations of the prior used for ā. For the CO J =3 → 2 map, we performed the additional step of postprocessing the maps to correct for the residual systematics, with different values for the ratio of scaling factors. We computed the standard deviation of the 200 different output maps we obtained in the process. For the CO J =1 → 0 and J = 2 → 1 maps, the systematic uncertainty is only nonzero in the regions away from the Galactic plane where the xGNILC pipeline used the GLS weights. In the case of the CO J =3 → 2 maps, the systematic uncertainty of the postprocessing step is larger than the GLS step, and we see more systematic uncertainty in the Galactic plane. However, for all three maps, the systematic uncertainties are still negligible compared to the statistical uncertainty.

3.7.3 Contamination

Our basic xGNILC maps are obtained by combining TYPE 1 and TYPE 2 maps. This generates the results with the lowest noise contamination, but the possibility of contamination by other foregrounds originating from the TYPE 2 maps cannot be ruled out.

We hence repeated the whole process using only TYPE 1 maps (i.e., three input maps instead of five). We then computed the difference between xGNILC products obtained in the two cases and we flagged for possible significant contamination, for each of the CO lines, map pixels where all of the following three conditions are met. First, a discrepancy between the two data products (xGNILC with or without TYPE 2 maps) was detected locally at more than 4σ in the xGNILC difference map. Second, the difference was larger than 15% of the local xGNILC CO emission (i.e., we tolerate discrepancies that do not change the CO emission by more than 15%). Third, the CO signal itself in the final xGNILC CO map was detected locally at more than 4σ (i.e., we do not flag regions where there is very low signal that is barely detectable in our final maps). For each of the three CO lines, a confidence mask was generated by excluding a circle of radius 30′ around each pixel flagged in this way. We concatenated these three masks into a single confidence mask.

4 Data products

The xGNILC CO data products contain, for each CO line, the following set of maps:

• xGNILC CO map at N_side = 1024.

• xGNILC CO null map at N_side = 1024 obtained by projecting the Planck null maps with xGNILC weights.

• Systematic uncertainty maps for the xGNILC data product.

• 100 realizations of the projected noise for each xGNILC CO map⁷.

• all the masks used in the analysis and confidence mask for the data product.

• the beam transfer function of xGNILC data products. The xGNILC CO J =1 → 0 map and projected null map are shown in Fig. 8 along with corresponding figures showing the Planck TYPE 1 J =1 → 0 data product. Both are for an angular resolution of 10′.

A visual inspection of the maps (CO map and null) show a spectacular improvement in the signal-to-noise ratio of the CO J =1 → 0 line emission in the new xGNILC data product. We note that the structure of the xGNILC null map is more complex, with a clear distinction between the behavior near the Galactic center, most regions in the Galactic plane, and outside the Galactic CO emission region. The noise behavior is largely driven by the number of foreground subspace dimensions estimated by GNILC in the last needlet band, except for the masked bright sources such as the Galactic center, for which noise properties are identical to those of the corresponding TYPE 1 input map. The immediate region near the Galactic center is more noisy than the rest of the sky as GNILC tends to find higher foreground subspace dimensions in regions near the Galactic center. From Table 2, we see the typical statistical uncertainty from the noise at high galactic latitudes is far superior than both TYPE 1 and TYPE 2 CO J =1 → 0 maps when comparing at the same resolution.

The xGNILC CO J = 2 → 1 map and projected null map are shown in Fig. 9, with the corresponding TYPE 1 map and null being reconvolved with the xGNILC beam transfer function. The statistical uncertainty for the xGNILC J = 2 → 1 map shows a significant improvement over the TYPE 1 J = 2 → 1 map and has lower noise than the cleanest Planck CO data product - the TYPE 2 J = 2 → 1 map (see Table 2). The TYPE 1 CO J = 2 → 1 map has better signal-to-noise ratio than other Planck TYPE 1 maps. Again, the xGNILC CO J = 2 → 1 map has a visibly better signal-to-noise level than the original TYPE 1 CO J = 2 → 1 map. The discussion above for noise in the different regions of the sky for xGNILC CO J =1 → 0 applies for the CO J = 2 → 1 map as well.

In Fig. 10, we compare the xGNILC CO J = 3 → 2 map and projected null map with those of the Planck TYPE 1 CO J = 3 → 2 map, smoothed for a comparison at 10′ resolution of the xGNILC map. The xGNILC J = 3 → 2 map is visibly less noisy, and the systematics around the Galactic ridge have been removed. The xGNILC noise is significantly lower than the uncertainty of the corresponding TYPE 1 map as shown in Table 2.

In Fig. 11, we show the estimate of systematic uncertainty due to the uncertainty in the choice of priors on the mixing vector. For the xGNILC CO J =1 → 0 and J = 2 → 1 data products, the systematic uncertainty is only nonzero in the region outside the Galactic plane where we carry out a GLS. For the CO J = 3 → 2 map, there is more contamination in the Galactic plane because of the postprocessing stage, which also assumes this prior. The systematic uncertainty is subdominant for all three xGNILC data products.

In the preprocessing step of the pipeline, we masked most of the thermal SZ clusters present in the TYPE 2 maps (which would contaminate the xGNILC data products otherwise). However, all input maps also have point source contaminations that are not removed and propagate to xGNILC data products.

Additionally, there are few regions of bright extra-galactic CO emission from our local neighborhood, such as the Large Magel-lenic Cloud, Small Magellenic Cloud, M81, and so on. Given the different use cases for the data products, we leave it to the end user to mitigate these contaminations as per requirement.

Fig. 8 Comparison of xGNILC CO J =1 → 0 data products with Planck TYPE 1. Top: xGNILC CO J =1 → 0 map at output resolution of 10′ and the corresponding jackknife noise map. Bottom: similar maps for Planck TYPE 1 CO J = 1 → 0 data product and jackknife noise of Planck TYPE 1 CO J = 1 → 0 map, both smoothed to 10′ resolution.

Fig. 9 Comparison of xGNILC CO J = 2 → 1 data products with Planck TYPE 1. xGNILC CO J = 2 → 1 map at output resolution of 10′ and jackknife noise map for xGNILC CO J = 2 → 1 map are shown in the top row. In the bottom row, the Planck TYPE 1 CO J = 2 → 1 map and the jackknife noise of the Planck TYPE 1 CO J = 2 → 1 map are both smoothed to 10′ resolution.

5 Validation and discussion

For validation, we compared the xGNILC maps with the existing Planck CO data products. Planck product validation was extensively discussed in Planck Collaboration XIII (2014), and we make no attempt to revalidate the Planck CO maps as compared to other external data sets.

5.1 Map-level validation

In Planck Collaboration XIII (2014), three specific regions of CO emission were studied for validation: Taurus (l = 173°, b = −15°), Orion (l = 210°, b = −14°), and Polaris (l = 124°, b = 26°). In Figs. 12–14, we compare the xGNILC maps and Planck TYpE 1 maps for these three regions. We also show the difference map (Planck TYpE 1 – xGNILC) for these regions. For the CO J =1 → 0 map comparison, we also show the Dame CO J =1 → 0 map. The CO J =1 → 0 and J = 3 → 2 line emissions’ map comparisons are at a 20′ resolution to suppress the noise in the Planck TYpE 1 maps.

The visual comparison for CO J =1 → 0 maps shows excellent agreement between the Dame and xGNILC maps, with a slight excess in the xGNILC maps, as expected given the relative calibration coefficients, and consistent with a contribution from CO isotopologues. The Planck TYpE 1 map remains noisy even after smoothing to 20′, but it shows features consistent with the xGNILC map. Finally, the difference of the Planck TYpE 1 and xGNILC map shows no difference that is visually different from noise.

In Fig. 13, we present a similar comparison for the CO J = 2 → 1 maps. The comparison for this line, for which the signal-to-noise ratio is higher than for the other two, is done at the 10^' resolution of the xGNILC map. We see good agreement between the features seen in the xGNILC and the Planck TYPE 1 maps. The difference map for CO J = 2 → 1 shows some small ringing features that are seen in the xGNILC map in the Orion region, but these features are small in comparison to the brightness of the CO signal in this region. The difference maps in the other two regions are consistent with noise.

We show the CO J = 3 → 2 map comparison in Fig. 14 at 20′ resolution to reduce the noise in the Planck map. The signal-to-noise level of the Planck map shown here is poor, particularly in the Polaris region. In the Taurus and Orion regions, we find good agreement between the Planck TYPE 1 and xGNILC. The difference maps for Taurus and Polaris are consistent with noise. The difference map in the Orion region shows small variations between the Planck TYPE 1 and xGNILC, correlated with the signal, but this difference is small compared to the CO emission.

We also investigate high-latitude Galactic CO emission seen in the xGNILC data products. The Planck TYPE 1 CO maps are too noisy to contribute much reconstructed CO signal at these high latitudes. Most of the reconstruction relies on the TYPE 2 maps and can have residual dust contamination. To trace the dust in regions of interest, we use the IRAS 100 µm map (Miville-Deschênes & Lagache 2005).

We focus on the north ecliptic sky for this study to make use of the Dame & Thaddeus (2022) data. For this purpose, we took the DT22+DHT masked datacube⁸ and integrated the cube over the velocity axis. The velocity integrated data have a beam of ~8.4′ with a 0.25° sampling. We interpolated the data to a 1′ grid and projected it to Healpix pixellization with N_side of 1024. Finally, we changed the beam size from 8.4′ to 10′.

In Fig. 15, we show the comparison for three samples of high latitude regions. We see good agreement between the xGNILC CO map and the results of Dame & Thaddeus (2022). All three regions have strong correlation with the dust, but the xGNILC maps do not show any obvious excess dust contamination. In the northern ecliptic sky where we have the new Dame & Thaddeus (2022) observations of faint CO clouds (excluding the Galactic region of the previous surveys), the xGNILC CO J = 1 → 0 map has a Pearson’s correlation coefficient of 0.91 for pixels brighter than 0.1 K_RJ km s⁻¹. This indicates excellent agreement between the two observations and gives us confidence that the xGNILC CO maps do not suffer from very significant dust contamination. In Fig. 15, we also show the Planck TYPE 1 CO J = 1 → 0 maps in the three regions. The noise levels of Planck TYPE 1 maps make it difficult to see anything but the brightest CO clouds. Detecting the faint CO clouds at high Galactic latitude is significantly easier with the new xGNILC maps.

In Fig. 16, we plot a scatter of xGNILC map pixels with corresponding Planck TYPE 1 map pixels at 10′ resolution, downgraded to HEALPix N_side = 256. Unlike similar plots for Planck TYPE 2 versus TYPE 1 maps discussed in Sect. 2.4, we see no disagreement between the xGNILC and Planck TYPE 1 maps, and near-perfect correlation. At low pixel values, we find slightly increased spread due to the noise in the Planck TYPE 1 maps.

Next, we take a look at the difference between the original Planck TYPE 1 and the xGNILC CO maps. We smooth both maps to 30′ resolution, and then produce difference maps, shown in Fig. 17. For all three CO line emissions, we see small differences in the Galactic plane. In the case of CO J =1 → 0, this difference is barely above the noise threshold, and for both CO J =1 → 0 and J = 2 → 1 emissions, the difference in the Galactic plane is small compared to the signal in the Galactic plane. In the case of CO J = 3 → 2 maps, the difference maps is largely a display of the systematic residuals that were visible in the Planck TYPE 1 CO J =3 → 2 map, and they were subtracted by the postprocessing step described in Sect. 3.6.

Fig. 10 Comparison of xGNILC CO J = 3 → 2 data products with Planck TYPE 1. xGNILC CO J = 3 → 2 map at output resolution of 10′ and jackknife noise for xGNILC CO J = 3 → 2 map are shown in the top row. In the bottom row, the Planck TYPE 1 CO J = 3 → 2 map and the corresponding jackknife noise map are smoothed to 10′ resolution.

Fig. 11 Maps of systematic uncertainty in xGNILC data products arising from uncertainty in priors used for CO scaling parameters. The uncertainty in the CO J =1 → 0 and J = 2 → 1 maps are only in the region where we do a GLS, and so it does not contribute to any uncertainty in the Galactic plane. For the CO J = 3 → 2 map, the dominant systematics uncertainty comes from the postprocessing stage.

Fig. 12 Images (12.5° × 12.5°) of important CO emission regions for CO J = 1 → 0: Taurus (left column), Orion (middle column), and Polaris (right column), at 20′ resolution. The first row shows the Dame et al. (2001) map (unobserved pixels are shown in gray), second row shows our xGNILC map, third row shows the Planck TYPE 1 map, and in the last row we show the difference map (Planck TYPE 1 - xGNILC). We smoothed these maps to 20′ to suppress the noise in the Planck TYPE 1 map. We highlight the strong consistency between xGNILC maps and Dame maps.

Fig. 13 Images (12.5° × 12.5°) of important CO emission regions for CO J = 2 → 1: Taurus (left column), Orion (middle column), and Polaris (right column), at 10′ resolution. The first row shows our xGNILC map, the second row shows the Planck TYPE 1 map, and the difference map (Planck TYPE 1 – xGNILC) is shown in the third row.

5.2 Power spectrum validation

We now compare map power spectra (Figs. 18 and 19). For comparison of xGNILC with Planck TYPE 1, we smoothed the TYPE 1 maps to 10′, while for comparison with the TYPE 2 maps we smoothed the xGNILC maps to 15′ . We used HEALpix anafast to compute the power spectra on the full sky, as well as the brightest regions masked in the Galactic ridge using the mask shown in Fig. 3. The spectrum was smoothed with a bin width of seven multipoles to average out the fluctuations. The null maps were used to estimate the noise power spectra, and we show the TYPE 1 and TYPE 2 power spectra after noise debias-ing. There is no need to debias the power spectra for the xGNILC maps, as those are strongly signal-dominated over the full range of ℓ multipoles.

Comparisons of xGNILC and Planck TYPE 1 power spectra are shown in Fig. 18. In the xGNILC maps, the signal dominates over the noise for the full 2 ≤ ℓ ≤ 1000 multipole range. An improvement by orders of magnitude of the noise level is achieved with the xGNILC postprocessing. Signal spectra for the xGNILC and Planck TYPE 1 maps agree well where the TYPE 1 spectra are above the noise, but deviate from one another at high ℓ, where noise dominates. The power in xGNILC maps is lower than the power in the original maps after noise-debiasing. This is not unexpected, as even for the xGNILC algorithm (as compared to standard GNILC) some of the signal at small scales, where the input TYPE 1 maps have low signal-to-noise levels, is filtered in the step of projection onto the signal subspace. It is plausible that, in addition, null maps for the TYPE 1 products do not capture all of the “noise” properties (and in particular contamination by point sources), and thus there is a small amount of correlated error in the half-ring maps, which contributes to the TYPE 1 “signal” spectra, but it is at least partially filtered by xGNILC. For the CO J = 3 → 2 line emission specifically, systematic residuals present in the TYPE 1 map have been removed from the xGNILC map, also reducing the total power on a large scale.

We next compare the power spectra of the xGNILC maps with the Planck TYPE 2 maps (Fig. 19). Differences seen in Fig. 2 between TYPE 1 and TYPE 2 spectra are also present here, confirming that the signal in the xGNILC maps is strongly constrained to match the TYPE 1 data products. The xGNILC CO J =1 → 0 map has lower power than the corresponding TYPE 2 map with and without masking, and the full sky power spectrum for CO J = 2 → 1 xGNILC map is higher than the TYPE 2 map spectrum. However, the power spectra for the two J = 2 → 1 maps with masking agree almost perfectly. This is expected because for the J = 2 → 1 line calibration, coefficients for the TYPE 1 and TYPE 2 maps are almost the same, and the TYPE 2 CO J = 2 → 1 map has the best signal-to-noise ratio. Hence, the GLS weights in the noisiest regions would basically prescribe dominant contribution from the TYPE 2 CO J = 2 → 1 map.

In Fig. 19, we also note that the xGNILC CO J =1 → 0 maps have lower noise than the corresponding TYPE 2 maps. This demonstrates that (as expected from the reprocessing pipeline by design) the noise level in the xGNILC maps is lower than in any of the original Planck CO maps, both TYPE 1 and TYPE 2.

Fig. 14 Images (12.5° × 12.5°) of important CO emission regions for CO J = 3 → 2: Taurus (left column), Orion (middle column), and Polaris (right column), at 20′ resolution. The first row shows our xGNILC map, the second row shows the Planck TYPE 1 map, and the third row shows the difference map (Planck TYPE 1 – xGNILC). We have additional smoothing for this comparison as the signal-to-noise level of the Planck map is low.

Fig. 15 Images (12.5° × 12.5°) of some high-latitude CO clouds and dust. The first row shows the Dame & Thaddeus (2022) CO map, the second row our xGNILC CO J =1 → 0 map, the third row the Planck TYPE 1 CO J =1 → 0 map smoothed with the 10′ beam, and in the last row we show the IRIS 100 µm map also smoothed to 10′ beam.

Fig. 16 Scatter plot with pixel values of Planck TYPE 1 map as abscissa and pixel values of xGNILC map as ordinate. Top: CO J = 1 → 0; middle: J = 2 → 1; bottom: J = 3 → 2. Planck TYPE 1 maps are smoothed to 10′ resolution, and the HEALPix pixel size is increased to Nside = 256 for this figure. In the inset windows in the upper left corner we show the scatter for low-brightness pixels. In the inset windows located in the bottom right corner, we show the scatter plot in log scale. The orange line indicates y = x in all plots.

Fig. 17 Difference maps (Planck TYPE 1 – xGNILC) for CO J = 1 → 0 (top), CO J = 2 → 1 (middle), and CO J = 3 → 2 (bottom) are shown in this figure. Differences are visible in the Galactic plane, but they are a small fraction of the total signal and lower than calibration uncertainties. For CO J = 3 → 2, the effect of the postprocessing of residual systematics shows up as a negative trough close to the Galactic plane and extended bright red regions above and below the Galactic ridge.

5.3 Limitations of the new maps

In the previous two sections we validate our xGNILC data products by comparing them to the Planck CO data products. It is reasonable to argue that our maps are very high signal-to-noise CO maps that have most of the good features of the Planck TYPE 1 CO data products. We reduce the contamination by noise and improve the signal-to-noise ratio on the small CO signal buried below the noise in the original TYPE 1 maps. In addition, wherever there are discrepancies between TYPE 1 and TYPE 2 maps in high signal-to-noise ratio regions, the xGNILC signal remains close to the original TYPE 1 signal, and hence it is presumably less contaminated by other Galactic foregrounds. xGNILC maps have also been cleaned from SZ clusters visible in TYPE 2 maps, as well as from residual systematic effects in the case of the J = 3 → 2 line.

In spite of these improvements, some limitations remain. Original calibration errors in the TYPE 1 maps, as well as contributions from CO isotopologues, are basically unchanged. The J = 3 → 2 map, originally very noisy, is much cleaner after reprocessing, but this improvement is achieved in combination with maps of the other two lines, and some of the signal-correlated differences between the xGNILC and the TYPE 1 maps, visible, for instance, in the Orion region, may be due to some level of line mixing.

Our final data products combine TYPE 1 and TYPE 2 maps. One should note that strong correlation between the noise in the three line maps exist after the xGNILC reprocessing.

Fig. 18 Power spectrum comparison for Planck TYPE 1 CO maps (blue) and xGNILC CO maps (orange) at 10′ resolution. We noise debiased the TYPE 1 map spectra shown in the figures. The top figure shows spectra comparison for the CO J =1 → 0 maps, the middle figure shows them for the CO J = 2 → 1 maps, and the bottom figure shows them for the CO J = 3 → 2 maps.

Fig. 19 Power spectrum comparison for Planck TYPE 2 CO maps (green) and xGNILC CO maps (orange) at 15′ resolution. The top figure shows a spectra comparison for the CO J =1 → 0 maps, and the bottom figure shows them for the CO J = 2 → 1 maps.

6 Conclusion

In this work, existing Planck space mission maps of full sky rotational lines of Galactic CO emission have been revisited and post-processed to enhance signal-to-noise ratio and reduce contamination by systematics and astrophysical confusion. Maps are produced at a common (apodized) 10′ angular resolution. The data products include, in addition to the three CO line full-sky maps, the pipeline’s propagated null maps and maps of systematic uncertainty, as well as processing and confidence masks. All data products are made available to the scientific community on a dedicated website⁹.

Acknowledgements

The authors acknowledge useful discussions with Francois Boulanger, Jean-Loup Puget, and Reijo Keskitalo. MR acknowledges support from the CSIC program Ayuda a la Incorporación de Científicos Titulares provided under the project 202250I159. Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) python package (Zonca et al. 2019).

Appendix A Statistical uncertainties

The noise in the xGNILC data products depends on the sky pixel, is correlated between the different data products, and is correlated between pixels. For this reason, simple standard deviation maps do not fully capture the properties of the statistical uncertainty of the xGNILC data products.

Fig. A.1 Figure shows xGNILC null map (left) and 0th noise realization of the xGNILC noise simulation suite (right) for CO J=1 →0. The color scale spans −0.5 to 0.5 KRJ km s−1.

Fig. A.2 Figure shows xGNILC null map (top) and 0th noise realization of the xGNILC noise simulation suite (bottom) for CO J = 2→ 1. The color scale spans −0.25 to 0.25 KRJ km s−1.

Therefore, we provide the end user with 100 simulated noise realizations for each xGNILC CO line emission data product, which are generated as follows.

1. For each of the input Planck TYPE 1 and TYPE 2 inputs, we generated maps of white noise, which were scaled in each pixel using the standard deviation map for the corresponding CO data product.

2. We smoothed or deconvolved those five maps to put them at a common 10′ angular resolution.

3. We computed the multivariate power spectrum C_w of these maps (for each ℓ, a 5 × 5 matrix, with vanishing off-diagonal terms).

4. Similarly, we computed the multivariate power spectrum C_n of the null maps. We set cross-spectra terms between TYPE 1 maps, which were produced using independent data sets, to zero, but we kept nonzero values for cross-spectra between TYPE 1 and TYPE 2 maps, as well as between the maps for the two TYPE 2 lines, since the generation of each of the two TYPE 2 maps makes use of Planck observations also used to generate the other CO products.

5. We obtained simulated noise by multiplying the a_ℓms obtained in step 2 by $C_n^{1/2}{C}_w^{-1/2}$. In practice, the filter is smoothed using a box-cart average.

Fig. A.3 Figure shows xGNILC null map (top) and 0th noise realization of the xGNILC noise simulation suite (bottom) for CO J = 3 → 2. The color scale spans −0.25 to 0.25 KRJ km s−1.

Fig. A.4 Power spectrum comparison for xGNILC null maps (dotted) and mean of 100 realizations of simulated xGNILC noise (solid). For the cross-spectra, we plot the absolute values. We compute the spectra with the brightest CO mask shown in Fig. 3. The simulated noise agrees well with the actual xGNILC null maps.

The simulated noise maps, shown in Figs. A.1–A.3, appear consistent with xGNILC null maps. Our method of simulating the noise works well for 98% of the sky with a sub-percent-level difference in the noise level. On the brightest 2% of the sky, corresponding to the mask shown in Fig. 3, differences of 8%, 4%, and 12% for the CO J = 1 → 0, J = 2 → 1, and J = 3 → 2 maps, respectively, are probably due to residual systematic errors that are not fully captured by the noise simulation process. In Fig. A.4, we compare the noise auto and cross-power spectra of those simulations to the auto and cross-spectra of xGNILC null maps, with 2% sky masking (mostly toward the Galactic ridge), and we find excellent agreement.

References

All Tables

All Figures

Fig. 1 Scatter plot with pixel values of TYPE 1 CO map as abscissa and pixel values of TYPE 2 CO maps as ordinate. We smooth all maps to 30 arcmin resolution and increase the HEALPix pixel size to Nside = 256 for this figure. The top figure shows the scatter for CO J = 1 → 0 Planck maps, and the bottom figure shows it for CO J = 2 → 1 maps. For large signals, the signal in TYPE 2 maps is higher than in TYPE 1 maps for the J = 1 → 0 line, and lower for the J = 2 → 1 line. In the inset window, we show that the low-brightness pixels do not show obvious inconsistency. This is in agreement with the comparison made in Planck Collaboration XIII (2014).

In the text

Fig. 2 Power spectra comparison of TYPE 1 and TYPE 2 CO J =1 → 0 maps (top) and J = 2 → 1 maps (bottom), with and without an attenuation mask to reduce the amplitude of emission in the brightest 2% of pixels (shown in the top panel of Fig. 3). The TYPE 1 J =1 → 0 have lower power compared to the TYPE 2 map. For the CO J = 2 → 1 line emission, the TYPE 1 map has higher power compared to the TYPE 2 map. Both of them have higher power than the Commander CO J = 2 → 1 map. The discrepancy between TYPE 1 and TYPE 2 CO J = 2 → 1 power spectra largely goes away when masking the regions of strongest CO.

In the text

Fig. 3 Masks used in this work. Top: mask used to compute power spectrum after excluding the brightest CO emission regions in the Galactic ridge. Bottom: mask used for fitting global CO mixing coefficients at low Galactic latitude the Galactic plane, obtained from the sky fraction observed in the Dame et al. (2001) survey, excluding regions with CO emission larger than 3 KRJ km s−1. This results in the brightest 0.5% of pixels (in the Galactic ridge) being masked.

In the text

	Fig. 4 TYPE 1 CO J = 3 → 2 data product from Planck, shown here after smoothing to 30′ resolution. We are showing a region ±90° in Galactic longitude and ±30° in Galactic latitude about the Galactic center. Ringing effects and systematic residuals are visible above and below the Galactic ridge.
In the text

	Fig. 5 Needlet bandpass windows, $h_{\ell }^j$, showing the selection weights in harmonic space (plotted in black). The “apodized” beam transfer function of the output map resolution of 10′ is also shown in red. The end of the beam transfer function is modified to smoothly go to zero at ℓmax = 2048.
In the text

Fig. 6 Maps of calibration obtained by fitting Planck TYPE 1 (top row, from J = 1 → 0 to J = 3 → 2, left to right) and TYPE 2 (bottom row, from J =1 → 0 to J = 2 → 1, left to right) maps against Dame J = 1 → 0 map in Nside = 16 super pixels. We perform a fit in any super pixel if there are 50 or more original pixels at Nside = 512 with a high signal-to-noise ratio in the Dame map. For all super-pixels where this condition is not met, we assume the mean calibration given in Table 1. These maps illustrate the variation over the sky of the ratio between the CO detected in Planck TYPE 1 and TYPE 2 maps and the CO detected in the Dame J = 1 → 0 map.

In the text

Fig. 7 Maps of dimension of CO subspace as determined by GNILC for the 60′ needlet band (top) and the 10′ needlet band (bottom). At smaller angular scales, larger parts of the high-latitude sky becomes noise-dominated (dark blue in the figure). The standard GNILC algorithm assigns those regions a zero-dimensional signal subspace, which results in varying effective resolution in standard GNILC map. In our processing, in those needlet domains, we produce the CO maps using the GLS weights. With this modification to the GNILC algorithm, final xGNILC maps are produced with a uniform angular resolution.

In the text

	Fig. 8 Comparison of xGNILC CO J =1 → 0 data products with Planck TYPE 1. Top: xGNILC CO J =1 → 0 map at output resolution of 10′ and the corresponding jackknife noise map. Bottom: similar maps for Planck TYPE 1 CO J = 1 → 0 data product and jackknife noise of Planck TYPE 1 CO J = 1 → 0 map, both smoothed to 10′ resolution.
In the text

	Fig. 9 Comparison of xGNILC CO J = 2 → 1 data products with Planck TYPE 1. xGNILC CO J = 2 → 1 map at output resolution of 10′ and jackknife noise map for xGNILC CO J = 2 → 1 map are shown in the top row. In the bottom row, the Planck TYPE 1 CO J = 2 → 1 map and the jackknife noise of the Planck TYPE 1 CO J = 2 → 1 map are both smoothed to 10′ resolution.
In the text

	Fig. 10 Comparison of xGNILC CO J = 3 → 2 data products with Planck TYPE 1. xGNILC CO J = 3 → 2 map at output resolution of 10′ and jackknife noise for xGNILC CO J = 3 → 2 map are shown in the top row. In the bottom row, the Planck TYPE 1 CO J = 3 → 2 map and the corresponding jackknife noise map are smoothed to 10′ resolution.
In the text

	Fig. 11 Maps of systematic uncertainty in xGNILC data products arising from uncertainty in priors used for CO scaling parameters. The uncertainty in the CO J =1 → 0 and J = 2 → 1 maps are only in the region where we do a GLS, and so it does not contribute to any uncertainty in the Galactic plane. For the CO J = 3 → 2 map, the dominant systematics uncertainty comes from the postprocessing stage.
In the text

Fig. 12 Images (12.5° × 12.5°) of important CO emission regions for CO J = 1 → 0: Taurus (left column), Orion (middle column), and Polaris (right column), at 20′ resolution. The first row shows the Dame et al. (2001) map (unobserved pixels are shown in gray), second row shows our xGNILC map, third row shows the Planck TYPE 1 map, and in the last row we show the difference map (Planck TYPE 1 - xGNILC). We smoothed these maps to 20′ to suppress the noise in the Planck TYPE 1 map. We highlight the strong consistency between xGNILC maps and Dame maps.

In the text

	Fig. 13 Images (12.5° × 12.5°) of important CO emission regions for CO J = 2 → 1: Taurus (left column), Orion (middle column), and Polaris (right column), at 10′ resolution. The first row shows our xGNILC map, the second row shows the Planck TYPE 1 map, and the difference map (Planck TYPE 1 – xGNILC) is shown in the third row.
In the text

Fig. 14 Images (12.5° × 12.5°) of important CO emission regions for CO J = 3 → 2: Taurus (left column), Orion (middle column), and Polaris (right column), at 20′ resolution. The first row shows our xGNILC map, the second row shows the Planck TYPE 1 map, and the third row shows the difference map (Planck TYPE 1 – xGNILC). We have additional smoothing for this comparison as the signal-to-noise level of the Planck map is low.

In the text

	Fig. 15 Images (12.5° × 12.5°) of some high-latitude CO clouds and dust. The first row shows the Dame & Thaddeus (2022) CO map, the second row our xGNILC CO J =1 → 0 map, the third row the Planck TYPE 1 CO J =1 → 0 map smoothed with the 10′ beam, and in the last row we show the IRIS 100 µm map also smoothed to 10′ beam.
In the text

Fig. 16 Scatter plot with pixel values of Planck TYPE 1 map as abscissa and pixel values of xGNILC map as ordinate. Top: CO J = 1 → 0; middle: J = 2 → 1; bottom: J = 3 → 2. Planck TYPE 1 maps are smoothed to 10′ resolution, and the HEALPix pixel size is increased to Nside = 256 for this figure. In the inset windows in the upper left corner we show the scatter for low-brightness pixels. In the inset windows located in the bottom right corner, we show the scatter plot in log scale. The orange line indicates y = x in all plots.

In the text

Fig. 17 Difference maps (Planck TYPE 1 – xGNILC) for CO J = 1 → 0 (top), CO J = 2 → 1 (middle), and CO J = 3 → 2 (bottom) are shown in this figure. Differences are visible in the Galactic plane, but they are a small fraction of the total signal and lower than calibration uncertainties. For CO J = 3 → 2, the effect of the postprocessing of residual systematics shows up as a negative trough close to the Galactic plane and extended bright red regions above and below the Galactic ridge.

In the text

	Fig. 18 Power spectrum comparison for Planck TYPE 1 CO maps (blue) and xGNILC CO maps (orange) at 10′ resolution. We noise debiased the TYPE 1 map spectra shown in the figures. The top figure shows spectra comparison for the CO J =1 → 0 maps, the middle figure shows them for the CO J = 2 → 1 maps, and the bottom figure shows them for the CO J = 3 → 2 maps.
In the text

	Fig. 19 Power spectrum comparison for Planck TYPE 2 CO maps (green) and xGNILC CO maps (orange) at 15′ resolution. The top figure shows a spectra comparison for the CO J =1 → 0 maps, and the bottom figure shows them for the CO J = 2 → 1 maps.
In the text

	Fig. A.1 Figure shows xGNILC null map (left) and 0th noise realization of the xGNILC noise simulation suite (right) for CO J=1 →0. The color scale spans −0.5 to 0.5 KRJ km s−1.
In the text

	Fig. A.2 Figure shows xGNILC null map (top) and 0th noise realization of the xGNILC noise simulation suite (bottom) for CO J = 2→ 1. The color scale spans −0.25 to 0.25 KRJ km s−1.
In the text

	Fig. A.3 Figure shows xGNILC null map (top) and 0th noise realization of the xGNILC noise simulation suite (bottom) for CO J = 3 → 2. The color scale spans −0.25 to 0.25 KRJ km s−1.
In the text

	Fig. A.4 Power spectrum comparison for xGNILC null maps (dotted) and mean of 100 realizations of simulated xGNILC noise (solid). For the cross-spectra, we plot the absolute values. We compute the spectra with the brightest CO mask shown in Fig. 3. The simulated noise agrees well with the actual xGNILC null maps.
In the text