A 3D model for CO molecular line emission as a potential CMB polarization contaminant

TL;DR

This paper develops a 3D simulation model for CO molecular line emission in galaxies, calibrates it with Planck data, and forecasts its potential to contaminate CMB polarization measurements, especially B-modes.

Contribution

The paper introduces the MCMole3D Python package for simulating CO emission and assesses its impact on CMB polarization studies, a novel approach for contamination forecasting.

Findings

Calibrated the model using Planck data on the Galactic plane.

Forecasted CO contamination level equivalent to r ≈ 0.02 in CMB B-mode polarization.

Provided a publicly available simulation tool for future CMB analysis.

Abstract

We present a model for simulating Carbon Monoxide (CO) rotational line emission in molecular clouds, taking account of their 3D spatial distribution in galaxies with different geometrical properties. The model implemented is based on recent results in the literature and has been designed for performing Monte-Carlo simulations of this emission. We compare the simulations produced with this model and calibrate them, both on the map level and on the power spectrum level, using the second release of data from the Planck satellite for the Galactic plane, where the signal-to-noise ratio is highest. We use the calibrated model to extrapolate the CO power spectrum at low Galactic latitudes where no high sensitivity observations are available yet. We then forecast the level of unresolved polarized emission from CO molecular clouds which could contaminate the power spectrum of Cosmic Microwave Background (CMB) polarization B-modes away from the Galactic plane. Assuming realistic levels of the polarization fraction, we show that the level of contamination is equivalent to a cosmological signal with $r \lesssim 0.02$. The Monte-Carlo MOlecular Line Emission (MCMole3D) Python package, which implements this model, is being made publicly available.

Figures

Tables

Table 1. Table 1 : List of parameters used in MCMole3D  simulations. † only for LogSpiral .

MCMole3D Default parameters
$N_{clouds}$	40,000
$R_{ring}$ [kpc]	5.3
$L_{min}$ [pc]	0.3
$L_{max}$ [pc]	60
${\sigma }_{z,0}$ [pc]	100
$h_R$ [kpc]	9
$R_{bar}^{†}$ [kpc]	3
$i^{†}$ [deg]	-12
${ϵ}_c\left[K_{\mathrm{RJ}}\mathrm{km}{s}^{-1}\right]$	240
$R_{em}$ [kpc]	6.6
$L_0$ [pc]	[5,50] Default: 20
${\sigma }_{ring}$ [kpc]	[1,5] Default:2.5

Table 2. Table 2 : Summary table of best fit parameters obtained using the two different Planck CO maps.

	$L_0$ [pc]	${\sigma }_{ring}$ [kpc]	$A_{CO}$	${\tilde{\chi }}^2$	dof	$p$-value	${\rho }_{L\sigma }$
Type 1	$14.50\pm 0.58$	$2.76\pm 0.19$	$0.69\pm 0.06$	$1.48$	11	0.13	0.74
Type 2	$11.59\pm 1.09$	$3.11\pm 0.32$	$1.37\pm 0.19$	1.95	11	0.03	0.92

Table 3. Table 3 : Best-fit parameters for the Axisymmetric   MCMole3D  model.

	$L_0$ [pc]	${\sigma }_{ring}$ [kpc]	$A_{CO}$	${\tilde{\chi }}^2$	dof	$p$-value	${\rho }_{L\sigma }$
Type 1	$19.47\pm 12.68$	$2.12\pm 0.23$	$1.00\pm 0.12$	$7.35$	11	0.00	0.99
Type 2	$16.24\pm 17.56$	$2.12\pm 0.30$	$2.25\pm 0.35$	$18.08$	11	0.00	0.99

Equations

$$z_{1/2}(R)\propto {\sigma }_z(R)={\sigma }_{z,0}\cosh \left(\frac{R}{h_R}\right),$$

$$z(R)=\frac{1}{\sqrt{2\pi }{\sigma }_z(R)}\exp \left[-{\left(\frac{z}{\sqrt{2}{\sigma }_z(R)}\right)}^2\right].$$

$${ϵ}_0(R)={ϵ}_c\exp \left(R/R_{em}\right),$$

$$\xi(L)=\frac{dn}{dL}\propto\left\{\begin{array}[]{ll}L^{0.8}&\mbox{if }L_{min}<L<L_{0},\\ L^{-\alpha_{L}}&\mbox{if }L_{0}<L<L_{max},\end{array}\right.$$

$$\mathcal{P}(<L)={\int }_{L_{min}}^{L}d{L}^{\prime }\xi (L^{\prime }).$$

$${\Theta }^{CO}=\{$$

$${\sigma }_{z,0},h_R,L_{min},L_{max},L_0\}.$$

$$\varphi (R)=A\log R+B,$$

$$\Delta {\tilde{C}}_{\ell }=\sqrt{\frac{2}{\nu }}(C_{\ell }+N_{\ell })$$

$$I^X(l)=\int db{I}^X(l,b),$$

$$I_{tot}^X=\int dldb{I}^X(l,b),$$

$$f=\frac{I_{tot}^{observ}}{I_{tot}^{model}}.$$

$$L_0=14.50\pm 0.58\mathrm{pc},$$

$${\sigma }_{ring}=2.76\pm 0.19\mathrm{Kpc},$$

$$A_{CO}=0.69\pm \mathrm{0.06.}$$

$${\mathcal{C}}_{\ell }^{high,CO}={\mathcal{C}}_{\ell }^{Gal}\frac{\mathrm{var}(high)}{\mathrm{var}(Gal)}.$$

$$Q(\hat{n}{)}^{CO}$$

$$U(\hat{n}{)}^{CO}$$

Full text

††pubyear: 2016††pagerange: A 3D model for CO molecular line emission as a potential CMB polarization contaminant–A 3D model for CO molecular line emission as a potential CMB polarization contaminant

1 Introduction

The Carbon monoxide (CO) molecule is one of the most interesting molecules present in molecular clouds within our Galaxy. Although the most abundant molecule in Galactic molecular clouds is molecular hydrogen ($\mathrm{H_{2}}$), it is inconvenient to use the emission from that as a tracer since it is difficult to detect because of having a low dipole moment and so being a very inefficient radiator. We therefore need to resort to alternative techniques for tracing molecular clouds using rotational or vibrational transitions of other molecules such as CO. Observations of CO emission are commonly used to infer the mass of molecular gas in the Milky Way by assuming a linear proportionality between the CO and $H_{2}$ densities via the CO-to-$\mathrm{H}_{2}$ conversion factor, $X_{CO}$. A commonly accepted value for $X_{CO}$ is $2\times 10^{20}\,\mathrm{molecules\cdot cm^{-2}\,\left(K\,km\,s^{-1}\right)^{-1}}$, although this could vary with position in the Galactic plane, particularly in the outer Galaxy (Balser et al., 2011).

The most intense CO rotational transition lines are the $J=1\rightarrow 0,\,2\rightarrow 1,\,3\rightarrow 2$ transitions at sub-millimetre wavelengths (115, 230 and 345 GHz respectively). These can usually be observed in optically thick and thermalized regions of the interstellar medium. Traditionally, the observations of standard $\mathrm{{}^{12}CO}$ emission are complemented by measurements of $\mathrm{{}^{13}CO}$ lines. Being less abundant (few percent), this isotopologue can be exploited for inferring the dust extinction in nearby clouds and hence providing a better constraint for measuring the $\mathrm{H_{2}}$ abundance (Bally et al., 1987; Jackson et al., 2006). However, there is growing evidence that $\mathrm{{}^{13}CO}$ regions could be associated with colder and denser environments, whereas $\mathrm{{}^{12}CO}$ emission originates from a diffuse component of molecular gas (Roman-Duval et al., 2016).

The spatial distribution of the CO line emission reaches a peak in the inner Galaxy and is mostly concentrated close to or within the spiral arms, in a well-defined ring, the so-called molecular ring between about $4-7$ kpc from the Galactic centre. This property is not unique to the Milky Way but is quite common in barred spiral galaxies (see Regan et al. (2002) for further references). The emission in the direction orthogonal to the Galactic plane is confined within a Gaussian slab with roughly 90 pc full width half maximum (FWHM) in the inner Galaxy getting broader towards the outer Galactic regions, reaching a FWHM of several hundred parsecs outside the solar circle. In the centre of the Galaxy, we can also identify a very dense CO emission zone ,rich in neutral gas and individual stars, stretching out to about 700 light years (ly) from the centre and known as the Central Molecular Zone.

Since the 1970s, many CO surveys of the Galactic plane have been carried out with ground-based telescopes, leading to accurate catalogues of molecular clouds (Dame et al., 2001; Mizuno and Fukui, 2004). Usually these surveys have observed a strip of $|b|\lesssim 5\deg$ around the Galactic plane. At higher Galactic latitudes ($|b|>30\deg$), the low opacity regions of both gas and dust, together with a relatively low stellar background which is useful for spotting extinction regions, complicate the observation of CO lines making this very challenging. In fact, only $\approx 100$ clouds have been detected so far in these regions.

The Planck satellite team recently released CO emission maps of the lowest rotational lines, $J=1-0,\,2-1\,,3-2$ observed in the 100, 217, 353 GHz frequency channels of the High Frequency Instrument (HFI) (Planck Collaboration et al., 2014a, 2016a). These were sensitive enough to map the CO emission even though the widths of these lines are orders of magnitude narrower than the bandwidth of the Planck frequency channels. These single frequency maps have been processed with a dedicated foreground cleaning procedure so as to isolate this emission. The Planck maps were found to be broadly consistent with the data from other CO surveys (Dame et al., 2001; Heyer and Dame, 2015), although they might still be affected by residual astrophysical emissions and instrumental systematics. In Figure 1, we show the so called Type 1 Planck map of the CO $J:1-0$ line Planck Collaboration et al. (2014a)111http://pla.esac.esa.int/pla which will be used in the following.

Many current and future CMB polarization experiments222For a complete list of the operating and planned probes see e.g. lambda.gfsc.nasa.gov are designed to exploit the faint B-mode signal of CMB polarization as a cosmological probe, in particular to constrain the physics of large scale structure formation or the inflationary mechanism in the early universe (Seljak and Zaldarriaga, 1997; Hu and White, 1997). One of the main challenges in the way of achieving these goals is the contamination of the primordial CMB signal by diffuse Galactic emission. In this respect, the synchrotron and thermal dust emission are known to be potentially the most dangerous contaminants, because they are intrinsically polarized. In fact, several analyses conducted on Planck and Wilkinson Microwave Anisotropy Probe (WMAP) data from intermediate and high Galactic latitudes at high (Planck Collaboration et al., 2016b) and low frequencies (Krachmalnicoff et al., 2015; Planck Collaboration et al., 2016c) showed that these emissions are dangerous at all microwave frequencies and locations on the sky (even if far from the galactic plane), confirming early studies using the WMAP satellite (Gold et al., 2011; Page et al., 2007; Baccigalupi, 2003).

Appropriate observations and theoretical investigations and modelling of polarized foreground emission for all emissions at sub-mm frequencies are therefore crucial for the success of future experiments. As these will observe at frequencies overlapping with the CO lines, unresolved CO line emission could significantly contaminate these measurements as well.

CO lines are in fact expected to be polarized at the percent level or below (Goldreich and Kylafis, 1981) because of interaction of the magnetic moment of the molecule with the Galactic magnetic field. This causes the so-called Zeeman splitting of the rotational quantum levels $J$ into the magnetic sub-levels $M$ which are intrinsically polarized. Moreover, if molecular clouds are somehow anisotropic (e.g when in the presence of expanding or collapsing envelopes in star formation regions) or are asymmetric, population imbalances of the $M$ levels can arise. This leads to different line intensities depending on the directions (parallel or perpendicular to the magnetic field) and to a net linearly polarized emission. Greaves et al. (1999) detected polarization in five star-forming regions near to the Galactic Centre while observing the CO lines $J=2-1,\,3-2$ and the $J=2-1$ of the isotopologue ${}^{13}\mathrm{CO}$. The degree of polarization ranged from $0.5\mathrm{\,to\,}2.5\%$. Moreover, the deduced magnetic field direction was found to be consistent with previous measurements coming from dust polarimetry, showing that the polarized CO emission could become a sensitive tracer of small-scale Galactic magnetic fields.

The goal of this paper is to propose a statistical 3D parametric model of CO molecular cloud emission, in order to forecast the contamination of CMB signal by this, including in the polarization. Being able to perform statistical simulation of this emission is crucial for assessing the impact of foreground residual uncertainties on cosmological constraints coming from the CMB. In addition, the capability of modeling the Galactic foreground emission in its full complexity taking into account line-of-sight effects is becoming necessary in light of the latest experimental results and the expected level of sensitivity for the future experiments (Tassis and Pavlidou, 2015; Planck Collaboration et al., 2017). In section 2 we present the assumptions made for building the model and the simulation pipeline for its implementation. In section 3 we describe the methodology for calibrating the CO simulations to match Planck observations. In subsection 3.4 we show how the parameters describing molecular cloud distribution shape the angular power spectrum of CO emission. Finally, in section 4 we forecast the expected level of polarized CO contaminations for the CMB  B-modes at high Galactic latitudes using our calibrated simulation of section 3 to infer statistically the emission at high Galactic latitude, where current observations are less reliable.

2 Building a statistical 3D CO emission model

In order to build an accurate description of CO emission in the Galaxy, we collected the most up to date astrophysical data present in the literature concerning the distribution of molecular gas as a function of the Galactic radius ($R$) and the vertical scale of the Galactic disk ($z$) as well as of the molecular size and the mass function. The model has been implemented in a Python package named MCMole3D333https://github.com/giuspugl/MCMole3D which is being made publicly available, and we present details of it in this Section444In the following we will refer to this model as the MCMole3D model for the sake of clarity.. The model builds on and extends the method proposed by Ellsworth-Bowers et al. (2015) who conducted a series of analyses distributing statistically a relative large number of molecular cloud objects according to the axisymmetric distribution of $\mathrm{H}_{2}$ observed in the Galaxy (Wolfire et al., 2003).

2.1 CO cloud spatial distribution

As mentioned in the introduction, the CO emission is mostly concentrated around the molecular ring. We have considered and implemented two different spatial distributions of the molecular clouds: an axisymmetric ring-shaped one and one with 4 spiral arms, as shown in Figure 3(b) and (a) respectively. The first is a simplified model and is parametrized by $R_{ring}$,and $\sigma_{ring}$ which are the radius and the width of the molecular ring respectively. On the other hand, the spiral arm distribution is in principle closer to the symmetry of our Galaxy and is therefore more directly related to observations. The distribution is described by two more parameters than for the axisymmetric case: the arm width and the spiral arm pitch angle. For the analysis conducted in the following sections, we fixed the value of the pitch angle to be $i\sim-13\deg$ following the latest measurements of Davis et al. (2012); Bobylev and Bajkova (2013) and fixed the arm half-width to be $340$ pc (Vallée, 2014).

Bronfman et al. (1988) found that the vertical profile of the CO emissivity can be optimally described by a Gaussian function of $z$ centred on $z_{0}$ and having a half-width $z_{1/2}$ from the Galactic plane at $z=0$. Both of the parameters $z_{0}$ and $z_{1/2}$ are in general functions of the Galactic radius $R$ (see Heyer and Dame (2015) for recent measurements). Since we are interested in the overall distribution of molecular clouds mainly in regions close to the Galactic plane, where data are more reliable, we adopted this parametrization but neglected the effects of the mid-plane displacement $z_{0}$ and set it to a constant value $z_{0}=0$, following Delabrouille et al. (2013). The vertical profile is then parametrized just by $z_{1/2}$ and mimics the increase of the vertical thickness scatter that is observed when moving from the inner Galaxy towards the outer regions:

[TABLE]

where $\sigma_{z,0}=0.1$ and $h_{R}=9$ kpc corresponds to the radius where the vertical thickness starts increasing. The half-width $z_{1/2}$ is related to $\sigma_{z}$ through the usual relation $z_{1/2}=\sqrt{2\ln 2}\sigma_{z}$. The final vertical profile is then:

[TABLE]

2.2 CO cloud emission

The key ingredients for modeling the molecular cloud emission are the dimension of the cloud and its typical emissivity. We assume an exponential CO emissivity profile which is a function of the Galactic radius following Heyer and Dame (2015); Roman-Duval et al. (2016):

[TABLE]

where $\epsilon_{c}$ is the typical emissivity of a particular CO line observed towards the centre of the Galaxy and $R_{em}$ the scale length over which the emissivity profile changes. Clouds observed in the outer Galaxy are in fact dimmer.

We then assume the distribution of cloud size $\xi(L)$ defined by their typical size scale, $L_{0}$, the range of sizes $\left[L_{min},L_{max}\right]$ and two power-laws with spectral indices (Roman-Duval et al., 2010)

[TABLE]

with $\alpha_{L}=3.3,3.9$ for clouds inside or outside the solar circle respectively. From the cloud size function $\xi(L)$ we derive the corresponding probability $\mathcal{P}(L)$ of having clouds with sizes smaller than $L$:

[TABLE]

The probability functions for different choices of the spectral index $\alpha_{L}$ are shown in Figure 2. We then inverted Equation 5 to get the cloud size associated with a given probability $L(p)$. The cloud sizes are drawn from a uniform distribution in $\left[0,1\right]$. The histograms of the sizes generated following this probability function are shown in the top panel of Figure 2 and are peaked around the most typical size $L_{0}$. In the analysis presented in the following $L_{0}$ is considered as a free parameter.

Finally, we assume a spherical shape for each of the simulated molecular clouds once they are projected on the sky. However, we implemented different emissivity profiles that are function of the distance from the cloud center, such as Gaussian or cosine profiles. These are particularly useful because, by construction, they give zero emissivity at the boundaries555For the Gaussian profile, we set $\sigma$ in order to have the cloud boundaries at $6\sigma$, i.e. where the Gaussian function is zero to numerical precision. and the maximum of the emissivity in the centre of the projected cloud on the sky. This not only mimics a decrease of the emission towards the outer regions of the cloud, where the density decreases, but also allows to minimize numerical artifacts when computing the angular power spectrum of the simulated maps (see section 3). An abrupt top-hat transition at the boundary of each cloud would in fact cause ringing effects that could bias the estimate of the power spectrum.

2.3 Simulation procedure

The model outlined in the previous Section enables statistical simulations of CO emission in our Galaxy to be performed for a given set of free parameters $\Theta^{CO}$ that can be set by the user:

[TABLE]

The values chosen for our analysis are listed in Table 1.

For each realization of the model, we distribute by default 40,000 clouds within our Galaxy. This number is adopted for consistency with observations when observational cuts are applied (for further details see Ellsworth-Bowers et al. (2015)). The product of each simulation is a map, similar to the one in Figure 4, in the Hierarchical Equal Area Latitute Pixelization (HEALPIX, Górski et al. (2005) ) 666http://healpix.sourceforge.net pixelization scheme including all the simulated clouds as seen by an observer placed in the solar system. This map can be smoothed to match the resolution of a specific experiment and/or convolved with a realistic frequency bandwidth. When we compare with the Planck maps described in section 3, we convolve the simulated maps to the beam resolution of the 100 GHz channel ($\sim 10$ arcmin).

The procedure implemented for each realization is the following:

assign the $(R_{gal},\phi,z)$ Galacto-centric positions. In particular:

• •

$R_{gal}$ is extracted from a Gaussian distribution defined by the $R_{ring}$ and $\sigma_{ring}$ parameters. However, the $\sigma_{ring}$ is large enough to give non-zero probability at $R_{gal}\leq 0$. All of the negative values of $R_{gal}$ are either automatically set to $R_{gal}=0$ (axisymmetric case), or recomputed extracting new positive values from a normal distribution centred at $R=0$ and with the r.m.s given by the scale of the Galactic bar (spiral-arm case). This choice allows us to circumvent not only the issue of negative values of $R_{gal}$ due to a Gaussian distribution, but also to produce the high emissivity of the Central Molecular zone (see Ellsworth-Bowers et al. (2015) for a similar approach).

• •

the $z$-coordinate is drawn randomly from the distribution in Equation 2.

• •

the azimuth angle $\phi$ is computed from a uniform distribution ranging over $[0,2\pi)$ in the case of the axial symmetry. Conversely, in the case of spiral arms, $\phi$ follows the *logarithmic spiral * polar equation

[TABLE]

where $A=(\tan i)^{-1}$ and $B=-\log R_{bar}$ are, respectively, functions of the mean pitch angle and the starting radius of the spiral arm. In our case we set $i=-12$ deg, $R_{bar}=3$ kpc; 2. 2.

assign cloud sizes given the probability function $\mathcal{P}(L)$ (Equation 5); 3. 3.

assign emissivities to each cloud from the emissivity profile (see Equation 3); 4. 4.

convert $(R_{gal},\phi,z)$ positions into the *heliocentric * coordinate frame $(\ell,b,d_{\odot})$.

In Figure 3 we show an example of the 3D distribution of the emission as well as the distribution of the location of the simulated clouds using both of the geometries implemented in the package.

2.4 Simulation results

In Figure 4 we show two typical realizations of maps of CO emission for the Axisymmetric and LogSpiral geometries prior to any smoothing. As we are interested in the statistical properties of the CO emission, we report a few examples of the angular power spectrum $\mathcal{C}_{\ell}$ corresponding to different distributions of CO emission in Figure 5. In the ones shown subsequently the spectra are $\mathcal{D}_{\ell}$ encoding a normalization factor $\mathcal{D}_{\ell}=\ell(\ell+1)\mathcal{C}_{\ell}/2\pi$.

We can observe two main features in the morphology of the power spectrum: a bump around $\ell\sim 100$ and a tail at higher $\ell$. We interpret both of these features as the projection of the distribution of clouds from a reference frame off-centred (on the solar circle).

The bump reflects the angular scale ($\sim 1$ deg) related to the clouds which have the most likely size, parametrized by the typical size parameter, $L_{0}$, and which are close to the observer. On the other hand, the tail at $\ell\gtrsim 600$ (i.e. the arcminute scale) is related to the distant clouds which lie in the diametrically opposite position with respect to the observer. This is the reason why the effect is shifted to smaller angular scales. The $L_{0}$ and $\sigma_{ring}$ parameters modify the power spectrum in two different ways. For a given typical size, if the width of the molecular ring zone $\sigma_{ring}$ increases, the peak around $\ell\sim 100$ shifts towards lower multipoles, i.e. larger angular scales, and its amplitude increases proportionally to $\sigma_{ring}$, see for instance the bottom right panel in Figure 5. This can be interpreted as corresponding to the fact that the larger is $\sigma_{ring}$, the more likely it is to have clouds closer to the observer at the solar circle with a typical size given by $L_{0}$. On the other hand, if we choose different values for the size parameter (left panels in Figure 5) the tail at small angular scales moves downwards and flattens as $L_{0}$ increases. Vice versa, if we keep $L_{0}$ constant (Figure 5 bottom right panel), all of the tails have the same amplitude and an $\ell^{2}$ dependency. In fact, if $L_{0}$ is small, the angular correlation of the simulated molecular clouds looks very similar to the one of point sources which has Poissonian behaviour. Conversely if the typical size increases, the clouds become larger and they behave effectively as a coherent diffuse emission and less as point sources.

Far from the Galactic plane, the shape of the power spectrum is very different. In Figure 6 we show an example of the average power spectrum of 100 MC realizations of CO emission at high Galactic latitudes, i.e. $|b|>30$ deg, for both the Axisymmetric and LogSpiral geometries. For this run we choose the so-called best fit values for the $L_{0}$ and $\sigma_{ring}$ parameters discussed later in subsection 3.4. In addition to the different shape depending on the assumed geometry, one can notice a significant amplitude difference with respect to the power spectrum at low latitudes. Moreover, this is in contrast with the trend observed in the galactic plane, where the LogSpiral geometry tends to predict a power spectrum of higher amplitude. In both cases, however, the model suppresses the emission in these areas, as shown in Figure 4. In the LogSpiral case, the probability of finding clouds in regions in between spiral arms is further suppressed and could explain this feature. The emission is dominated by clouds relatively close to the observer for both geometries, and so the angular correlation is mostly significant at large angular scales (of the order of a degree or more) and is damped rapidly at small angular scales.

3 Comparison with Planck data

3.1 Dataset

The Planck collaboration released three different kinds of CO molecular line emission maps, described in Planck Collaboration et al. (2014a, 2016a). We decided to focus our analysis on the so-called Type 1 CO maps which have been extracted exploiting differences in the spectral transmission of a given CO emission line in all of the bolometer pairs relative to the same frequency channel. Despite being the noisiest set of maps, Type 1 are in fact the cleanest maps in terms of contamination coming from other frequency channels and astrophysical emissions. In addition, they have been obtained at the native resolution of the Planck frequency channels, and so allow full control of the effective beam window function for each map.

For this study we considered in particular the CO $1-0$ line, which has been observed in the 100 GHz channel of the HFI instrument. This channel is in fact the most sensitive to the CO emission in terms of signal-to-noise ratio (SNR) and the 1-0 line is also the one for which we have the most detailed external astrophysical observations. However, the Planck frequency bands were designed to observe the CMB and foreground emissions which gently vary with frequency and, thus, they do not have the spectral resolution required to resolve accurately the CO line emission. To be more quantitative, the Planck spectral response at 100 GHz is roughly 3 GHz, which corresponds to $\sim 8000\,\mathrm{km\,s^{-1}}$, i.e. about 8 orders of magnitude larger than the CO rotational line width (which can be easily approximated as a Dirac delta). Therefore, the CO emission observed by Planck  along each line of sight is integrated over the whole channel frequency band. Further details about the spectral response of the HFI instrument can be found in Planck Collaboration et al. (2014b).

3.2 Observed CO angular power spectrum

Since one of the goals of this paper is to understand the properties of diffuse CO line emission, we computed the angular power spectrum of the Type 1 $1-0$ CO map to compare qualitatively the properties of our model with the single realization given by the emission in our Galaxy. We distinguish two regimes of comparison, low Galactic latitudes ($|b|\leq 30\deg$) and high Galactic latitude ($|b|>30\deg$). While at low Galactic latitudes the signal is observed with high sensitivity, at high latitudes it is substantially affected by noise and by the fact that the emission in this region is faint due to its low density with respect to the Galactic disk.

In Figure 7 we show the angular power spectra of the first three CO rotational line maps observed by Planck  as well as the expected noise level at both high and low Galactic latitudes computed using a pure power spectrum estimator X2PURE (Grain et al., 2009). This is a pseudo power spectrum method (Hivon et al., 2002) which corrects the so called E-to-B-modes leakage in the polarization field that arises in the presence of incomplete sky coverage (Smith and Zaldarriaga, 2007; Bunn et al., 2003; Lewis et al., 2001). Although this feature is not strictly relevant for the analysis of this section, because we are considering the unpolarized component of the signal, it is important for the forecast presented in section 4. We estimated the noise as the mean of 100 MC Gaussian simulations based on the the diagonal pixel-pixel error covariance included in the Planck maps. One may notice how the noise has a level comparable to that of the CO power spectrum at high Galactic latitude. However, we note that the released Type 1 maps are obtained from the full mission data from Planck, and not from subsets of the data (e.g. using the so called half-rings or half-mission splits). Thus, it was not possible to test whether the observed flattening of the power spectrum at large angular scale is due to additional noise correlation not modelled by the Gaussian uncorrelated model discussed above. We notice that, if these maps were present, we could have had an estimate of this correlation using the noise given by the difference between the map auto-spectra and the noise-bias free signal obtained from the cross-spectra of the maps from data subsets. Since even for the 1-0 line, the noise becomes dominant on scales $\ell\approx 20$ we decided to limit the comparison at low Galactic latitude where the signal to noise ratio is very high.

We note that in the following we considered the error bars on the power spectrum as coming from the gaussian part of the variance, i.e., following Hivon et al. (2002)

[TABLE]

where $\nu$ is the number of degrees of freedom taking into account the finite number of modes going into the power spectrum calculation in each $\ell$ mode and the effective sky coverage. $N_{\ell}$ represents the noise power spectrum and the $C_{\ell}$ is the theoretical model describing the CO angular power spectrum with the tilde denoting measured quantities. Because we do not know the true CO theoretical power spectrum we assumed that $C_{\ell}+N_{\ell}=\tilde{C}_{\ell}$. The gaussian approximation however underestimates the error bars. The CO field is in fact a highly non-gaussian field with mean different from zero. As such, its variance should contain contributions coming from the expectation value of its 1 and 3 point function in the harmonic domain that are zero in the gaussian approximation. These terms are difficult to compute and we considered the gaussian approximation sufficient for the level of accuracy of this study.

As can be seen in Figure 7, all of the power spectra of CO emission at low Galactic latitudes have a broad peak around the multipole $100\div 300$, i.e. at the $\approx 1\deg$ angular scale. The signal power starts decreasing up to $\ell\sim 600$ and then grows again at higher $\ell$ due to the Planck instrumental noise contamination. Such a broad peak suggests that there is a correlated angular scale along the Galactic plane. This can be understood with a quick order of magnitude estimate. If we assume that most of the CO emission is localized at a distance of 4 kpc (in the molecular ring) and molecular clouds have a typical size of 30 pc, we find that each cloud subtends a $\sim 0.5$ deg area in the sky. This corresponds to a correlated scale in the power spectrum at an $\ell$ of the order of a few hundred but the detail of this scale depends on the width of the molecular ring zone.

3.3 Galactic plane profile emission comparison

As a first test we compared the profile of CO emission in the Galactic plane predicted by the model and the one observed in the data. Since we are mostly interested in a comparison as direct as possible with the Planck observed data, we convolved the MCMole3D maps with a Gaussian beam of 10 arcmin FWHM, corresponding to the nominal resolution of the 100 GHz channel of HFI, prior to any further processing.

In order to compare the data and the simulations, we constrained the total flux of the simulated CO maps with the one observed in the Planck data. This is necessary, otherwise the emission would be directly proportional to the number of clouds distributed in the simulated Galaxy. Such a procedure also breaks possible parameter degeneracies with respect to the amplitude of the simulated power spectra (see next section). Following Bronfman et al. (1988), we therefore computed the integrated flux of the emission along the two Galactic latitudes and longitudes $\left(l,b\right)$ defined as

[TABLE]

where $X$ refers both to the model and to the observed CO map. We then rescaled the simulated maps, dividing by the factor $f$ defined as:

[TABLE]

We estimated the integrals in Equation 7 and Equation 8 by considering a narrow strip of Galactic latitudes within $\left[-2,2\right]$ degrees. We found that the value of $f$ is essentially independent of the width of the Galactic latitude strip used to compute the integrals because most of the emission comes from a very thin layer along the Galactic plane of amplitude $|b|\lesssim 2$ deg.

In Figure 8 we show the comparison between $I^{observ}(l)$ and the $I^{model}(l)$ as defined in Equation 7 computed as the mean of 100 MC realizations of galaxies populated by molecular clouds for both the Axisymmetric and LogSpiral models as well as their typical standard deviation. In particular, we chose for these simulations the default parameters in Table 1.

The emission profiles are quite consistent in the regions from which most of the CO emission comes, i.e. in the inner Galaxy, the I and the IV quadrants (longitude in $\left[-90,90\right]$ deg777We stress that the definition of quadrants comes from the Galactic coordinates centred on the Sun. The I and IV quadrants are related to the inner Galaxy, while the II and the III ones look at its outer regions.). On the contrary, the emission in the other two quadrants looks to be under-estimated but within the scatter of the simulations. In fact, the observed emissions in both the II and III quadrants come mainly from the closer and more isolated system of clouds. These are actually more difficult to simulate because in that area (at Galactic longitudes $|l|>100$ deg) the presence of noise starts to be non-negligible (see shaded blue in Figure 8).

In addition, we note that the bump in the profile at $l\simeq 60-70$ deg, where we see a lack of power in both the Axisymmetric and LogSpiral cases, corresponds to the complex region of Cygnus-X, which contains the very well known X-ray source Cyg-X1, massive protostars and one of the most massive molecular clouds known, $3\times 10^{6}\mathrm{M}_{\odot}$, 1.4 kpc distant from the Sun. Given the assumptions made in section 2, these large and closer clouds are not easy to simulate with MCMole3D especially where they are unlikely to be found, as in inter-spiral arm regions. Despite of this, one can notice an overall qualitative better agreement with observations for the LogSpiral model than for the Axisymmetric one. The latter reconstructs the global profile very well, but the former contains more peculiar features such us the central spike due to the Central Molecular Zone within the bar, or the complex of clouds at longitudes around $\sim-140,\,-80,\,120$ deg. We will perform a more detailed comparison of the two geometries in the following section and in Appendix A.

3.4 Constraining the MCMole3D model with Planck data

After comparing the CO profile emission we checked whether the MCMole3D model is capable of reproducing the characteristic shape of the Planck CO angular power spectrum. Given the knowledge we have on the shape of the Milky Way, we decided to adopt the LogSpiral geometry as a baseline for this comparison, and to fix the parameters for the specific geometry to the values describing the shape of our Galaxy (see subsection 2.3). For sake of completeness we reported the results of the same analysis adopting an Axisymmetric geometry in Appendix A.

We left the typical cloud size $L_{0}$ and $\sigma_{ring}$ (the width of the molecular ring) as free parameters of the model. While the former is directly linked to the observed angular size of the clouds, the role of the second one is not trivial, especially if we adopt the more realistic 4 spiral arms distribution. Intuitively, it changes the probability of observing more clouds closer to the observer and affects more the amplitude of the power on the larger angular scales.

We defined a large interval, reported in Table 1, where $L_{0}$ and $\sigma_{ring}$ are allowed to vary. Looking at the series of examples reported in Figure 5 we can see that suitable parameter ranges which yield power spectra close to the Planck observations are $L_{0}=10\div 30$ pc and $\sigma_{ring}=2\div 3$ kpc. It is interesting to note that these are in agreement with estimates available in the literature (see e.g. (Ellsworth-Bowers et al., 2015; Roman-Duval et al., 2010)).

We then identified a set of values within the intervals just mentioned for which we computed the expected theoretical power spectrum of the specific model. Each theoretical model is defined as the mean of the angular power spectrum of 100 MC realizations of the model computed with X2PURE. For each realization of CO distribution we rescaled the total flux following the procedure outlined in the previous section before computing its power spectrum.

Once the expected angular power spectra for each point of the parameter domain had been computed, we built the hyper-surface $\mathcal{F}(\ell;\sigma_{ring},L_{0})$ which for a given set of values $(\sigma_{ring},L_{0})$ retrieved the model power spectrum, by interpolating it from its value at the closest grid points using splines. We checked that alternative interpolation methods did not impact significantly our results. We then computed the best-fit parameters of the MCMole3D model by performing a $\chi^{2}$ minimization with the Planck CO power spectrum data. For this procedure we introduced a further global normalization parameter $A_{CO}$ to take account of the Planck bandpass effects or other possible miscalibration of the model. These might come either from variations from the scaling laws employed in the model (that are thus not captured by the total flux normalization described earlier), or calibration differences between the Planck data and the surveys used to derive the scaling laws themselves. The bandpass effects tend to decrease the overall amplitude of the simulated signal because each line gets diluted over the width of the Planck frequency band.

Since the theoretical model has been estimated from Monte Carlo simulations, we added linearly to the sample variance error of the Planck data an additional uncertainty budget corresponding to the uncertainty of the mean theoretical power spectrum estimated from MC. We note that when we compute the numerator of the $f$ rescaling factor, we include not only the real flux coming from the CO lines but also an instrumental noise contribution. We therefore estimated the expected noise contribution to $f$ by computing the integral of Equation 8 on the Planck error map and found it to be equal to 10%. We propagated this multiplicative uncertainty to the power spectrum level, rescaling the mean theoretical MC error bars by the square of this factor.

We limit the range of angular scales involved in the fit to $\ell\leq 400$ in order to avoid the regions that display an unusual bump at scales of around $\ell\approx 500$ that is not captured by any realization of our model (see next section). The best-fit parameters are reported in Table 2

[TABLE]

The values are within the ranges expected from the literature. As can be seen in Figure 9 the power spectrum corresponding to the model with the best fit parameters, describes the Planck data reasonably well. The minimum $\chi^{2}$ obtained by the minimization process gives $1.48$ that corresponds to a $p$-value of $13\%$. We note, however,that all of the parameters are highly correlated. This is somewhat expected as the larger is $\sigma_{ring}$, the closer the clouds get to the observer placed in the solar circle. This effect can be compensated by an overall decrease of the typical size of the molecular cloud as shown in Figure 5(d).

Finally, we note that $A_{CO}\lesssim 1$ suggests that, despite the rescaling procedure constraining quite well the overall power spectrum amplitude, the spatial distribution seems to be more complex than the one implemented in the model. This might partially be explained by the fact that we do not model explicitly any realistic bandpass effect of the Planck channel or the finite width of the CO line. Additional sources of signal overestimation could be residual contamination of ${}^{13}CO$ 1-0 line or thermal dust in the map or variations of the emissivity profile in Equation 3.

3.5 Consistency checks on other maps

The Planck collaboration released multiple CO maps extracted using different component separation procedures. We can test the stability of our results by using CO maps derived with these different approaches, in particular the so-called Type 2 maps. These have been produced exploiting the intensity maps of several frequencies (multi-channel approach) to separate the CO emission from the astrophysical and CMB signal (Planck Collaboration et al., 2014a). The maps are smoothed at a common resolution of $15$ arcmin and have better S/N ratio than the Type 1 ones. However, the CO is extracted by assuming several simplifications which may leak into contamination due to foreground residuals and systematics, as explained in Planck Collaboration et al. (2016a, section 5.5.3).

We repeated the procedure outlined in subsection 3.3 and subsection 3.4 using the Type 2 $1-0$ map. The values of the best fit parameters are summarized in Table 2 and we show in Figure 9 the best-fit model power spectrum together with the power spectrum of Type 2 map data. We found that the values of $A_{CO}$ obtained for Type 2 are inconsistent with the one obtained for the Type 1 maps. However, this discrepancy is consistent with the overall inter calibration difference between the two maps reported in Planck Collaboration et al. (2016a). Such differences are mainly related to a combination of bandpass uncertainties in the Planck observations and presence of a mixture of ${}^{12}CO$ and ${}^{13}CO$ (emitted at $110$ GHz) lines for the Type 1 maps. While $\sigma_{ring}$ is consistent between the two maps, the Type 2 $L_{0}$ parameters are in slight tension at $2.7\sigma$ level. The overall correlation of the parameters is increased and the overall agreement between data and the MCMole3D mode is reduced although it remains acceptable. We cannot exclude however that this is a sign of additional systematic contamination in the Type 2 maps.

The Planck collaboration provided maps of the 2-1 line for both of the methods and we could use our model to constrain the relative amplitudes of the lines, while fixing the parameter of the cloud distribution. However, such analysis is challenging and might be biased by the presence of variations of local physical properties of the clouds (opacity and temperature) or by the red or blue shift of the CO line within the Planck bandpass induced by the motion of the clouds themselves (Planck Collaboration et al., 2016a). For these reasons, we decided to restrict our analysis only to the CO $1-0$ line, since it is the one for which the observational data are more robust.

We finally note that the observed angular power spectra of the Planck maps display an oscillatory behaviour at a scale of $\ell\geq 400$ with a clear peak at around $\ell\approx 500$. The fact that this feature is present in all of the lines and for all of the CO extraction methods means that it can reasonably be considered as a meaningful physical signature. Because a single cloud population produces an angular power spectrum with a characteristic peak scale, we speculate that this could be the signature of the presence of an additional cloud population with a different typical size or location. We however decided to leave the investigation of this feature for a future work.

3.6 Comparison with data at high Galactic latitudes

In Figure 10, we compare the Planck CO 1-0 power spectrum at high Galactic latitudes with the average power spectrum of 100 MC realizations of the MCMole3D model for the same region of the sky. We assumed for these runs the best-fit values of the $L_{0},\,\sigma_{ring}$ parameters reported in Equation 10 and a LogSpiral distribution. Because the Planck maps at these latitudes are dominated by noise, we subtracted our MC estimates of the noise bias data power spectrum so as to have a better estimate of the signal (blue circles).

As can be observed in Figure 10, some discrepancy arises when comparing the power spectrum expected from the simulation of LogSpiral MCMole3D at high Galactic latitudes with the noise debiased Planck data. This is rather expected because the model has larger uncertainties at high Galactic latitudes than in the Galactic mid-plane (where the best-fit parameters are constrained) given the lack of high sensitivity data. The discrepancy seems to point to an overestimation of the vertical profile parameters $\sigma_{z,0}$ and $h_{R}$ (see Equation 2) which gives a higher number of clouds close to the observer at high latitude. However, we also point out that, as explained in subsection 3.2, the error bars in Figure 10 might be underestimated especially at the largest angular scales where we are signal dominated. Therefore, discrepancies of order $\approx 3\sigma$ do not seem unlikely. Since we are mainly interested in using the model to forecast the impact of unresolved CO emission far from the Galactic plane ($|b|>30$), we investigated whether removing the few high Galactic latitude clouds in the simulation that appear close to the observer would improve the agreement with the Planck data. All of these clouds have, in fact, a flux exceeding the Planck CO map noise in the same sky area and they should have already been detected in real data. We will refer to this specific choice of cut as the High Galactic Latitudes (HGL) cut in the following. The power spectrum of the MCMole3D simulated maps obtained after the application of the HGL cut is shown in Figure 10. We found that the model calibrated at low latitudes and after the application of the HGL-cut agrees very well with the data on the angular scales where the signal slightly dominates, i.e. $\ell\lesssim 80$. We could not extend the comparison to smaller angular scales because the data become noise dominated and the residual increase of power observed on the power spectrum is dominated by a noise bias residual.

4 Polarization Forecasts

As noted in section 1, CO lines are polarized and could contaminate sensitive CMB polarization measurements together with other polarized Galactic emission (synchrotron and the thermal dust) at sub-millimeter wavelengths. Future experiments will preferentially perform observations at intermediate and high Galactic latitudes, to minimize contamination from strong Galactic emissions close to the plane. Since CO data at high Galactic latitudes are not sensitive enough to perform accurate studies of this emission, we provide two complementary estimates of the possible contamination from its polarized counterpart to the CMB B-mode power spectrum in this sky region.

4.1 Data-based order of magnitude estimate

Starting from the measured Planck power spectrum at low Galactic latitudes, one can extrapolate a very conservative value of the CO power spectrum at higher latitudes. Assuming that all of the variance observed in the high Galactic latitude region is distributed among the angular scales in the same way as in the Galactic plane, we can write

[TABLE]

This is a somewhat conservative assumption because we know that the bulk of the CO line emission is concentrated close to the Galactic disk and also because it assumes that the Planck noise at high Galactic latitudes is diffuse CO emission. The variance of the Planck CO map is $0.3\,\mathrm{K^{2}(km\,s^{-1})^{2}}$, at $|b|>30$ deg, while for $|b|<30$ deg we get a variance of $193.5$ $\mathrm{K^{2}\,(km\,s^{-1})^{2}}$. Taking $1\%$ as the polarization fraction, $p_{CO}$, of the CO emission and an equal power in E and B-modes of polarized CO, we can convert $\mathcal{C}_{\ell}^{COhigh}$ into its B-mode counterpart as $\mathcal{C}_{\ell}^{COhigh,EE}=\mathcal{C}_{\ell}^{COhigh,BB}=\mathcal{C}_{\ell}^{high,CO}p_{CO}^{2}/2$. We then apply the conversion factors of Planck Collaboration et al. (2014a) to convert the CO power spectrum into thermodynamic units (from $\mathrm{K_{RJ}kms^{-1}}$ to $\mu\mathrm{K}$). We can compare $\mathcal{C}_{\ell}^{COhigh,BB}$ to the amplitude of equivalent cosmological CMB inflationary B-modes with tensor-to-scalar ratio $r=1$ at $\ell={80}$. In terms of $\mathcal{D}^{BB}_{\ell}$, this is equal to $\sim 6.67\times 10^{-2}\mu\mathrm{K}^{2}$ for a fiducial Planck 2015 cosmology. We found that the amplitude of the extrapolated CO B-mode power spectrum is equal to a primordial B-mode signal having $r_{CO}=0.025$.

4.2 Simulation estimate

In order to verify and refine the estimate given in the previous Section, we used the model presented in section 2 to infer the level of contamination from unresolved polarized CO emission. For doing this, we first set the free parameter of the MCMole3D model to the best-fit value derived in Equation 10.

From the total unpolarized emission in each sky pixel of the simulation, $I^{CO}$ we can then extract its linearly polarized part by taking into account the global properties of the Galactic magnetic field. Following Delabrouille et al. (2013); Tassis and Pavlidou (2015) the Q and U Stokes parameters of each CO cloud can be related to the unpolarized emission as

[TABLE]

where $p_{CO}$ is the intrinsic polarization fraction of the CO lines, while $g_{d}$ is the geometric depolarization factor which accounts for the induced depolarization of the light when integrated along the line of sight. The polarization angle $\psi$ describes the orientation of the polarization vector and, for the specific case of Zeeman emission, it is related to the orientation of the component of the Galactic magnetic field orthogonal to the line of sight $B_{\bot}$. Following the findings of Greaves et al. (1999), we adopted a conservative choice of a constant $p_{CO}=1\%$ for each molecular cloud of the simulation. Because the polarized emission in molecular clouds is correlated with the polarized dust emission (Crutcher, 2012, see, e.g.), we used the $g_{d}$ and $\psi$ templates for the Galactic dust emission available in the public release of the Planck Sky Model suite888http://www.apc.univ-paris7.fr/~delabrou/PSM/psm.html (Delabrouille et al., 2013). These have been derived from 3D simulations of the Galactic magnetic field (including both a coherent and a turbulent component) and data of the WMAP satellite.

Since we assumed a constant polarization fraction, the geometrical depolarization effectively induces a change in the polarization fraction as a function of Galactic latitudes decreasing it when moving away from the poles. This effect has already been confirmed by Planck observations (Planck Collaboration et al., 2016a) of thermal dust, whose polarization fraction increases at high latitudes.

In order to forecast the contamination of unresolved CO polarized emission alone, we apply the HGL-cut as described in subsection 3.6 to each realization of the model for consistency.

Once the $Q^{CO}$ and $U^{CO}$ maps have been produced, we computed the angular power spectrum using X2PURE.

In Figure 12 we show the mean and standard deviation of the B-mode polarization power spectrum extracted from 100 MC realizations of the CO emission following the procedure just outlined. Even though in subsection 3.4 we showed that our model tends to slightly overestimate the normalization of the power spectrum, we decided not to apply the best-fit amplitude $A_{CO}$ to the amplitude of the B-mode power spectrum in order to provide the most conservative estimates of the signal.

As could be seen from Figure 12, there is a significant dispersion compared to the results of the MC simulations at low Galactic latitudes (see Figure 5). This simply reflects the fact that the observations, and hence our model, do not favour the presence of molecular clouds at high Galactic latitudes. Therefore their number can vary significantly between realizations. We repeated this test using the Axisymmetric geometry and changing the parameter $\sigma_{ring}$. The result is stable with respect to these assumptions. We found that the spatial scaling of the average E and B-mode power spectrum can be approximated by a decreasing power-law $\mathcal{D}_{\ell}\sim\ell^{\alpha}$, with $\alpha=-1.78$.

Our simulations suggest that the level of polarized CO emission from unresolved clouds, despite being significantly lower than synchrotron or thermal dust, can nevertheless significantly bias the primordial B-mode signal if not taken into account. The signal concentrates mainly on large angular scales and at $\ell\sim 80$, $\mathcal{D}_{\ell}=(1.1\pm 0.8)\times 10^{-4}\mu\mathrm{K}^{2}$ where the uncertainty corresponds to the error in the mean spectra estimated from the 100 MC realizations. Therefore, the level of contamination is comparable to a primordial B-mode signal induced by tensor perturbations of amplitude $r_{CO}=0.003\pm 0.002$, i.e. below the recent upper limit $r<0.07$ reported by the BICEP2 Collaboration (BICEP2 Collaboration et al., 2016) but higher than the $r=0.001$ target of upcoming experiments (Abazajian et al., 2016; Matsumura et al., 2014; CORE Collaboration et al., 2016). The contamination quickly becomes sub-dominant on small angular scales ($\ell\approx 1000$) where the B-modes are mostly sourced by the gravitational lensing.

We finally note that these estimates are conservative since the assumed polarization fraction of $1\%$ of polarized is close to the high end of the polarization fractions observed in CO clouds.

5 Conclusions

In this work we have developed a parametric model for CO molecular line emission which takes account of the CO clouds distribution within our Galaxy in 3D with different geometries, as well as the most recent observational findings concerning their sizes, locations, and emissivity.

Despite most of the observations having so far been confined to the Galactic plane, we have built the model to simulate the emission over the full sky. The code implementing MCMole3D is being made publicly available.

We have compared the results of our simulations with Planck CO data on the map level and statistically (by matching angular power spectra). We found that:

the parameters of the size function, $L_{0}$, and the width of the Galactic radial distributions $\sigma_{ring}$ play a key role in shaping the power spectrum; 2. 2.

the choice of symmetries in the cloud distribution changes the profile of the integrated emission in the Galactic plane (Equation 7) but not the power spectrum morphology; 3. 3.

our model is capable of reproducing fairly well the observations at low Galactic latitudes (see Figure 8) and the power spectrum at high latitudes (Figure 10).

We used our model to fit the Planck observed CO power spectrum and to estimate the most relevant parameters of the CO distribution, such as the typical size of clouds and the thickness of the molecular ring, finding results in agreement with values reported in the literature. The model which we have developed could easily be generalized and extended whenever new data become available. In particular, its accuracy at high Galactic latitudes would greatly benefit from better sub-mm measurements going beyond the Planck sensitivity, as well as from better information about the details of the CO polarization properties.

Finally, we used the best-fit parameters obtained from comparing the MCMole3D model with Planck data to forecast the unresolved CO contamination of the CMB  B-mode power spectrum at high Galactic latitudes. We conservatively assumed a polarization fraction of $p_{CO}=1\%$, which corresponds to the high end of those observed at low latitudes,since no polarized CO cloud has yet been observed far from the Galactic plane due to the weakness of this emission.

We found that this signal could mimic a B-mode signal with tensor-to-scalar ratio $0.001\lesssim r\lesssim 0.025$. This level of contamination is indeed relevant for accurate measurements of CMB B-modes. It should therefore be inspected further in light of the achievable sensitivities of upcoming and future CMB experiments together with the main diffuse polarized foreground (thermal dust and synchrotron). From the experimental point of view, trying to find dedicated instrumental solution for minimizing the impact of CO emission lines, appears to be particularly indicated in the light of these results.

Acknowledgements. We would like to thank Françoise Combe for many useful comments and suggestions for the development of this study, as well as Alessandro Bressan, Luigi Danese, Andrea Lapi, Akito Kusaka, Davide Poletti and Luca Pagano for several enlightening discussions. We thank John Miller for his careful reading of this work. We thank Guillaume Hurier for several clarifications about the Planck CO products. This work was supported by the RADIOFOREGROUNDS grant of the European Union’s Horizon 2020 research and innovation programme (COMPET-05-2015, grant agreement number 687312) and by by Italian National Institute of Nuclear Physics (INFN) INDARK project. GF acknowledges support of the CNES postdoctoral programme. Some of the results in this paper have been derived using the HEALPIX (Górski et al., 2005) package.

Appendix A Best-fit with Axisymmetric geometry

In this appendix we present the results of the analysis described in subsection 3.4 to constraint the CO distribution using the MCMole3D model adopting an Axisymmetric geometry instead of the LogSpiral one. Following the procedure of subsection 3.4 we construct a series of $\mathcal{F}(\ell;\sigma_{ring},L_{0})$ hyper-surfaces sampled on an ensemble of specific values of the $L_{0}$ and $\sigma_{ring}$ parameters within the same ranges reported in subsection 3.4.

In Figure 13 we show the results of the fit of the axisymmetric MCMole3D model to the CO power spectrum of the Planck Type 1 and Type 2 CO maps in the Galactic plane. We summarize the best-fit values of these parameters in Table 3. As it can be seen from the results of the $\chi^{2}$ test in Table 3 the Axisymmetric model does not fit the data satisfactorily. Moreover one of the parameters of the model, the typical cloud size $L_{0}$, is in practice unconstrained. For this reason we decided to adopt the LogSpiral geometry as a baseline choice for our forecast presented in section 4. Nevertheless, we pushed the comparison between the two geometries in the high galactic latitude area for sake of completeness.

In Figure 14 we show the comparison between the Planck data for Type 1 maps and the MCMole3D axisymmetric best-fit model after the application of the HGL cut described in the paper. The Axisymmetric model describes the data similarly to the LogSpiral model at the larger scales. The difference in the signal amplitude is in fact less then 30% for angular scales $\ell\lesssim 100$ and the two models are compatible within the error bars. This seems to indicate that in this regime, the details of the CO distribution in the high galactic latitude region are mainly affected by the properties of the vertical profile rather than by the geometry of the distribution. Conversely, the difference between the two geometries becomes important at smaller angular scales reaching a level of $\approx 2$ at $\ell\approx 1000$.

We finally performed a series of polarized simulations as in subsection 4.2 to access the level of contamination to the CMB B-modes power spectrum with the best-fit Axisymmetric model and found $r_{CO}\lesssim 0.001$. Moreover, the slope of the BB power spectrum in Figure 14(b) is $-2.2$ similar to the one computed with the LogSpiral geometry.

Because the LogSpiral model describes the data both in the high and low galactic latitude area, we consider the upper limit derived with this setup more reliable and the reference estimate for the contamination to the cosmological signal due to the CO polarized emission.