A 3D model for CO molecular line emission as a potential CMB polarization contaminant
Giuseppe Puglisi Giulio Fabbian Carlo Baccigalupi



firstpage–lastpage
Antiferromagnetic textures and dynamics on the surface of a heavy metal
    Ricardo Zarzuela and Yaroslav Tserkovnyak
    December 30, 2023
=======================================================================


 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  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 ≲ 
0.02. The Monte-Carlo MOlecular Line Emission ()  
package, which implements this model, is being made publicly available.





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





§ 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 (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-H_2
conversion factor, X_CO. A commonly accepted value for X_CO is 2× 10^20 molecules· cm^-2 (K  km
 s^-1)^-1, although this could vary with position in  the
Galactic plane, particularly in the outer Galaxy <cit.>.



The most intense  CO rotational transition lines are the J = 1→0, 
2→1,  3→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 ^12CO emission are complemented by
measurements of  ^13CO 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 
H_2 abundance <cit.>.
However, there is growing evidence that  ^13CO regions could be
associated with  colder and denser environments, whereas ^12CO
emission originates from a diffuse component of molecular gas
<cit.>.

 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 <cit.>  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 <cit.>. Usually these surveys
have observed a strip of |b|≲ 5 around the Galactic plane. At
higher Galactic latitudes (|b|>30), 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  ≈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)
<cit.>. 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  frequency channels. These single frequency maps
have been processed with a dedicated foreground cleaning procedure so as to
isolate this emission. The  maps were found to be broadly consistent with
the data from other CO surveys <cit.>, although
they might still be affected by residual astrophysical emissions and
instrumental systematics.  In <ref>, we show the so called
Type 1  map of the CO J:1-0 line  <cit.>[http://pla.esac.esa.int/pla] which will be
used in the following.


Many current and future CMB polarization experiments[For 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
<cit.>.  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 <cit.> and low frequencies
<cit.> 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
<cit.>. 


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
<cit.> 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.  <cit.> 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
^13CO. The degree of polarization ranged from 0.5
 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 <cit.>. In <ref> we present the assumptions made for building the model and the simulation pipeline for its implementation. In
<ref> we describe the methodology for calibrating the CO
simulations to match  observations. 

In <ref> we show how the parameters describing molecular cloud
distribution shape the angular power spectrum of CO emission. Finally, in
<ref> we forecast the expected level of polarized CO
contaminations for  the  B-modes at high Galactic latitudes using our
calibrated simulation of <ref> to infer statistically the
emission at high Galactic latitude, where current observations are less
reliable.



§ 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 [https://github.com/giuspugl/MCMole3D]
which is being made publicly available, and we present  details  of it in this
Section[In the following we will refer to this model as the model for the  sake of clarity.]. The model builds on and extends the method
proposed by <cit.> who conducted a series of analyses
distributing statistically a relative large number of molecular cloud objects
according to the axisymmetric distribution of H_2 observed in the
Galaxy <cit.>.



 §.§ 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 <ref>(b) and (a) respectively.  
The first is a simplified model and is parametrized by R_ring,and σ_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 ∼ -13  following the latest measurements of <cit.> and fixed the arm half-width to be 340 pc <cit.>.


<cit.> 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 <cit.> 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 <cit.>. 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)∝σ_z(R)=σ_z,0cosh(R/h_R ),

where σ_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 σ_z through the usual relation z_1/2=√(2 ln 2)σ_z.
The final vertical profile is then:
 
    z(R) = 1/√(2 π)σ_z(R) exp[- (z/√(2)σ_z(R) )^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 <cit.>:  

    ϵ_0(R)=ϵ_c exp(R/R_em),

where ϵ_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 ξ(L) defined by their typical size scale, L_0, the range of sizes [L_min, L_max] and two power-laws with spectral indices <cit.>

    ξ(L)=dn/dL∝{[         L^0.8 L_min<L< L_0,;        L^-α_L L_0<L< L_max, ].

with α_L= 3.3, 3.9 for clouds inside or outside the solar circle respectively. From the cloud size function ξ(L) we derive the corresponding probability 𝒫(L) of having clouds with sizes smaller than L:

    𝒫(<L)=∫_L_min ^L d L'ξ(L') .
  
The probability functions for different choices of the spectral index α_L are shown in <ref>. 
We then inverted <ref> to get the cloud size associated with a given probability L(p). The cloud sizes are drawn from a uniform distribution in [ 0,1]. The histograms of the sizes generated following this probability function are shown in the top panel of <ref> 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  boundaries[For the Gaussian profile, we set σ in order to have the cloud boundaries at 6σ, 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 <ref>). 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. 





 §.§ 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  Θ^CO that can be set by the user:

    Θ^CO= {    N_clouds, ϵ_c, R_em,  R_ring, σ_ring,
       σ_z,0, h_R,   L_min,L_max, L_0}.

The values chosen for our analysis are listed in <ref>.

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 <cit.>). The product of each simulation is a map, similar to the one in <ref>, in the Hierarchical Equal Area Latitute Pixelization  (HEALPIX, <cit.> ) [http://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  maps described in <ref>, we convolve the simulated maps to the beam resolution of the 100 GHz channel (∼ 10 arcmin).
 
The procedure implemented for each realization is the following: 


  * assign the (R_gal,ϕ,z) Galacto-centric positions. In particular:


  * R_gal is extracted from a Gaussian distribution defined by  the R_ring and σ_ring  parameters. However, the σ_ring is large enough to give non-zero  probability at R_gal≤ 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 <cit.>  for a similar approach). 


  * the z-coordinate is drawn randomly from the distribution in <ref>.
 

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

ϕ(R)= A logR +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;
 

  * assign cloud sizes given  the probability function 𝒫(L) (<ref>);

  * assign emissivities to each cloud from the emissivity profile (see <ref>);

  * convert (R_gal,ϕ,z) positions into the heliocentric  coordinate frame (ℓ,b,d_⊙). 


In <ref> 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.





 §.§ Simulation results

In <ref> we show two typical realizations of maps of CO emission for the  and  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 𝒞_ℓ corresponding to different distributions of CO emission in <ref>. In the ones shown subsequently the spectra are  𝒟_ℓ encoding a normalization factor 𝒟_ℓ =ℓ (ℓ +1)𝒞_ℓ/2π. 


We can observe two main features in the morphology of the power spectrum: a bump around ℓ∼ 100 and a tail at higher ℓ.  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 (∼ 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 ℓ≳ 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 σ_ring parameters modify the power spectrum in two different ways. For a given typical size, if the width of the molecular ring zone σ_ring increases, the peak around ℓ∼ 100 shifts towards lower multipoles, i.e. larger angular scales, and its amplitude increases proportionally to σ_ring, see  for instance the bottom right panel in <ref>.  
This can be interpreted as corresponding to the fact that the larger is σ_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  <ref>) the tail at small angular scales moves downwards and flattens as L_0 increases.   Vice versa, if we keep L_0 constant (<ref> bottom right panel), all of the tails have the same amplitude and an ℓ^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 <ref> 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  and  geometries. For this run we choose  the so-called best fit values for the L_0 and σ_ring parameters discussed later in <ref>. 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  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 <ref>. In the  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.



§ COMPARISON WITH PLANCK DATA



 §.§ Dataset

The  collaboration released three different kinds of CO molecular line emission maps, described in <cit.>. 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  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  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  spectral response at 100 GHz is roughly 3 GHz,  which corresponds to ∼ 8000 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   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 <cit.>.



 §.§ 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|≤30) and high Galactic latitude (|b|>30). 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 <ref> we show the angular power spectra  of the  first three CO rotational line maps observed by   as well as the expected noise level at both high and low Galactic latitudes computed using a pure power spectrum estimator  <cit.>.  This is a pseudo power spectrum method <cit.> which corrects the so called E-to-B-modes leakage in the polarization field that arises in the presence of incomplete sky coverage <cit.>. 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 <ref>. We estimated the noise as the mean of 100 MC Gaussian simulations based on the the diagonal pixel-pixel error covariance included in the  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 ℓ≈ 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 <cit.>

    ΔC̃_ℓ = √(2/ν)(C_ℓ + N_ℓ)

where ν is the number of degrees of freedom taking into account the finite number of modes going into the power spectrum calculation in each ℓ mode and the effective sky coverage. N_ℓ represents  the noise power spectrum and the C_ℓ 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_ℓ + N_ℓ = C̃_ℓ. 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 <ref>, all of the power spectra of CO emission at low Galactic latitudes have a broad peak around  the multipole 100÷ 300, i.e. at the ≈1 angular scale. The signal power starts decreasing up to ℓ∼ 600 and then grows again at higher ℓ 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 ∼ 0.5 deg area in the sky. This corresponds to a correlated scale in the power spectrum at an ℓ of the order of a few hundred but the detail of this scale depends on the width of the molecular ring zone.





 §.§ 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  observed data, we convolved the  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  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 <cit.>, we therefore computed the integrated flux of the emission along the two Galactic latitudes and longitudes (l, b) defined as

    I^X(l)=∫ db I^X(l, b), 
       I_tot^X=∫ dl db I^X(l,b),

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=I_tot^observ/I_tot^model.

We estimated the integrals in <ref> and <ref> by considering a narrow strip of Galactic latitudes within [-2,2] 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|≲ 2 deg.


In <ref> we show the comparison between  I^observ(l) and the I^model(l) as defined in <ref> computed as the mean of 100 MC realizations of galaxies populated by molecular clouds  for both the  and  models as well as their typical standard deviation. In particular, we chose for these simulations the default parameters in <ref>. 


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 [  -90 ,90 ] deg[We 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 <ref>). 


In addition, we note that the bump in the profile at l ≃ 60 - 70 deg, where we see a lack of power in both the  and  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× 10^6 M_⊙, 1.4 kpc distant from the Sun. 
Given the assumptions made in <ref>, these  large and closer clouds are not easy to simulate with  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  model than for the  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 ∼ -140, -80, 120 deg. We will perform a more detailed comparison of the two geometries in the following section and in <ref>.



 §.§ Constraining the  model with Planck data



After comparing the CO profile emission we checked whether the  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  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 <ref>). For sake of completeness we reported the results of the same analysis adopting an  geometry in <ref>. 


We left  the typical cloud size L_0 and σ_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 <ref>, where L_0 and σ_ring are allowed to vary. Looking at the series of examples reported in <ref> we can see that suitable parameter ranges which yield power spectra close to the  observations are L_0=10 ÷30 pc and σ_ring=2 ÷ 3 kpc. It is interesting to note that these are in agreement with estimates available in the literature (see e.g. <cit.>).


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 . 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 ℱ(ℓ; σ_ring, L_0) which for a given set of values (σ_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  model by performing a χ^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  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  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  frequency band.

Since the theoretical model has been estimated from Monte Carlo simulations, we added linearly to the sample variance error of the  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 <ref> on the  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 ℓ≤400 in order  to avoid the regions that display an unusual bump at scales of around ℓ≈ 500 that is not captured by any realization of our model (see next section). The best-fit parameters are reported in <ref>

    L_0 =  14.50 ± 0.58 pc, 
       σ_ring= 2.76 ± 0.19 Kpc, 
       A_CO = 0.69 ± 0.06 .

The values are within the ranges expected from the literature. As can be seen in  <ref> the power spectrum corresponding to the model with the best fit parameters, describes the  data reasonably well. The minimum χ^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 σ_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 <ref>(d).


Finally, we note that A_CO≲ 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  channel or the finite width of the CO line. Additional sources of signal overestimation could be residual contamination of ^13CO 1-0 line or thermal dust in the map or variations of the emissivity profile in <ref>.



 §.§ Consistency checks on other maps

The  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 <cit.>. 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 <cit.>.


We repeated the procedure outlined in <ref> and <ref> using the Type 2 1-0 map. The values of the best fit parameters are summarized in <ref> and we show in <ref> 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 <cit.>.
Such differences are mainly related to a combination of bandpass uncertainties in the  observations and presence of a mixture of ^12 CO and ^13 CO  (emitted at 110 GHz) lines for the Type 1 maps. 
While σ_ring is consistent between the two maps, the Type 2 L_0 parameters are in slight tension at 2.7σ level. The overall correlation of the parameters is increased and the overall agreement between data and the  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  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  bandpass induced by the motion of the clouds themselves <cit.>. 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  maps display an oscillatory behaviour at a scale of ℓ≥ 400 with a clear peak at around ℓ≈ 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.



 §.§ Comparison with data at high Galactic latitudes


In <ref>, we compare the  CO 1-0 power spectrum at high Galactic latitudes with the average power spectrum of 100 MC realizations of the  model for the same region of the sky. We assumed for these runs the best-fit values of the L_0, σ_ring parameters reported in <ref> and a  distribution. Because the  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 <ref>,  some discrepancy arises when comparing the power spectrum expected from the simulation of   at high Galactic latitudes with the noise debiased  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 σ_z,0 and h_R (see <ref>) which gives  a higher number of clouds close to the observer at high latitude. 
However, we also point out that, as explained in <ref>, the error bars in <ref> might be underestimated especially at the largest angular scales where we are signal dominated. Therefore, discrepancies of  order ≈ 3σ 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  data. All of these clouds have, in fact, a flux exceeding the  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  simulated maps obtained after the application of the  HGL cut is shown in <ref>.
 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. ℓ≲ 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.





§ POLARIZATION FORECASTS


As noted in <ref>, 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.



 §.§ Data-based order of magnitude estimate

Starting from the measured  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

    𝒞_ℓ  ^high, CO= 𝒞_ℓ  ^Galvar(high)/var(Gal) .


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  noise at high Galactic latitudes is diffuse CO emission. The variance of the  CO map is 0.3  K^2 (km  s^-1)^2, at |b|>30 deg, while for |b|<30 deg we get a variance of 193.5 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 𝒞_ℓ^CO high  into its B-mode counterpart as 𝒞_ℓ ^CO high, EE=𝒞_ℓ ^CO high, BB= 𝒞_ℓ^high, COp_CO^2/2. We then apply the conversion factors of <cit.> to convert the CO power spectrum into thermodynamic units (from K_RJ km s^-1 to μK). We can compare 𝒞_ℓ ^CO high, BB to the amplitude of  equivalent cosmological CMB inflationary B-modes with tensor-to-scalar ratio r=1 at ℓ=80.  In terms of 𝒟^BB_ℓ, this  is equal to ∼ 6.67 × 10^-2μ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. 



 §.§ Simulation estimate

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


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  <cit.> the Q and U Stokes parameters of each CO cloud can be related to the unpolarized emission as


    Q(n̂)^CO   =p_CO g_d(n̂)I(n̂)^CO  cos(2 ψ(n̂)),
    
    U(n̂)^CO   =p_CO g_d(n̂)I(n̂)^CO  sin(2 ψ(n̂))..
 
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 ψ 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_.  Following the findings of <cit.>, 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 <cit.>, we used the g_d and ψ templates for the Galactic dust emission available in the public release of the Planck Sky  Model suite[<http://www.apc.univ-paris7.fr/ delabrou/PSM/psm.html>] <cit.>. 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  observations <cit.> 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 <ref> 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 .

In <ref> 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 <ref> 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 <ref>, there is a significant dispersion compared to the results of the MC simulations at low Galactic latitudes (see <ref>). 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  geometry and changing the parameter σ_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 𝒟_ℓ∼ℓ^α, with α = -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 ℓ∼ 80, 𝒟_ℓ=(1.1 ± 0.8)× 10^-4μ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 ± 0.002, i.e. below the recent upper limit r<0.07 reported by the BICEP2 Collaboration <cit.> but higher than the r=0.001 target of upcoming experiments <cit.>. The contamination quickly becomes  sub-dominant on small angular scales (ℓ≈ 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. 







§ 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  is being made publicly available. 


We have compared the results of our simulations with  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 σ_ring play a key role in shaping the power spectrum; 

  * the choice of symmetries in the cloud distribution changes the profile of the integrated emission in the Galactic plane (<ref>) but not the power spectrum morphology;

  * our model is capable of reproducing fairly well the observations at low Galactic latitudes (see <ref>) and the power spectrum at high latitudes (<ref>).

We used our model to fit the  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  model with  data to forecast the unresolved CO contamination of the  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≲ r≲ 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  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 <cit.> package.


plainnat









§ BEST-FIT WITH  GEOMETRY

In this appendix we present the results of the analysis described in <ref> to constraint the CO distribution using the  model adopting an  geometry instead of the  one. Following the procedure of <ref> we construct a series of ℱ(ℓ; σ_ring, L_0) hyper-surfaces sampled on an ensemble of specific values of the L_0 and σ_ring  parameters within the same ranges reported in <ref>.


 In <ref> we show the results of the fit of the axisymmetric  model to the CO power spectrum of the  Type 1 and Type 2 CO maps in the Galactic plane. We summarize the best-fit values of these parameters in <ref>. As it can be seen from the results of the χ^2 test in <ref> the  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  geometry as a baseline choice for our forecast presented in <ref>. Nevertheless, we pushed the comparison between the two geometries in the high galactic latitude area for sake of completeness.
 
In <ref> we show the comparison between the  data for Type 1 maps and the  axisymmetric best-fit model after the application of the HGL cut described in the paper. The  model describes the data similarly to the  model at the larger scales. The difference in the signal amplitude is in fact less then 30% for angular scales ℓ≲ 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 ≈ 2 at ℓ≈ 1000. 

We finally performed a series of polarized simulations as in <ref> to access the level of contamination to the CMB B-modes power spectrum with the best-fit  model and found r_CO≲0.001.
 Moreover, the slope of the BB power spectrum in <ref>(b) is -2.2 similar to the one computed with the  geometry.

 Because the  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.











