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

G. Puglisi,1,2 G. Fabbian,3,1,2 C. Baccigalupi1,2
1SISSA- International School for Advanced Studies, Via Bonomea 265, 34136
Trieste, Italy
2INFN-National Institute for Nuclear Physics, Via Valerio 2, I-34127 Trieste,
Italy
3Institut d’Astrophysique Spatiale, CNRS (UMR 8617), Univ. Paris-Sud,
Université Paris-Saclay, bât. 121, 91405 Orsay, France
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### 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.

###### keywords:

Interstellar Medium: molecules, magnetic fields, lines and bands Cosmology:
observations, cosmic background radiation

††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.

Figure 1:  Planck CO $1-0$ map (Planck Collaboration et al., 2016a). Note the
predominance of instrumental noise in regions far from the Galactic plane.

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:

$z_{1/2}(R)\propto\sigma_{z}(R)=\sigma_{z,0}\cosh\left(\frac{R}{h_{R}}\right),$
(1)

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:

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

Figure 2: (top) Histograms of $dN/dL$ computed by assigning the size of each
cloud with the probability function (bottom). The two spectral indices
$\alpha_{L}\approx 3.3\left(3.9\right)$ refer respectively to clouds in the
inner (outer) Galaxy.

(a)

(b)

(c)

(d)

Figure 3: _Top panels_ : Density contour plots of an MC galaxy population with
40,000 objects distributed following the (a) LogSpiral and (b) Axisymmetric
distributions. _Bottom panels_ : Probability Density Function (PDF) of 100 MC
realizations of 40,000 molecular clouds following the (c) LogSpiral and (d)
Axisymmetric geometry. The latter case is consistent with results in
Ellsworth-Bowers et al. (2015).

###  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):

$\epsilon_{0}(R)=\epsilon_{c}\exp\left(R/R_{em}\right),$ (3)

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)

$\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.$ (4)

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$:

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

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.

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}\,^{{\dagger}}$ [kpc] | 3
$i\,^{{\dagger}}$ [deg] | -12
$\epsilon_{c}\,\left[\mathrm{K_{RJ}\,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 1: List of parameters used in MCMole3D simulations. † only for
LogSpiral.

###  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:

$\displaystyle\Theta^{CO}=\\{$ $\displaystyle
N_{clouds},\,\epsilon_{c},\,R_{em},\,R_{ring},\,\sigma_{ring},$
$\displaystyle\sigma_{z,0},\,h_{R},\,L_{min},L_{max},L_{0}\\}.$

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

Figure 4: Two realizations of CO maps simulated with MCMole3D using the
distribution parameters given by the values 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:

  1. 1.

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

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

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.

(a)

(b)

(c)

(d)

Figure 5: Angular power spectra of CO emission in the Galactic plane computed
for 100 MC realizations of the MCMole3D model assuming different values of its
free parameters. The mean value of the simulation is shown by solid, dashed
and dot-dashed lines while the shaded area represents the measured variance of
the realizations. The top row shows the case of an Axisymmetric geometry while
the bottom panel displays results for a LogSpiral geometry. Results obtained
by varying the $L_{0}$ ($\sigma_{ring}$) parameters are shown on the left
(right) column.

###  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$.

Figure 6: Examples of the power spectra of CO emission at high Galactic
latitudes ($|b|>30$ deg) for Axisymmetric and LogSpiral geometries. For both
the geometries we assumed the best-fit values of the parameters describing the
CO distribution presented in subsection 3.4.

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)

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

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.

Figure 7: CO $1-0$ angular power spectrum (blue solid) estimated from the
Planck map at low (top) and high (bottom) Galactic latitudes. The shaded area
shows the error bar due to the sample and noise variance. The expected noise
level of the maps in the two regions is shown in (red dashed).

###  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

$\displaystyle I^{X}(l)=\int dbI^{X}(l,b),$ (7) $\displaystyle
I_{tot}^{X}=\int dldbI^{X}(l,b),$ (8)

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:

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

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.

(a)

(b)

Figure 8:  Comparison of 100 MC realization of MCMole3D simulated Galactic CO
emission profiles (dashed green) with Planck observations (solid blue). The
average profile integrals defined in (7) for the Axisymmetric (a) and
LogSpiral (b) geometries are shown with dashed lines. The shaded area displays
the standard deviation of all MCs in each longitude bin (green), or the noise
level of Planck (solid blue) estimated from the _Type 1_ null map.

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.

(a)

(b)

(c)

(d)

Figure 9: _Top panels_ : angular power spectra (solid thin blue) of the Planck  _Type 1_ (left) and _Type 2_ (right) maps. The shaded area correspond to the $1\sigma$ (dark blue) and $2\sigma$ (light blue) error bars including statistical and systematic uncertainties. The MCMole3D model CO angular power spectrum assuming the best-fit parameters of Equation 10 is shown in thick solid magenta. The Planck noise power spectrum is shown in red dashed, in the top right panel the noise level is about one order of magnitude smaller than the one in the top left panel. _Bottom panels_ : best-fit parameters of the MCMole3D model describing the Planck CO angular power spectrum and their correlations. |  $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 2: Summary table of best fit parameters obtained using the two different
Planck CO maps.

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

$\displaystyle L_{0}=14.50\pm 0.58\,\mathrm{pc},$
$\displaystyle\sigma_{ring}=2.76\pm 0.19\,\mathrm{Kpc},$ (10) $\displaystyle
A_{CO}=0.69\pm 0.06.$

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.

Figure 10:  CO 1-0 power spectrum at high Galactic latitudes of the LogSpiral
MCMole3D model (dotted green), using the parameters in (10). (thick solid
orange) We show the power spectrum for the LogSpiral geometry, the same
parameters in (10) and with the HGL-cut of clouds at $|b|>30$ deg whose flux
exceeds the Planck noise. The Planck  _Type 1_ CO power spectrum before and
after noise bias subtraction is shown with the blue solid line and filled
circles respectively; the error bars account for both Planck data statistical
uncertainties and systematics from the MCMole3D simulations. The noise bias is
shown with the dashed red line.

##  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

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

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

$\displaystyle Q(\hat{n})^{CO}$
$\displaystyle=p_{CO}\,g_{d}(\hat{n})I(\hat{n})^{CO}\,cos(2\psi(\hat{n})),$
(12) $\displaystyle U(\hat{n})^{CO}$
$\displaystyle=p_{CO}\,g_{d}(\hat{n})I(\hat{n})^{CO}\,sin(2\psi(\hat{n}))..$
(13)

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.

(a)

(b)

Figure 11: (a) Polarization angle and (b) geometrical depolarization maps used
for simulating polarized CO emission in this work. Figure 12:  B-mode power
spectra of polarized CO emission lines at high Galactic latitudes estimated
using the best-fit parameters of the LogSpiral MCMole3D model (see Equation
10). The expected Planck 2015 $\Lambda$CDM cosmological signal including the
gravitational lensing contribution is shown in black. Potential contributions
from inflationary B-modes for tensor-to-scalar ratios of $r=0.1,0.01,0.001$
are shown with solid, dot-dashed and dashed green lines respectively. The red
arrow indicates the upper limits obtained in subsection 4.1

##  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:

  1. 1.

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.

## References

  * Abazajian et al. [2016] K. N. Abazajian, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, K. S. Arnold, C. Baccigalupi, J. G. Bartlett, N. Battaglia, B. A. Benson, C. A. Bischoff, J. Borrill, V. Buza, E. Calabrese, R. Caldwell, J. E. Carlstrom, C. L. Chang, T. M. Crawford, F.-Y. Cyr-Racine, F. De Bernardis, T. de Haan, S. di Serego Alighieri, J. Dunkley, C. Dvorkin, J. Errard, G. Fabbian, S. Feeney, S. Ferraro, J. P. Filippini, R. Flauger, G. M. Fuller, V. Gluscevic, D. Green, D. Grin, E. Grohs, J. W. Henning, J. C. Hill, R. Hlozek, G. Holder, W. Holzapfel, W. Hu, K. M. Huffenberger, R. Keskitalo, L. Knox, A. Kosowsky, J. Kovac, E. D. Kovetz, C.-L. Kuo, A. Kusaka, M. Le Jeune, A. T. Lee, M. Lilley, M. Loverde, M. S. Madhavacheril, A. Mantz, D. J. E. Marsh, J. McMahon, P. D. Meerburg, J. Meyers, A. D. Miller, J. B. Munoz, H. N. Nguyen, M. D. Niemack, M. Peloso, J. Peloton, L. Pogosian, C. Pryke, M. Raveri, C. L. Reichardt, G. Rocha, A. Rotti, E. Schaan, M. M. Schmittfull, D. Scott, N. Sehgal, S. Shandera, B. D. Sherwin, T. L. Smith, L. Sorbo, G. D. Starkman, K. T. Story, A. van Engelen, J. D. Vieira, S. Watson, N. Whitehorn, and W. L. Kimmy Wu.  CMB-S4 Science Book, First Edition.  _ArXiv e-prints_ , October 2016.
  * Baccigalupi [2003] C. Baccigalupi.  Cosmic microwave background polarisation: foreground contrast and component separation.  _New Astron. Rev._ , 47:1127–1134, December 2003.  doi: 10.1016/j.newar.2003.09.038.
  * Bally et al. [1987] J. Bally, W. D. Langer, A. A. Stark, and R. W. Wilson.  Filamentary structure in the Orion molecular cloud.  _ApJ_ , 312:L45–L49, January 1987.  doi: 10.1086/184817.
  * Balser et al. [2011] D. S. Balser, R. T. Rood, T. M. Bania, and L. D. Anderson.  H II Region Metallicity Distribution in the Milky Way Disk.  _ApJ_ , 738:27, September 2011.  doi: 10.1088/0004-637X/738/1/27.
  * BICEP2 Collaboration et al. [2016] BICEP2 Collaboration, Keck Array Collaboration, P. A. R. Ade, Z. Ahmed, R. W. Aikin, K. D. Alexander, D. Barkats, S. J. Benton, C. A. Bischoff, J. J. Bock, R. Bowens-Rubin, J. A. Brevik, I. Buder, E. Bullock, V. Buza, J. Connors, B. P. Crill, L. Duband, C. Dvorkin, J. P. Filippini, S. Fliescher, J. Grayson, M. Halpern, S. Harrison, G. C. Hilton, H. Hui, K. D. Irwin, K. S. Karkare, E. Karpel, J. P. Kaufman, B. G. Keating, S. Kefeli, S. A. Kernasovskiy, J. M. Kovac, C. L. Kuo, E. M. Leitch, M. Lueker, K. G. Megerian, C. B. Netterfield, H. T. Nguyen, R. O’Brient, R. W. Ogburn, A. Orlando, C. Pryke, S. Richter, R. Schwarz, C. D. Sheehy, Z. K. Staniszewski, B. Steinbach, R. V. Sudiwala, G. P. Teply, K. L. Thompson, J. E. Tolan, C. Tucker, A. D. Turner, A. G. Vieregg, A. C. Weber, D. V. Wiebe, J. Willmert, C. L. Wong, W. L. K. Wu, and K. W. Yoon.  Improved Constraints on Cosmology and Foregrounds from BICEP2 and Keck Array Cosmic Microwave Background Data with Inclusion of 95 GHz Band.  _Physical Review Letters_ , 116(3):031302, January 2016.  doi: 10.1103/PhysRevLett.116.031302.
  * Bobylev and Bajkova [2013] V. V. Bobylev and A. T. Bajkova.  Estimation of the pitch angle of the Galactic spiral pattern.  _Astronomy Letters_ , 39:759–764, November 2013.  doi: 10.1134/S1063773713110017.
  * Bronfman et al. [1988] L. Bronfman, R. S. Cohen, H. Alvarez, J. May, and P. Thaddeus.  A CO survey of the southern Milky Way - The mean radial distribution of molecular clouds within the solar circle.  _The Astrophysical Journal_ , 324:248, jan 1988.  ISSN 0004-637X.  doi: 10.1086/165892.  URL http://adsabs.harvard.edu/abs/1988ApJ...324..248B.
  * Bunn et al. [2003] E. F. Bunn, M. Zaldarriaga, M. Tegmark, and A. de Oliveira-Costa.  E/B decomposition of finite pixelized CMB maps.  _Phys. Rev. D_ , 67(2):023501, January 2003.  doi: 10.1103/PhysRevD.67.023501.
  * CORE Collaboration et al. [2016] CORE Collaboration, F. Finelli, M. Bucher, A. Achúcarro, M. Ballardini, N. Bartolo, D. Baumann, S. Clesse, J. Errard, W. Handley, M. Hindmarsh, K. Kiiveri, M. Kunz, A. Lasenby, M. Liguori, D. Paoletti, C. Ringeval, J. Väliviita, B. van Tent, V. Vennin, F. Arroja, M. Ashdown, A. J. Banday, R. Banerji, J. Baselmans, J. G. Bartlett, P. de Bernardis, M. Bersanelli, A. Bonaldi, J. Borril, F. R. Bouchet, F. Boulanger, T. Brinckmann, Z.-Y. Cai, M. Calvo, A. Challinor, J. Chluba, G. D’Amico, J. Delabrouille, J. María Diego, G. De Zotti, V. Desjacques, E. Di Valentino, S. Feeney, J. R. Fergusson, S. Ferraro, F. Forastieri, S. Galli, J. García-Bellido, R. T. Génova-Santos, M. Gerbino, J. González-Nuevo, S. Grandis, J. Greenslade, S. Hagstotz, S. Hanany, D. K. Hazra, C. Hernández-Monteagudo, E. Hivon, B. Hu, E. D. Kovetz, H. Kurki-Suonio, M. Lattanzi, J. Lesgourgues, J. Lizarraga, M. López-Caniego, G. Luzzi, B. Maffei, C. J. A. P. Martins, E. Martínez-González, D. McCarthy, S. Matarrese, A. Melchiorri, J.-B. Melin, A. Monfardini, P. Natoli, M. Negrello, F. Oppizzi, E. Pajer, S. P. Patil, M. Piat, G. Pisano, V. Poulin, A. Ravenni, M. Remazeilles, A. Renzi, D. Roest, L. Salvati, A. Tartari, G. Tasinato, J. Torrado, N. Trappe, M. Tucci, J. Urrestilla, P. Vielva, and R. Van de Weygaert.  Exploring Cosmic Origins with CORE: Inflation.  _ArXiv e-prints_ , December 2016.
  * Crutcher [2012] R. M. Crutcher.  Magnetic Fields in Molecular Clouds.  _ARA &A_, 50:29–63, September 2012.  doi: 10.1146/annurev-astro-081811-125514.
  * Dame et al. [2001] T. M. Dame, D. Hartmann, and P. Thaddeus.  The Milky Way in Molecular Clouds: A New Complete CO Survey.  _ApJ_ , 547:792–813, February 2001.  doi: 10.1086/318388.
  * Davis et al. [2012] B. L. Davis, J. C. Berrier, D. W. Shields, J. Kennefick, D. Kennefick, M. S. Seigar, C. H. S. Lacy, and I. Puerari.  Measurement of Galactic Logarithmic Spiral Arm Pitch Angle Using Two-dimensional Fast Fourier Transform Decomposition.  _ApJS_ , 199:33, April 2012.  doi: 10.1088/0067-0049/199/2/33.
  * Delabrouille et al. [2013] J. Delabrouille, M. Betoule, J.-B. Melin, M.-A. Miville-Deschênes, J. Gonzalez-Nuevo, M. Le Jeune, G. Castex, G. de Zotti, S. Basak, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, J.-P. Bernard, F. R. Bouchet, D. L. Clements, A. da Silva, C. Dickinson, F. Dodu, K. Dolag, F. Elsner, L. Fauvet, G. Faÿ, G. Giardino, S. Leach, J. Lesgourgues, M. Liguori, J. F. Macías-Pérez, M. Massardi, S. Matarrese, P. Mazzotta, L. Montier, S. Mottet, R. Paladini, B. Partridge, R. Piffaretti, G. Prezeau, S. Prunet, S. Ricciardi, M. Roman, B. Schaefer, and L. Toffolatti.  The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths.  _A &A_, 553:A96, May 2013.  doi: 10.1051/0004-6361/201220019.
  * Ellsworth-Bowers et al. [2015] T. P. Ellsworth-Bowers, E. Rosolowsky, J. Glenn, A. Ginsburg, N. J. Evans, II, C. Battersby, Y. L. Shirley, and B. Svoboda.  The Bolocam Galactic Plane Survey. XII. Distance Catalog Expansion Using Kinematic Isolation of Dense Molecular Cloud Structures with 13CO(1-0).  _ApJ_ , 799:29, January 2015.  doi: 10.1088/0004-637X/799/1/29.
  * Gold et al. [2011] B. Gold, N. Odegard, J. L. Weiland, R. S. Hill, A. Kogut, C. L. Bennett, G. Hinshaw, X. Chen, J. Dunkley, M. Halpern, N. Jarosik, E. Komatsu, D. Larson, M. Limon, S. S. Meyer, M. R. Nolta, L. Page, K. M. Smith, D. N. Spergel, G. S. Tucker, E. Wollack, and E. L. Wright.  Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Galactic Foreground Emission.  _ApJS_ , 192:15, February 2011.  doi: 10.1088/0067-0049/192/2/15.
  * Goldreich and Kylafis [1981] P. Goldreich and N. D. Kylafis.  On mapping the magnetic field direction in molecular clouds by polarization measurements.  _The Astrophysical Journal_ , 243:L75, 1981.  ISSN 0004-637X.  doi: 10.1086/183446.  URL http://adsabs.harvard.edu/abs/1981ApJ...243L..75G.
  * Górski et al. [2005] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann.  HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere.  _ApJ_ , 622:759–771, April 2005.  doi: 10.1086/427976.
  * Grain et al. [2009] J. Grain, M. Tristram, and R. Stompor.  Polarized CMB power spectrum estimation using the pure pseudo-cross-spectrum approach.  _Phys. Rev. D_ , 79(12):123515, June 2009.  doi: 10.1103/PhysRevD.79.123515.
  * Greaves et al. [1999] J. S. Greaves, W. S. Holland, P. Friberg, and W. R. F. Dent.  Polarized CO Emission from Molecular Clouds.  _ApJ_ , 512:L139–L142, February 1999.  doi: 10.1086/311888.
  * Heyer and Dame [2015] Mark Heyer and T.M. Dame.  Molecular Clouds in the Milky Way.  _Annual Review of Astronomy and Astrophysics_ , 53(1):583–629, aug 2015.  ISSN 0066-4146.  doi: 10.1146/annurev-astro-082214-122324.
  * Hivon et al. [2002] E. Hivon, K. M. Górski, C. B. Netterfield, B. P. Crill, S. Prunet, and F. Hansen.  MASTER of the Cosmic Microwave Background Anisotropy Power Spectrum: A Fast Method for Statistical Analysis of Large and Complex Cosmic Microwave Background Data Sets.  _ApJ_ , 567:2–17, March 2002.  doi: 10.1086/338126.
  * Hu and White [1997] Wayne Hu and Martin J. White.  A CMB polarization primer.  _New Astron._ , 2:323, 1997.  doi: 10.1016/S1384-1076(97)00022-5.
  * Jackson et al. [2006] J. M. Jackson, J. M. Rathborne, R. Y. Shah, R. Simon, T. M. Bania, D. P. Clemens, E. T. Chambers, A. M. Johnson, M. Dormody, R. Lavoie, and M. H. Heyer.  The Boston University-Five College Radio Astronomy Observatory Galactic Ring Survey.  _ApJS_ , 163:145–159, March 2006.  doi: 10.1086/500091.
  * Krachmalnicoff et al. [2015] N. Krachmalnicoff, C. Baccigalupi, J. Aumont, M. Bersanelli, and A. Mennella.  Characterization of foreground emission at degree angular scale for CMB B-modes observations. Thermal Dust and Synchrotron signal from Planck and WMAP data.  _Astronomy & Astrophysics_, page 10, 2015.  URL http://arxiv.org/abs/1511.00532.
  * Lewis et al. [2001] Antony Lewis, Anthony Challinor, and Neil Turok.  Analysis of cmb polarization on an incomplete sky.  _Phys. Rev. D_ , 65:023505, Dec 2001.  doi: 10.1103/PhysRevD.65.023505.  URL http://link.aps.org/doi/10.1103/PhysRevD.65.023505.
  * Matsumura et al. [2014] T. Matsumura, Y. Akiba, J. Borrill, Y. Chinone, M. Dobbs, H. Fuke, A. Ghribi, M. Hasegawa, K. Hattori, M. Hattori, M. Hazumi, W. Holzapfel, Y. Inoue, K. Ishidoshiro, H. Ishino, H. Ishitsuka, K. Karatsu, N. Katayama, I. Kawano, A. Kibayashi, Y. Kibe, K. Kimura, N. Kimura, K. Koga, M. Kozu, E. Komatsu, A. Lee, H. Matsuhara, S. Mima, K. Mitsuda, K. Mizukami, H. Morii, T. Morishima, S. Murayama, M. Nagai, R. Nagata, S. Nakamura, M. Naruse, K. Natsume, T. Nishibori, H. Nishino, A. Noda, T. Noguchi, H. Ogawa, S. Oguri, I. Ohta, C. Otani, P. Richards, S. Sakai, N. Sato, Y. Sato, Y. Sekimoto, A. Shimizu, K. Shinozaki, H. Sugita, T. Suzuki, A. Suzuki, O. Tajima, S. Takada, S. Takakura, Y. Takei, T. Tomaru, Y. Uzawa, T. Wada, H. Watanabe, M. Yoshida, N. Yamasaki, T. Yoshida, and K. Yotsumoto.  Mission Design of LiteBIRD.  _Journal of Low Temperature Physics_ , 176:733–740, September 2014.  doi: 10.1007/s10909-013-0996-1.
  * Mizuno and Fukui [2004] A. Mizuno and Y. Fukui.  Physical properties of molecular clouds as revealed by NANTEN CO survey: from the galactic center to the galactic warp.  In D. Clemens, R. Shah, and T. Brainerd, editors, _Milky Way Surveys: The Structure and Evolution of our Galaxy_ , volume 317 of _Astronomical Society of the Pacific Conference Series_ , page 59, December 2004.
  * Page et al. [2007] L. Page, G. Hinshaw, E. Komatsu, M. R. Nolta, D. N. Spergel, C. L. Bennett, C. Barnes, R. Bean, O. Doré, J. Dunkley, M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, H. V. Peiris, G. S. Tucker, L. Verde, J. L. Weiland, E. Wollack, and E. L. Wright.  Three-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Polarization Analysis.  _ApJS_ , 170:335–376, June 2007.  doi: 10.1086/513699.
  * Planck Collaboration et al. [2014a] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. I. R. Alves, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, and et al.  Planck 2013 results. XIII. Galactic CO emission.  _A &A_, 571:A13, November 2014a.  doi: 10.1051/0004-6361/201321553.
  * Planck Collaboration et al. [2014b] Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. J. Banday, and et al.  Planck 2013 results. IX. HFI spectral response.  _A &A_, 571:A9, November 2014b.  doi: 10.1051/0004-6361/201321531.
  * Planck Collaboration et al. [2016a] Planck Collaboration, R. Adam, P. A. R. Ade, N. Aghanim, M. I. R. Alves, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, and et al.  Planck 2015 results. X. Diffuse component separation: Foreground maps.  _A &A_, 594:A10, September 2016a.  doi: 10.1051/0004-6361/201525967.
  * Planck Collaboration et al. [2016b] Planck Collaboration, R. Adam, P. A. R. Ade, N. Aghanim, M. Arnaud, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al.  Planck intermediate results. XXX. The angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes.  _A &A_, 586:A133, February 2016b.  doi: 10.1051/0004-6361/201425034.
  * Planck Collaboration et al. [2016c] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. I. R. Alves, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, and et al.  Planck 2015 results. XXV. Diffuse low-frequency Galactic foregrounds.  _A &A_, 594:A25, September 2016c.  doi: 10.1051/0004-6361/201526803.
  * Planck Collaboration et al. [2017] Planck Collaboration, N. Aghanim, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, K. Benabed, J.-P. Bernard, M. Bersanelli, P. Bielewicz, A. Bonaldi, L. Bonavera, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, A. Bracco, C. Burigana, E. Calabrese, J.-F. Cardoso, H. C. Chiang, L. P. L. Colombo, C. Combet, B. Comis, B. P. Crill, A. Curto, F. Cuttaia, R. J. Davis, P. de Bernardis, A. de Rosa, G. de Zotti, J. Delabrouille, J.-M. Delouis, E. Di Valentino, C. Dickinson, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, E. Falgarone, Y. Fantaye, F. Finelli, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, M. Giard, J. González-Nuevo, K. M. Górski, A. Gregorio, A. Gruppuso, J. E. Gudmundsson, F. K. Hansen, G. Helou, D. Herranz, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, E. Keihänen, R. Keskitalo, T. S. Kisner, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, A. Lähteenmäki, J.-M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, F. Levrier, M. Liguori, P. B. Lilje, M. López-Caniego, P. M. Lubin, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, M. Maris, P. G. Martin, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, A. Melchiorri, A. Mennella, M. Migliaccio, S. Mitra, M.-A. Miville-Deschênes, D. Molinari, A. Moneti, L. Montier, G. Morgante, A. Moss, P. Naselsky, H. U. Nørgaard-Nielsen, C. A. Oxborrow, L. Pagano, D. Paoletti, B. Partridge, L. Patrizii, O. Perdereau, L. Perotto, V. Pettorino, F. Piacentini, S. Plaszczynski, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, M. Rossetti, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, C. Sirignano, G. Sirri, L. Stanco, A.-S. Suur-Uski, J. A. Tauber, M. Tenti, L. Toffolatti, M. Tomasi, M. Tristram, T. Trombetti, J. Valiviita, F. Vansyngel, F. Van Tent, P. Vielva, B. D. Wandelt, I. K. Wehus, A. Zacchei, and A. Zonca.  Planck intermediate results. L. Evidence of spatial variation of the polarized thermal dust spectral energy distribution and implications for CMB B-mode analysis.  _A &A_, 599:A51, February 2017.  doi: 10.1051/0004-6361/201629164.
  * Regan et al. [2002] M. W. Regan, K. Sheth, P. J. Teuben, and S. N. Vogel.  Inner Molecular Rings in Barred Galaxies: BIMA Survey of Nearby Galaxies CO Observations.  _ApJ_ , 574:126–133, July 2002.  doi: 10.1086/340793.
  * Roman-Duval et al. [2010] J. Roman-Duval, J. M. Jackson, M. Heyer, J. Rathborne, and R. Simon.  Physical Properties and Galactic Distribution of Molecular Clouds Identified in the Galactic Ring Survey.  _ApJ_ , 723:492–507, November 2010.  doi: 10.1088/0004-637X/723/1/492.
  * Roman-Duval et al. [2016] J. Roman-Duval, M. Heyer, C. M. Brunt, P. Clark, R. Klessen, and R. Shetty.  Distribution and Mass of Diffuse and Dense CO Gas in the Milky Way.  _ApJ_ , 818:144, February 2016.  doi: 10.3847/0004-637X/818/2/144.
  * Seljak and Zaldarriaga [1997] U. Seljak and M. Zaldarriaga.  Signature of Gravity Waves in the Polarization of the Microwave Background.  _Physical Review Letters_ , 78:2054–2057, March 1997.  doi: 10.1103/PhysRevLett.78.2054.
  * Smith and Zaldarriaga [2007] Kendrick M. Smith and Matias Zaldarriaga.  General solution to the $e\mathrm{\text{-}}b$ mixing problem.  _Phys. Rev. D_ , 76:043001, Aug 2007.  doi: 10.1103/PhysRevD.76.043001.  URL http://link.aps.org/doi/10.1103/PhysRevD.76.043001.
  * Tassis and Pavlidou [2015] K. Tassis and V. Pavlidou.  Searching for inflationary B modes: can dust emission properties be extrapolated from 350 GHz to 150 GHz?  _MNRAS_ , 451:L90–L94, July 2015.  doi: 10.1093/mnrasl/slv077.
  * Vallée [2014] J. P. Vallée.  Catalog of Observed Tangents to the Spiral Arms in the Milky Way Galaxy.  _ApJS_ , 215:1, November 2014.  doi: 10.1088/0067-0049/215/1/1.
  * Wolfire et al. [2003] M. G. Wolfire, C. F. McKee, D. Hollenbach, and A. G. G. M. Tielens.  Neutral Atomic Phases of the Interstellar Medium in the Galaxy.  _ApJ_ , 587:278–311, April 2003.  doi: 10.1086/368016.

##  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.

Figure 13:  Angular power spectra of Planck CO $1-0$ line (blue) for _Type 1_ (left) and _Type 2_ (right) maps together with the results of the MCMole3D best-fit model adopting an Axisymmetric geometry (magenta). |  $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
Table 3:  Best-fit parameters for the Axisymmetric MCMole3D model.

Figure 14:  Left: CO 1-0 power spectrum at high Galactic latitudes of the
Axisymmetric MCMole3D model (dotted green), for the best-fit parameters
reported in Table 3. Thick orange solid line shows the power spectrum for the
Axisymmetric geometry with the HGL-cut applied. The Planck  _Type 1_ CO power
spectrum before and after noise bias subtraction is shown with the blue solid
line and dots respectively. The Planck noise bias is shown with the dashed red
line. Right: B-mode power spectra of polarized CO emission lines at high
Galactic latitudes estimated using the best-fit parameters of the Axisymmetric
MCMole3D model, see Figure 12 for a comparison with the LogSpiral geometry.