© ESO, 2017

1. Introduction

The structural, dynamical and thermodynamical evolution of interstellar clouds is often probed with observations of the H i 21-cm line for the atomic gas (Kalberla et al. 2005, 2010; Peek et al. 2011) and of CO rotational emission lines for the molecular gas (Heyer & Dame 2015). However, measuring gas column densities in the dense media, both atomic and molecular, is hampered by the opacity of the interstellar medium (ISM) to these radiations and we have no efficient means to correct for line opacities in the absence of absorption data against background point sources. Complementary information can be gained from the sub-millimeter and infrared thermal emission of large dust grains mixed with the gas (Planck Collaboration XI 2014), and from γ rays produced in interactions of high-energy cosmic rays (CRs) with interstellar gas nucleons (Lebrun et al. 1982). Dust grains and CRs trace all chemical and thermodynamical forms of the gas to large depths into the clouds (Planck Collaboration Int. XXVIII 2015), but they trace only integrals of the gas column densities, N_H, along sight lines, and they bear no kinematic information. To use them, we rely on critical assumptions, namely a uniform dust-to-gas mass ratio and uniform grain emissivity across a cloud and, likewise, a uniform CR flux through a cloud. Those assumptions still need to be tested in a variety of clouds and across complex phase changes, especially on consideration that dust grains and their radiation properties vary from the diffuse ISM to molecular clouds (Stepnik et al. 2003; Martin et al. 2012; Roy et al. 2013a; Ysard et al. 2013; Flagey et al. 2009; Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014; Planck Collaboration Int. XXVIII 2015). One can also test the uniformity of the γ-ray emissivity spectrum of the gas in well-sampled nearby clouds to verify the smooth penetration of CRs with energies above ~ 1 GeV (Skilling & Strong 1976; Cesarsky & Volk 1978; Morlino & Gabici 2015).

The recent progress in H i, CO, dust, and γ-ray observations provides excellent opportunities to confront these observations at a resolution of a few parsecs inside nearby clouds to probe biases in the different tracers and to search for missing gas in the census. This comparison can shed light on the highly opaque atomic gas forming the cold neutral medium (CNM), where the H ibrightness temperature probes the excitation (spin) temperature of the gas rather than the column density (Murray et al. 2015, and references therein). One can also probe how the ¹²CO, J = 1 → 0, emission rapidly saturates at high column density and how it fades out in the diffuse molecular gas where CO lines are too weakly excited and/or CO molecules are too efficiently photodissociated to map the corresponding molecular hydrogen (van Dishoeck & Black 1988; Wolfire et al. 2010). Large H iopacities and the lack of CO emission both conspire to hide large gas masses at the H–H₂ interface. This dark neutral medium (DNM) can be jointly revealed by cosmic rays and dust (Grenier et al. 2005).

A more precise census of the gas column densities is key to establishing and quantifying variations in dust opacity (optical depth per gas nucleon) across cloud structures. Mapping those variations as a function of grain temperature, T, and opacity spectral index, β gives important diagnostics for grain evolution through mantle accretion, coagulation, and ice coating (Köhler et al. 2015).

With an approach similar to the multi-tracer analysis of the Chamaeleon complex (Planck Collaboration Int. XXVIII 2015), we present an analysis of the anticentre region, which exhibits more massive clouds, with higher levels of star-formation activity, so we can test the properties of the gas tracers to larger column densities and molecular fractions. The analysed region spans a few hundred parsecs in distance in the local ISM (see below). It encompasses the well-known Taurus, Auriga, and Perseus molecular clouds (Ungerechts & Thaddeus 1987; Lombardi et al. 2010; Schnee et al. 2008), the massive complex associated with the California nebula (Lada et al. 2009), and a nearby, diffuse, hardly molecular cloud, which we refer to as Cetus.

Our goal is to test the dust and γ-ray emission as tracers of the total gas against linear combinations of the emission line intensities recorded in the H iand CO observations, a new free-free emission map of the ionised gas, plus the DNM gas which escapes the radio and millimeter spectroscopic surveys at the atomic-molecular interface. We have separated the H iand CO clouds in space (i.e., direction, velocity, and distance) to study their individual contributions to the total dust and γ-ray signals. We have used the model results to derive scaling factors such as the CO-to-H₂ conversion factor, X_CO, in each cloud, and the average dust opacity in the different gas phases for each cloud. We have then used the multi-wavelength data to follow dust evolution with increasing N_H in several molecular clouds with masses spanning nearly two orders of magnitude. A companion paper will present the results on the transitions between gas phases, and on the mass fractions and contribution of CO-dark H₂ to the molecular phase in individual regions.

The paper is structured as follows: the data are presented in Sect. 2 and the models are described in Sect. 3. The model results and the detection of gas in the DNM, CO-saturated, and ionised phases are detailed in Sect. 4. The CR content of the clouds, the X_CO ratios, and the dust opacities are discussed in Sects. 5, 6, and 7, respectively. Additionally, we present the H iand CO component separation into nearby cloud complexes in Appendix A, the derivation of the free-free intensities in Appendix B, the dependence of the goodness of fit on the H ispin temperature in Appendix C, the best-fit coefficients of the dust and γ-ray models in Appendix D, and a table with historical γ-ray measurements of X_CO factors in E.1.

2. Data

The analysis region extends from 139° to 191° in Galactic longitude and from − 56° to − 3° in Galactic latitude. We have masked regions with large column densities from the background Galactic disc, at − 9°<b< − 3° for l> 183° and l< 152° and regions with small column densities at declinations below − 5°where the H imaps have lower angular resolution. All maps in Fig. 1 have been projected onto the same 0.̊125-spaced Cartesian grid as that of the CfA CO survey (Dame et al. 2001). We have checked that the use of the Cartesian projection, instead of equal-area bins, does not significantly bias the results. The largest differences are less than a fifth of the uncertainties found in the jackknife tests discussed in Sect. 4.2.

2.1. HI and CO data

For this analysis we have used the 4′-resolution GALFA-HI survey where available (Peek et al. 2011) and we have complemented it with the 10.8′-resolution EBHIS survey elsewhere (Winkel et al. 2016). As shown in Fig. 1, the GALFA survey covers 74% of the region. GALFA data-cubes¹ were re-sampled into the 0.̊125-spaced Cartesian grid. We have used the narrow-band cubes with their original velocity resolution of 0.18 km s^-1 in the local standard of rest (LSR). The EBHIS frequency sampling is coarser, with a velocity resolution of 1.44 km s^-1. In the analysed region covered by both survey we have found a tight correlation in column densities between the two: (1)In order to trace the molecular gas we have used ¹²CO J = 1 → 0 observations at 115 GHz from the moment-masked CfA CO survey of the Galactic plane (Dame et al. 2001; Dame & Thaddeus 2004). Most of the survey is based on a 0.̊125-spaced Cartesian grid except for the high-latitude clouds, at b ≲ − 50°, which have been interpolated from 0.̊25 to 0.̊125.

2.2. HI and CO velocity separation

We have decomposed the H iand CO spectra into individual lines and we have used this information to identify six nearby entities that are coherent in velocity, position (in Galactic coordinates), and distance, in addition to the emission from the background Galactic disc. The H ivelocity resolution in the region covered by the EBHIS survey is coarser, but still adequate to separate the main structures because of the limited confusion along sight lines at medium-to-high latitudes. The separation process is three-fold and is described by Planck Collaboration Int. XXVIII (2015) and in Appendix A. The choice of six main entities has been guided by the structure in (l, b, v) density and by distance information obtained from stellar reddening (Schlafly et al. 2014). Details on the (l, b, v) cuts are given in Table A.1. The main structures found in the local ISM, namely Cetus, Main Taurus, South Taurus, North Taurus, Perseus, and California, together with the background Galactic structures are depicted in N_HI column densities and in W_CO intensities in Fig. 1. The set of H imaps includes all of the H iemission observed in this region except for the high velocity clouds (Wakker et al. 2008). In order to investigate the effect of the unknown H ioptical depth, we have derived all the N_HI maps for a sample of uniform spin temperatures (150, 200, 300, 400, 500, 600, 700, and 800 K) and for the optically thin case.

Distance estimates to the local clouds have been compiled from the photometric measurements of PanSTARRS-1 and the detection of reddening fronts toward stellar groups (Schlafly et al. 2014). We have selected the stars along lines of sight intersecting each (and only one) cloud within its W_CO contour at 0.5 K km s^-1 and we have calculated the average distance in each sample. For the reddening fronts to the stars toward Cetus, Main Taurus, South Taurus, North Taurus, Perseus, and California we find average distances of 190 ± 30 pc, 140 ± 30 pc, 160 ± 10 pc, 190 ± 50 pc, 270 ± 20 pc, and 410 ± 20 pc, respectively. They compare reasonably with the geometry inferred from the modelling of the far-UV continuum emission (Lim et al. 2013). The North Taurus component is dominated by the Auriga cloud near the northern edge of our analysis region; it also includes MBM 12 molecular cloud from Magnani et al. (1985) catalogue. The South Taurus component includes MBM 6, MBM 12 and MBM 18.

2.3. Dust data

In order to trace the dust column density, we have used the optical depth, τ₃₅₃, obtained at 353 GHz by Planck Collaboration XI (2014). They have modelled the spectral energy distributions of the thermal emission of the large grains with modified black-body curves in specific intensity, I_ν = τ_ν0B_ν(T_d)(ν/ν₀)^β. The optical depth τ_ν0 at frequency ν₀, the dust colour temperature T_d, and the spectral index β were fitted to the data from Planck and the Infrared Astronomical Satellite (IRAS) across the sky at 30′ resolution. The fits were then repeated at 5′ resolution by fixing β as obtained in the first step. This procedure limited the noise impact on the T_d–β degeneracy in faint regions. The modelled intensities were checked to be consistent with the data at all frequencies (see Fig. 11 in Planck Collaboration XI 2014). The contamination from CO line emission in the 353 GHz filter band amounts to a few per cent of the total signal. It was not removed so as to avoid increasing the noise in directions away from CO clouds.

The optical depth, τ_ν, and opacity, σ_ν, of the dust at frequency ν follow the relations: (2)We have also calculated the total power radiated by the dust, also referred to as radiance: (3)Figure 2 shows the dust properties thus measured across the analysis region.

The corresponding analysis of the Chamaeleon region (Planck Collaboration Int. XXVIII 2015) compared different dust tracers, namely the optical depth and an estimate of the dust extinction, A_VQ, empirically corrected for starlight intensity. The latter distribution best followed the gas column densities in the Chamaeleon region, but we have checked that this is not true in the anticentre region in particular toward regions of dust temperatures exceeding 23 K, so we present the results of our analyses using only the τ₃₅₃ optical depth. Understanding why the empirical correction used for A_VQ is less efficient in the anticentre region will be investigated separately.

2.4. Ionised gas

Ionised gas is visible in Hα emission (Finkbeiner 2003) in the California nebula around l = 160.̊1 and b = − 12.̊4, in the G159.6-18.5 H iiregion, and along the Eridanus loop. The Califonia Nebula, alias NGC 1499 or S 220, is ionised by ξ Persei, an O7.5III star located at a distance of 380 ± 70 pc (van Leeuwen 2007) that has run away from the Perseus OB 2 association (Shinohara & Ishida 1976). The 1.̊2-diameter region G159.6-18.5 has been seen by IRAS and the Wide-Field Infrared Survey Explorer (WISE) (Andersson et al. 2000; Anderson et al. 2014). It is excited by another run-away O9.5V star, HD 278942, at a distance of 190 ± 40 pc (van Leeuwen 2007), which would place it in front rather than behind the Perseus cloud (Bally et al. 2008).

The composite Hα map varies in angular resolution from 6′ to 1° across the region of interest and Hα emission can be efficiently absorbed behind dense clouds (see for instance the dusty absorption lane across G159.6-18.5 in Fig. 8 of Bally et al. 2008). We have therefore preferred measures of the free-free intensity at mm wavelengths to trace the ionised gas.

The Planck LFI map at 70 GHz (release 2.01 data, Planck Collaboration I 2016) shows free-free emission in spatial coincidence with the Hα emission. The angular resolution of 14′ of LFI better suits our analysis than the 1°-resolution of the free-free map inferred at 22 GHz from the 9-year WMAP data (Bennett et al. 2013). To separate the free-free intensity in the 70 GHz data, we have successively removed the contributions from the Cosmological Microwave Background (CMB), from dust emission, and from point sources. The method we have used to extract the free-free emission is described in Appendix B.

The resulting free-free dominated map is displayed in Fig. 1. It closely resembles the Hα intensity map of the region (Finkbeiner 2003), but for a uniform resolution of 14′ and a reduced contrast, by a factor 2 to 8, between the intensities of G159.6-18.5 and NGC 1499. because the Perseus clouds are more transparent to mm waves.

The integral of the free-free intensity along a sight line through ionised gas with temperature T_e and electron and ion number densities, n_e and n_i, scales at 70 GHz as: (4)where the emission measure is expressed in cm^-6 pc (Karzas & Latter 1961; Carson 1988). If the observed intensity results from discrete and homogeneous nebulae with comparable uniform electron densities, , I_ff roughly scales with the ion column density, (5)so we have used the I_ff map as a template for the N_HII distribution in the gas model and we discuss the scaling ratio obtained in the fit between N_HII and I_ff in Sect. 4.4.

2.5. Gamma-ray data

We have used six years of the Pass 8 photon data from Fermi-LAT, the associated instrument response functions (IRFs, P8R2_CLEAN_V6), and the corresponding isotropic spectrum for the extragalactic and residual instrumental backgrounds (Atwood et al. 2013; Acero et al. 2015). We have applied tight rejection criteria (CLEAN class selection, photon arrival directions within <100° of the Earth zenith and in time intervals when the LAT rocking angle was inferior to 52°) in order to reduce the contamination by residual CRs and Earth atmospheric γ rays in the photon data (see Nolan et al. 2012, for details). We have sampled the IRFs, the exposure map, the γ-ray emissivity spectrum, q_LIS, of the local interstellar gas (Casandjian 2015), and the spectrum of the isotropic background in 14 energy bins, 0.2 dex in width and centred from 10^2.3 to 10^4.9 MeV.

To ensure photon statistics robust enough to follow details in the spatial distributions of the different interstellar components, we have analysed the data in broad and independent energy bands, bounded by 10^2.6, 10^2.8, 10^3.2, 10^3.6, and 10⁵ MeV. We have also analysed the data in the integrated 10^2.6–10⁵ MeV band. The low-energy threshold was set to take advantage of the best LAT point-spread function (PSF), which strongly degrades at low energy. The analysis of the data in two lower-energy bands (10^2.2–10^2.4 and 10^2.4–10^2.6 MeV) proved to be less reliable to separate the gas phases in compact clouds. To optimise the angular resolution, we have selected all detected photons above 1.6 GeV, but only those at lower energy that converted to pairs in the front section of the tracker (Atwood et al. 2009; Ackermann et al. 2012a). Given the broad ranges of the energy bands, we have not corrected the fluxes for the energy resolution, which is better than 14% at these energies (Atwood et al. 2013). To account for the spill-over of emission produced outside the analysis region, but reconstructed inside it, we have modelled point sources and interstellar contributions in a region 4° wider than the analysis region. This choice corresponds to the 99.5% containment radius of the PSF in the lowest energy band. For the H iemission outside the analysis region not covered by EBHIS and GALFA surveys we have used the Leiden/Argentine/Bonn (LAB) survey (Kalberla et al. 2005).

The positions and the flux spectra of the γ-ray sources in the field are provided by the Fermi-LAT Third Source Catalog (Acero et al. 2015). The observed γ-ray emission also includes a contribution from the large-scale Galactic inverse Compton (IC) emission emanating from the interactions of CR electrons with the ISRF. The GALPROP² parameter file 54-LRYusifovXCO4z6R30-Ts150-mag2 has been used to generate an energy-dependent template of the Galactic IC emission across the analysis region (Ackermann et al. 2012e).

3. Models and analyses

To first order, the dust and γ-ray emissions should both trace the total gas column density, N_H. In order to compare them in Fig. 3, we have convolved the dust optical depth with the LAT PSF on the one hand, and we have subtracted non-gaseous γ-ray emissions obtained by the fitting (Sect. 4) from the γ-ray data on the other hand. Figure 3 shows strong similarities in the spatial distributions of both tracers, but it also reveals differences in their dynamical range in several places.

To detect neutral gas unaccounted for in N_HI and W_CO, we have used the fact that it is permeated by both cosmic rays and dust grains. We have therefore extracted the significant γ-ray and dust residuals over the N_HI, W_CO, and N_HII expectations and we have used the spatial correlation between those residuals to infer the additional gas column densities (see Sect. 3.3). We have separated the residuals for two types of environments: in regions of weak or no CO intensity, below 7 K km s^-1, which correspond to the DNM at the H-H₂ transitions; and in regions toward large CO intensities, above 7 K km s^-1, to capture additional H₂ gas where ¹²CO emission saturates and rarer isotopologues such as ¹³CO or C¹⁸O should be used. We refer to this saturated-CO molecular component as “CO_sat”.

3.1. Dust model

In the case of a uniform dust-to-gas mass ratio and uniform mass emission coefficient of the grains, the dust optical depth linearly scales with the total N_H. We have therefore modelled τ₃₅₃(l,b) in each direction as a linear combination of the gaseous contributions from the different phases (H ii , H i, DNM, CO-bright, CO_sat), with free normalisations to be fitted to the data, as in Planck Collaboration Int. XXVIII (2015). We have added a free isotropic term, y_iso, to account for the residual noise and the uncertainty in the zero level of the dust data (Planck Collaboration XI 2014). The τ₃₅₃(l,b) model can be expressed as: (6)where N_HI,i(l,b), W_CO,i(l,b), and I_ff(l,b) respectively denote the N_HI, W_CO, and free-free maps of the clouds depicted in Fig. 1. and stand for the column densities in the DNM and CO_sat components deduced from the coupled analyses of the γ-ray and dust data (see Sect. 3.3).

The y model parameters have been estimated using a χ² minimization. We expect the model uncertainties to exceed the measurement errors in τ₃₅₃(l,b) because of potential variations in grain properties through the clouds and because of the limitations of the gas tracers (survey sensitivities, emission saturation, self-absorption, etc.). As we cannot precisely determine the model uncertainties, we have set them to a fractional value of the data and we have determined this fraction to be 19% by reaching a reduced χ² of unity. This fraction is larger than the 3% to 9% error in the measurement of τ₃₅₃ across this region (Planck Collaboration XI 2014).

3.2. Gamma-ray model

Earlier studies have indicated that the bulk of the Galactic CRs radiating at 0.4–100 GeV have diffusion lengths far exceeding typical cloud dimensions and that they permeate all the H i-bright, DNM, and CO-bright gas phases. The observed γ-ray emission can therefore be modelled, to first order, by a linear combination of the same gaseous components as in the dust model. We have assumed that the emissivity spectrum of the gas follows the average one obtained in the local ISM (q_LIS(E), Casandjian 2015), but we have left a free normalisation in each energy band to account for possible deviations in CR density and spectrum. The model includes other radiation components such as the Galactic IC radiation, I_IC(l,b,E), the isotropic intensity mentioned above, I_iso(E), and point sources with individual flux spectra S_j(E). We have verified that the soft emission from the Earth limb is not detected in the present energy range for the choice of maximum zenith angle. The soft and transient emission from Sun and Moon is not expected to be detected as the number of γ-ray photons they emit over 6 years is negligible compared to those of the ISM components in the energy range studied. The γ-ray intensity I(l,b,E), expressed in cm^-2 s^-1 sr^-1 MeV^-1, can thus be modelled as: (7)with the and maps extracted from the coupled dust and γ-ray analyses (see Sect. 3.3).

The input q_LIS spectrum was based on four years of LAT data and on the correlation between the γ radiation and the N_HI column densities derived from the LAB survey, for a spin temperature of 140 K, at latitudes between 7° and 70° (Casandjian 2015). The q_HI,i scale factors in the model can therefore compensate for differences in the H idata (calibration, angular resolution, spin temperature) and potentially for cloud-to-cloud variations in CR flux. Such differences will affect the normalizations equally in all energy bands whereas a change in CR penetration in a specific cloud will show as an energy-dependent correction. For each cloud, the average γ-ray emissivity spectrum per H atom in the atomic phase is estimated from the product of the q_LIS spectrum and the best-fit q_HI,i normalization. This emissivity can be used to estimate the gas mass present in the other DNM, CO, and CO_sat parts of the cloud if one assumes a uniform CR flux across the whole structure.

The model includes 126 points sources from the 3FGL catalogue (Acero et al. 2015) inside the analysis region. Their flux spectra, S_j(E), have been computed with the spectral characteristics provided in the catalogue. Their q_Sj(E) normalization in the model allows for possible changes due to the longer exposure and different photon reconstruction dataset used in the present analysis (six years of Pass 8 data instead of four years of Pass 7 reprocessed data for 3FGL) and due to the use of a different interstellar background for source detection in 3FGL. The sources have been fitted simultaneously. The sources present in the 4°-wide peripheral band around the analysis region have been merged into a single map, S_ext(l,b,E), and its normalization, q_Sext, has been left free in each energy band. The IC emission has been derived from a GALPROP model (Ackermann et al. 2012e) and isotropic emission has been determined over the whole sky with a different interstellar background model (Casandjian 2015), so we have left their scaling free in each energy band.

In order to compare with the LAT photon data in the different energy bands, we have multiplied by the exposure and convolved by the LAT PSF each model component. Figure 4 gives the photon yields in the total energy band obtained by the fitting (Sect. 4) from those components. It shows that the emission originating from the gas dominates over other types of emission and that the LAT angular resolution allows the spatial separation of the various clouds, and of the different gas phases within the clouds. The q model parameters have been fitted to the LAT photon data in each energy band using a binned maximum-likelihood with Poisson statistics.

3.3. Analysis iterations

An important aspect of our analysis is the iterative coupling of the dust and γ-ray models in order to extract the common DNM and CO_sat gas components which are present in both datasets, but for which we have no a priori templates. Both components show up as positive residuals over the expectations from the H i-bright, CO-bright, and free-free-bright gas components. To extract them, we have built maps of the positive residuals between the data (dust or γ rays) and the best-fit contributions from the N_HI, W_CO, free-free, and ancillary (other than gas) components. We have kept only the positive residuals above the noise (see below). We have separated the residuals into DNM and CO_sat maps according to the W_CO intensity in each direction (below and above 7 K km s^-1 for the DNM and CO_sat, respectively). The DNM and CO_sat templates estimated from the dust are provided to the γ-ray model ( and in Eq. (7)); conversely, the DNM and CO_sat column densities derived from the γ rays are provided to the dust model ( and in Eq. (6)). We have started the iterative process by fitting the dust optical depth with the H i, CO, free-free, and isotropic components to build the first DNM and CO_sat maps for the γ-ray model. We have then iterated between the γ-ray and dust models until a saturation in the log-likelihood of the fit to γ-ray data is reached (from the third to the fourth iteration).

The estimates of the q and y model coefficients and the DNM and CO_sat templates improve at each iteration since there is less and less need for the other components, in particular the H iand CO ones, to compensate for the missing gas. They still do at some level because the DNM and CO_sat templates provided by the γ rays or dust emission have limitations (e.g., dust emissivity variations, limited γ-ray sensitivity).

Care must be taken in the extraction of the positive residuals because of the Gaussian noise around zero. A simple cut at zero would induce an offset bias, so we have de-noised the residual maps using the multi-resolution support method implemented in the MR filter software (Starck & Pierre 1998). The stability of the iterative analysis and the results of the dust and γ-ray fits are discussed in the following sections.

4. Results

4.1. H I optical depth correction

The unknown level of H iopacity in the different maps induces systematic uncertainties on the N_HI column densities, therefore on the H icontributions to the models. The γ rays can help constrain the average level of H ioptical-depth correction applicable to the whole region by comparing the T_S-dependent contrast of the N_HI maps with the structure of the γ-ray flux emerging from the H igas. We have not tested different spin temperatures for each cloud complex. The maximum log-likelihood value of the γ-ray fit peaks for an H ispin temperature near 400 K (see Fig. C.1). In the following, we present the results obtained for this temperature, unless otherwise mentioned.

Spin temperatures are generally measured from pairs of H iemission and absorption spectra against bright radio sources (Heiles & Troland 2003; Mohan et al. 2004; Kanekar et al. 2011; Roy et al. 2013b; Murray et al. 2015). Line-of-sight harmonic means assume a single spin temperature along a sight line over the whole H iline or per velocity channel, as in the present work (monophasic T_S). The average value near 400 K that we find is consistent with the sparse distribution of previous monophasic spin measurements in this region (Mohan et al. 2004; Kanekar et al. 2011) and with the temperature span of 200–600 K found at high Galactic latitudes for the range of N_HI column densities dominating the present maps. But the monophasic assumption is known to bias the spin temperatures to higher values and to characterise the mixture of the cold and warm neutral mediums (CNM and WNM) along the line of sight, rather than the physical temperatures of individual structures (Murray et al. 2015). Monophasic temperatures tend to decrease with increasing N_HI column densities because of the increasing proportion of CNM along the line of sight (Kanekar et al. 2011).

4.2. Best fits and jackknife tests

The best dust and γ-ray fits that could be achieved with the models described by Eqs. (6)and (7)include all the free-free, H i, DNM, CO, and CO_sat templates for the gaseous components. Table D.1 gives the corresponding best-fit coefficients and Sect. 4.5 discuss the goodness of fit with these models. Sections 4.3 and 4.4 further discuss how the models have improved by adding gas templates other than H iand CO. We focus here first on the determination of the uncertainties associated with each parameter in the models.

The large spatial extents of the maps and the existence of tight correlations between the gas, dust, and γ-ray distributions yield small statistical errors on the best-fit coefficients. They have been inferred from the information matrix of the fit (Strong 1985) and they include the effect of correlations between parameters. The gas parameters of the local clouds have precisions of 1–5% and 1–4% in dust and γ rays, respectively. Only the small contributions of the H iiregions and of the Cetus CO clouds to the γ-ray data have respective uncertainties of 10% and 14%.

We have checked the magnitude of systematic uncertainties in our linear modelling approximations, hence of spatial changes in the model and/or in the mean level of H iand CO self-absorption. We did so by repeating the last dust and γ-ray fits a thousand times over subsets of the analysis region, namely after masking out 20% of the pixels with a sum of 2.625°-wide, randomly selected squares. In γ rays, the jackknife tests have been performed only for the total 0.4–100 GeV energy band.

Figures 5 and 6 show the distributions thus obtained for the best-fit coefficients that are significantly detected. Most of them show Gaussian-like distributions, thereby indicating that the results presented in Table D.1 are statistically stable and that the average coefficients that describe our models are not driven by subset regions in each cloud complex. Several distributions exhibit long, non-Gaussian tails when the corresponding clouds subtend small solid angles (e.g., the H iiregions or the CO clouds from the Galactic disc or Perseus). The long tails reflect the indeterminacy of the parameter when a large fraction of their maps are masked.

The standard deviations found in the jackknife distributions amount to 1–3% and 1–4% for the extended H iand DNM components in dust and γ rays, respectively. The deviations are slightly larger, respectively 1–5% and 2–10%, for the more compact CO clouds in the local ISM.

We have quadratically added the 1σ fitting error and the standard deviation of the jackknife distributions to give the statistical errors listed in Table C.1. The jackknife tests have been carried out for the total energy band. We have assumed the same relative deviations, σ_q/q, for the individual energy bands as for the total one since the jackknife samples test potential non-uniformities in the models at larger scales than the angular resolution of the data.

The results indicate that the small contribution of the faint Cetus CO clouds to the total γ-ray photon counts is only detected below 4 GeV. Similarly, the faint CO clouds from the edge of the Galactic disc are not detected against the much larger γ-ray contribution from the Galactic atomic gas (see Fig. 4). In the dust fit, all components are significantly detected.

We have checked the convergence of the iterative scheme described in Sect. 3.3 by monitoring the increase in the maximum log-likelihood value of the γ-ray fit as we progressed into the iterations until saturation. The log-likelihood ratio 2Δln(L) of 1180 between the last and first iterations indicates a large improvement in the quality of the fits because of a better separation of the different gas structures. The initial fits yielded too-large H iand CO parameters to compensate for the missing gas. As the construction of the DNM and CO_sat maps improved through the iterations, the bias on the H iand CO parameters relaxed. Between the first and last iterations, the H iand CO γ-ray emissivities dropped by 5–21% and 4–22%, respectively, depending on the cloud complex and its spatial overlap with DNM and CO_sat structures. Similarly, the dust opacities of the H iand CO clouds decreased by 10–42% and 36–87%, respectively.

4.3. Detection of DNM and CO_sat gas

Prior to building the complete models, we have checked for the presence of substantial amounts of gas that are not traced by the H iand CO line intensities or free-free emission, but that are perceptible in parallel in the independent γ-ray and dust data sets. To do so, we have performed the dust and γ-ray fits with only the H i, CO, and free-free templates for the gaseous components. We find that the data exceed the model predictions in several extended regions (see the positive residuals in the left column of Fig. 7) and that the excesses correlate in space in the independent dust and γ-ray data. Adding the DNM and CO_sat templates to the models significantly improves the quality of the fits (see the middle column of Fig. 7). The log-likelihood ratio of 6427 between the γ-ray fits with and without the DNM and CO_sat components implies that those components are detected at a formal confidence level greater than 80σ (Neyman & Pearson 1933). Hadronic interactions between cosmic rays and dust grains, or IC up-scattering of the dust thermal radiation would yield too few γ rays to explain such a correlation (Grenier et al. 2005). A coincident variation both in CR density and in dust-to-gas ratio is very improbable, so the CRs and dust jointly reveal gas in excess of the H i, CO, and ionised expectations in those directions. The spectrum of the γ-ray emission associated with both components is consistent with that found in the other gas phases over the 0.4–100 GeV energy range (see Sect. 5). This gives further support to a gaseous origin of the γ-ray emission in the DNM and CO_sat structures and to a smooth penetration of the CRs from the atomic envelopes to the dense parts of the clouds with saturated CO lines.

A log-likelihood ratio of 1404 between the models with different or with equal emissivities for the DNM and CO_sat components strongly supports their separation into distinct components. The DNM map traces column densities of dense H iand/or diffuse H₂ at the atomic-molecular transition of the clouds whereas the CO_sat map reveals dense H₂ in excess of that linearly traced by W_CO because of the saturation of the CO lines (W_CO remains constant as N_H2 keeps increasing). Without kinematic information, we cannot separate the DNM and CO_sat contributions from the different cloud complexes when they overlap in direction. To extract the column densities displayed in Fig. 8, we have assumed that the DNM and CO_sat gas is pervaded by the average CR density measured in the atomic phase of the different clouds in the region (Grenier et al. 2005; Planck Collaboration Int. XXVIII 2015). We show in Sect. 5 that the small cloud-to-cloud dispersion in γ-ray emissivity per gas nucleon justifies the use of the average. We defer the discussion of the DNM and CO_sat gas masses and their relation to the H i-bright, ¹²CO-bright, and ¹³CO-bright phases to a companion paper.

4.4. Detection of ionised gas

We have also checked the significance of the addition of the free-free template to the final γ-ray and dust models, and we have tested the gaseous origin of the corresponding γ-ray signal. We develop both points in this section.

The log-likelihood ratio of 116 between the models with and without the free-free template supports a γ-ray detection of the H iiregions at a formal confidence level greater than 10σ. The close-up view in Fig. 7 shows the clear detection of a γ-ray excess toward the bright NGC 1499 region and the more marginal gain in adding the small column densities of ionised hydrogen from the fainter G159.6-18.5 region. By replacing the free-free map by a broad Gaussian source located toward the centre of NGC 1499, we have checked that the fit significantly requires an extended excess (7σ), with a FWHM comprised between 1.4° and 3.6° at the 95% confidence level. The Gaussian source, however, yields a poorer fit than the more elongated free-free emission.

The two H iiregions, NGC 1499 and G159.6-18.5, are respectively excited by a giant and a main-sequence O star. Part of their intense stellar light is upscattered to γ rays by the local CR electrons (Orlando & Strong 2007). Given the distance, luminosity, and effective temperature of each star, we have calculated the IC γ-ray flux produced in its vicinity and we have integrated the flux over the angular sizes of the H iiregion, in the 0.4–100 GeV energy band. We find IC fluxes of ~ 2.9 × 10^-11 cm^-2 s^-1 and ~ 5.9 × 10^-12 cm^-2 s^-1, respectively, that cannot account for the γ-ray fluxes of (2.5 ± 0.4) × 10^-8 cm^-2 s^-1 and (2.3 ± 0.4) × 10^-9 cm^-2 s^-1 seen in correlation with the free-free emission. Hence, the γ-ray excesses seen toward NGC 1499 and partially toward G159.6-18.5 are more likely due to hadronic interactions of the local CRs in the ionised gas.

The weak correlation between the dust optical depth τ₃₅₃ and the free-free template highlighted by the jackknife test probably reflects the difficulty of estimating τ₃₅₃ toward regions of strong spectral gradients (see Fig. 2). In particular because the β index and the grain temperature have been derived at different angular resolutions (Planck Collaboration XI 2014).

Assuming a uniform CR flux in the H iand H iigas phases, we have used the average γ-ray emissivity of the atomic gas in the region (see Sect. 5) and the best-fit q_ff parameter obtained for the 0.4–100 GeV band to convert the free-free intensities to hydrogen column densities and compare them with the values expected from Eq. (5). The electron temperature of H iiregions is known to vary with Galactocentric radius because of the metallicity gradient (Alves et al. 2012). Adopting a mean value in the solar neighbourhood between 7000 K and 8000 K, we find an average electron density of (4.3 ± 0.6) cm^-3 in the H iiregions sampled here.

At the distance d ≈ 380 pc to ξ Persei inside NGC 1499 (van Leeuwen 2007), the 80′  radius of the ionised region is L = 9 pc. The fractional thickness of the shell, relative to L, is estimated to be l = 0.38 (Shinohara & Ishida 1976) and the nebula subtends a solid angle 4πσ with at the exciting star. The H_α flux is corrected for the interstellar absorption as F_α0 = 10^0.292A_VF_α, with A_V the visual extinction up to 500 pc derived from the 3D extinction map of Green et al. (2015). According to Shinohara & Ishida (1976) for equal electron and ion volume densities, the mean electron density of a shell-like nebula such as NGC 1499 can be expressed as with the effective volume (8)and the excitation parameter (9)The electron temperature T_e of the nebula can be estimated by comparing the free-free and H_α emissions. The H_α intensity in Rayleigh units is expressed as a function of the emission measure as: (10)Combining equations (4)and (10), we obtain (11)Given this temperature, we can estimate the mean electron density.

If we apply those calculations to NGC 1499 within a low contour at I_ff = 2 × 10⁴ Jy sr^-1, we find a mean temperature of (6100 ± 2700) K and a mean electron density cm^-3. If we adopt a more restrictive contour where the γ-ray residuals exceed 0.8σ in the best-fit model without the free-free template, we find K and a mean density cm^-3 in good agreement with the gas density estimate inferred from the correlation between the γ rays and the nebular free-free emission. This agreement gives further weight to a hadronic origin of the γ rays associated with NGC 1499.

4.5. Final model residuals

The residuals obtained between the dust and γ-ray data and the best fits that include all gaseous components are presented in the middle column of Fig. 7. They show that the linear model provides an excellent fit to the γ-ray data in the overall energy band, as well as in the separate energy bands which are not shown. The residuals are consistent with noise at all angular scales except, marginally, toward the brightest CO peaks where the model tends to over-predict the data. Significant positive residuals remain at small angular scales in the dust fit. They closely follow the spatial distribution of the DNM and they are likely due to the limitations in angular resolution and in sensitivity of the γ-ray DNM template compared to its dust homologue. Small-scale clumps in the residual structure can also reflect localised variations in dust properties per gas nucleon that are not accounted for in the linear models, in particular toward the brightest CO regions. These effects are discussed in Sect. 7.

A wide, 17°-diameter shell, centred on l = 158° and b = − 24° is apparent in the initial dust residuals, and also marginally in the γ-ray residuals (see the left column of Fig. 7). The shell is rich in DNM gas (see Fig. 8). We have found no coincident structure at other wavelengths that could explain this shell in terms of a nearby H iiregion or old supernova remnant.

5. γ-ray emissivity of the gas

The γ-ray emissivity spectra of the local gaseous components are presented in Fig. 9. We note that all the q_HI normalizations exceed one because the present fits preferred a larger H ispin temperature than the value of 140 K adopted for the derivation of q_LIS across the whole sky off the Galactic plane (Casandjian 2015). We find no significant spectral variations except in Perseus for which we observe opposite trends for q_HI decreasing and q_CO increasing with energy. This behaviour can be attributed to the cross talk between the very compact set of H iand CO clouds in Perseus as the LAT PSF degrades at low energies.

The results indicate that the CR population permeating the various phases of the different clouds has the same energy distribution as the average in the local ISM (Casandjian 2015). We find no spectral signature that would betray exclusion or concentration processes in a particular cloud or with increasing gas density, from the diffuse atomic gas in Cetus up to the densities of 10^{3 − 4} cm^-3 sampled by 2.6-mm CO line emission. Together with similar findings in the Chamaeleon clouds (Planck Collaboration Int. XXVIII 2015), these results support the theoretical predictions that CRs with energies between a few hundreds of MeV and a few tens of GeV smoothly diffuse through the H i-bright and CO-bright parts of the ISM (Skilling & Strong 1976; Cesarsky & Volk 1978; Padovani & Galli 2011), up to column densities of 10^{21 − 22} cm^-2 (Morlino & Gabici 2015). The volume densities of the gas covered by H iand CO observations are modest, but they hold the bulk of the interstellar gas mass, γ rays should reliably trace most of the interstellar gas.

We can compare the γ-ray emissivities measured in the anticentre clouds for a spin temperature of 150 K with the average value found in the local ISM for T_spin = 140 K (Casandjian 2015), as shown in Fig. 10. The 35 ± 1% lower emissivity that we find for the atomic gas in the Galactic-disc background is not useful for comparison as it integrates all distances beyond the local ISM, through regions of decreasing CR density along lines of sight outward through the Galaxy (Abdo et al. 2010; Ackermann et al. 2011) and to high altitudes above the Galactic disc (Tibaldo et al. 2015). Within a few hundred parsecs, we find small cloud-to-cloud differences in γ-ray emissivity. The Cetus, South Taurus, and North Taurus clouds compare well with the average in the local ISM whereas the California, Main Taurus and Perseus clouds appear to be 15 ± 3%, 1 ± 2%, and 36 ± 6% more emissive than the average, respectively. Except for Perseus, these variations are commensurate with previous measurements. Figure 10 shows a cloud-to-cloud dispersion of ± 9% that does not relate to the cloud altitude above the Galactic disc, nor to the cloud location with respect to the local spiral arm (Galactocentric distance). This dispersion largely stems from uncertainties in the derivation of the N_HI column densities. We note indeed that the q_HI emissivity is 10% to 20% larger in the Cepheus-Polaris, Main Taurus, California, RCrA and Chamaeleon clouds where N_HI column densities in the 80th percentile exceed 10²¹ cm^-2, a level at which H ioptical depth corrections become significant.

The emissivity in Perseus is at variance with the other nearby clouds, but less so at high energy where the LAT PSF allows a firmer separation between the compact atomic and molecular phases of the cloud. The origin of the high emissivity will be investigated using photon selections with improved angular resolution as soon as the photon statistics permit the analysis. For the anticentre region under study, we derive very consistent averages regardless of whether we include the Perseus data to the other five measurements or not. We have used this average emission rate in the 0.4–100 GeV band, s^-1 sr^-1 H^-1, in the atomic gas to infer column densities in the other gas phases.

6. X_CO factors

The X_CO factor is often applied to convert W_CO intensities to hydrogen-equivalent N_H2 column densities in the molecular phase. In γ rays, under the assumption of a uniform CR flux in the H i-bright and CO-bright phases, it is derived as X_COγ = q_CO/ (2q_HI), independently in each energy band. This ratio also assumes that helium is uniformly mixed by number with hydrogen in the ISM. Likewise, assuming a uniform dust-to-gas mass ratio and a uniform emission coefficient κ₃₅₃ of the grains at 353 GHz, the X_CO factor is derived in the dust fit as X_COτ = y_CO/ (2y_HI). The results of the two methods for the six CO clouds selected in our analysis are given in Table 1. The γ-ray results in the individual energy bands are displayed in Fig. 11.

6.1. X_CO measurements in γ rays

We find X_CO values for the present clouds that compare well with previous γ-ray measurements in other nearby clouds. For X_CO factors obtained with Fermi-LAT at comparable angular resolution, we can cite the values of in the Cepheus-Polaris complex (Abdo et al. 2010), of in the RCrA cloud (Ackermann et al. 2012c), of 1.21 ± 0.02 in the Orion clouds (Ackermann et al. 2012b), and of in the Chameleon clouds (Planck Collaboration Int. XXVIII 2015), in units of 10²⁰ cm^-2 K^-1 km^-1 s. We note that the X_CO factors in those well-resolved clouds are all close to or lower than 10²⁰ cm^-2 K^-1 km^-1 s. Table E.1 lists historical γ-ray estimates of X_CO. It shows that the estimates in nearby clouds have remained stable or only slightly decreased as the angular resolutions of the H i, CO, and γ-ray surveys have improved. The lower level of cross-correlation between the H iand CO structures on the one hand, a better separation of the γ radiation produced in the atomic and molecular phases on the other hand, and a lower contamination from unresolved γ-ray point sources has primarily improved the precision of the measurements. The latest analyses with <10.8′  FWHM resolution in H i, 8.5′  resolution in CO, and Fermi-LAT data consistently yield low X_CO values, ranging between 0.5 ×10²⁰ cm^-2 K^-1 km^-1 s and 1.2 ×10²⁰ cm^-2 K^-1 km^-1 s. Systematic uncertainties on the N_HI column densities due to H ioptical depth corrections impact the estimation of the γ-ray emissivities per gas nucleon. Analyses for different H ispin temperatures show that the resulting systematic uncertainty on the X_CO factors amount to (0.1–0.2) ×10²⁰ cm^-2 K^-1 km^-1 s.

The values listed in Table E.1 highlight a recurrent discrepancy between the X_CO measurements at small scale in nearby clouds and the two-to-three times larger values found at kiloparsec scales when averaging over a population of molecular complexes in spiral arms or in Galactocentric rings. We return to this point at the end of this section.

In Fig. 11 we take advantage of the energy-dependent PSF of the LAT to probe various linear scales within the present nearby clouds. The scale is derived using the half-width at the half-maximum of the PSF integrated over the energy band for the q_LIS spectral shape. We find no evidence of X_CO changes at parsec scales except in Perseus.

Because the q_HI and q_CO variations in Perseus are inversely coupled in energy, the resulting change in X_CO is likely due to an increased level of cross-correlation between the compact H iand CO phases as the LAT PSF degrades (Figs. 4 and 9). If we force the γ-ray emissivity of the H igas to be the same as the average found among the other clouds, we obtain a larger X_CO value of (0.68 ± 0.04)×10²⁰ cm^-2 K^-1 km^-1 s, close to that measured above 4 GeV with the best LAT angular resolution. However, the use of this value implies a significant (4.6σ) degradation of the fit quality at lower energies. The likelihood analysis significantly supports a lower X_CO factor in Perseus. A low value is also indicated by the average of 0.3 ×10²⁰ cm^-2 K^-1 km^-1 s found at 0.4-pc resolution in the dust, H i, and CO study of Lee et al. (2014).

Several observations and numerical models have suggested that the X_CO factor changes with the density, turbulent velocity, and metallicity of the gas, and with the penetration of the CO-dissociating radiation inside a cloud (Polk et al. 1988; Bell et al. 2006; Sheffer et al. 2008; Liszt et al. 2010; Glover & Mac Low 2011; Shetty et al. 2011a; Liszt & Pety 2012; Lee et al. 2014; Bertram et al. 2016). One therefore expects intrinsic X_CO gradients in clouds. Calculations by Bell et al. (2006) show that X_CO should dramatically drop inward from the outskirts of a molecular cloud as the photo-dissociating radiation gets absorbed. The X_CO factor drops down to a minimum value at A_V around 2 mag. The X_CO ratio then rises again in the denser cores because of the saturation of the optically thick CO lines. One goal of our analysis is to explore whether these intrinsic variations can change the mean X_CO value averaged over a cloud, depending on how diffuse it is. We cannot yet measure the steepest X_CO gradients at the periphery of clouds (A_V≪ 1 mag) near the sensitivity threshold of CO surveys because of the complex 3D interface with the abundant atomic phase. We cannot map X_CO gradients pixel by pixel in γ rays because of the photon statistics, but our sample includes more or less translucent clouds of the local ISM and we can explore changes in the mean value of X_CO for clouds in different states.

Figure 12 compares γ-ray measurements of X_CO in nearby clouds from this work and Planck Collaboration Int. XXVIII (2015). We restrict the sample to estimates obtained in the Chamaeleon and anti-centre clouds because they are based on the same analysis method, the same γ-ray energy bands to ensure the same angular resolution of the LAT, and the same sampling resolution in the H iand CO data. They both include the DNM gas in the models. These clouds subtend large solid angles in the sky and their column density structures are well resolved.

The distributions in Fig. 12 show that the average X_CO factor per cloud does not depend on the H₂ mass mapped in CO, nor on the overall dynamics of the cloud characterised by the velocity dispersion of the CO lines. But X_CO appears to depend on the cloud structure and in particular on its diffuseness. To reflect changes in the latter, we have explored several characteristics:

• the mean W_CO intensity, , in a cloud, taken above1 K km s^-1 to avoid noise fluctuations;

• the surface fraction of dense regions with large W_CO intensity within a cloud, SF_dense = S_WCO> 7 K km s^-1/S_WCO> 1 K km s^-1 (with S a solid angle), which gauges the relative weight of diffuse and dense molecular regions in the determination of X_CO;

• the mean visual extinction toward a cloud, , taken from the NICER M2a 12′-resolution A_J map (Juvela & Montillaud 2016), translated into A_V using a colour ratio of 3.55 according to the extinction law of Cardelli et al. (1989). Only values in directions with W_CO> 1 K km s^-1 and A_V> 0.8 mag. have been retained in the average to avoid noise contamination;

• the average CO line width, , in a cloud since more CO line photons can escape the cloud and ultimately be detected as the line width grows. The resulting increase in W_CO intensity implies a decrease in X_CO with (Shetty et al. 2011b).

Figure 12 shows that, for all these diagnostics, X_CO tends to decrease from diffuse to more compact clouds. The change in X_CO is approximately a factor of two. The decrease remains significant if we use the Perseus X_CO value found when forcing the H iγ-ray emissivity to the local average, or if we exclude the Perseus data point. A uniform X_CO factor is strongly rejected (> 30σ) for all the diffuse diagnostics, whether we modify or exclude the Perseus data point or not.

The trend is consistent with the expectation that X_CO decreases from the diffuse envelopes, exposed to the photodissociating radiation field, to the shielded cores of the CO clouds (Bell et al. 2006; Glover & Mac Low 2011; Shetty et al. 2011a; Bertram et al. 2016). Hence, the average X_CO factor per cloud in our sample is seen to vary with the relative solid angle subtended by the diffuse and dense parts of the clouds, reflected in SF_dense and . The turbulence level and the degree of virialization of the cloud (e.g., the ratio of kinetic to gravitational energy) also influences the average X_CO factor by allowing the ISRF to penetrate more or less deeply into the turbulent layers of the cloud (Bell et al. 2006; Shetty et al. 2011b; Bertram et al. 2016). The mean X_CO factor per cloud in Fig. 12 does decrease with the average line width of the CO lines. The present measurements also give weight to the idea that the mass or size of a cloud does not play an important role in determining the X_CO factor (Shetty et al. 2011b).

According to X_CO calculations as a function of visual extinction A_V (Bell et al. 2006), X_CO should drop by more than two orders of magnitude with increasing A_V, down to a minimum at 1 ≲ A_V ≲ 3 mag. For long-lived clouds, exposed to the local cosmic-ray ionization rate of 1.4 × 10^-17 s^-1 (Grenier et al. 2015) and to the local ISRF, the minimum X_CO value reaches below 10²⁰ cm^-2 K^-1 km^-1 s for hydrogen densities larger than 10³ cm^-3 and/or for turbulent velocities larger than 1 km s^-1. Yet, the measured CO lines have widths below 0.7 km s^-1 and the clouds in our sample are not massive enough to sustain such large average densities through most of their volume. Moreover, models of dense clouds show a steady increase in X_CO at all N_H2 column densities (Shetty et al. 2011a), at variance with the declining trends found in Fig. 12.

We have also compared the mean visual extinction, , mean CO intensity, , and X_CO factor in each cloud with the simulation results of Glover & Mac Low (2011). The trend seen in Fig. 12 roughly compares in slope with the simulation expectations, but not in absolute values, even though the efficiency of the radiation screening compares well in the observations and simulations. We show in Sect. 7.1 that the A_V/N_H ratio measured in the clouds compares well with the Bohlin et al. (1978) value adopted in the simulations. For a given in a cloud, the simulations under-predict the mean intensity by a factor 20, and overpredict X_CO by a factor 40. The comparison should be taken with care because both the observation and simulation samples are sparse, because the observations span a smaller dynamical range in than the simulations, and because the and averages depend on the threshold applied to define CO cloud edges, which differs in both datasets. Nevertheless, the comparison highlights that, for initial volume densities ≤ 300 cm^-3 applicable to the set of observed clouds, the simulations produce underluminous CO clouds at low densities compared to the observations. This was noted by Levrier et al. (2012). The CO deficit in their simulation reached a factor of ten in column density in the peripheral regions exposed to the UV radiation, where N_H2≲ 2 × 10²⁰ cm^-2. At larger column densities, the simulations and data compare much better. So, there seems to be a general problem with our understanding of chemistry at low gas density in the photo-dominated regions (PDR), either in the treatment of the UV attenuation or in the absence of warm chemistry driven by intermittent energy deposition in regions where the interstellar turbulence is dissipated (Godard et al. 2014).

The list of historical X_CO measurements given in Table E.1 and the values we are adding with the present sample of anti-centre clouds confirm a recurrent discrepancy that exists between the mean X_CO factors measured with parsec resolution in nearby clouds and the X_CO averages obtained at the scale of spiral arms (Local and Perseus arms) or of the Galactic disk. The former are two to three times lower than the large-scale values, which are close to 2 ×10²⁰ cm^-2 K^-1 km^-1 s (see Table E.1). Simulations of molecular clouds have warned that X_CO factors can vary by up to a factor of five in the Milky Way because of the complexity of the overlap and blending in position and velocity of the CO clouds emitting along the lines of sight (Shetty et al. 2011b). Smith et al. (2014) have estimated an average value of 2.2 ×10²⁰ cm^-2 K^-1 km^-1 s, close to the large-scale γ-ray measurements, when including all regions where CO intensities are potentially observable (i.e. wherever >0.1 K km s^-1). But the spatial distribution of H₂ column densities and CO intensities in their Galactic simulation implies significant deviations from this average, depending for instance on the location of cloud complexes with respect to the spiral arms.

Because of intrinsic X_CO gradients across clouds, X_CO averages in large-scale observations respond to the mix of clouds in different states present in the region under study. At least three potential reasons can explain why the large-scale X_CO factor can exceed the local values. A first possibility is that a large fraction of the molecular mass lies in diffuse, gravitationally unbound clouds where X_CO factors can reach very large values of 10²¹⁻²² cm^-2 K^-1 km^-1 s (Bell et al. 2006; Glover & Mac Low 2011). A second possibility is that the average is driven by a large number of giant molecular clouds which are expected to exhibit X_CO values of a few 10²⁰ cm^-2 K^-1 km^-1 s because the gas is predominantly in the optically-thick CO regime. A third influence is the level of shearing in the clouds. For instance, clouds sheared out by Galactic differential rotation are more exposed to photo-dissociation and should exhibit larger X_CO values (Smith et al. 2014).

In parallel, the change in X_CO estimates with scale may be due to determination biases induced by the sampling resolution in the gas maps, such as the increased difficulty at large distance to separate the clumpy CNM clouds and the DNM envelopes from their CO-bright phase. Misattributing a small part of the dense atomic gas or diffuse H₂ would significantly bias the X_CO factor upward. For instance, confusing 10–20% of the CNM with the CO-related signal could explain the observed change in X_CO by a factor of two between the local ISM and spiral arms (Planck Collaboration Int. XXVIII 2015). In nearby galaxies, Sandstrom et al. (2013) have indeed measured a systematically larger X_CO factor in highly inclined galaxies than in face-on ones where the pile-up along sight lines is reduced. More tests are needed to disentangle the origin of the X_CO changes with scale. With larger photon statistics at high γ-ray energies, we should be able to investigate how the angular resolution and confusion between gas phases affect the calibration of X_CO beyond the solar neighbourhood.

6.2. X_CO measurements with dust

The results of the dust model provide an independent measure of X_CO. Figure 13 shows that the dust-derived X_CO factors are systematically 30% to 130% larger than the γ-ray values. Systematic discrepancies had already been noted in previous studies (Planck Collaboration XIX 2011; Bolatto et al. 2013), but on the basis of different H iand CO surveys, different correlation methods, and different angular resolutions. In the present study, as in Planck Collaboration Int. XXVIII (2015), the dust and γ-ray models are fitted to the same N_HI and W_CO maps, at the same resolution. Bolatto et al. (2013) proposed that the inclusion of the DNM gas in γ-ray analyses and its absence in dust analyses could explain the lower X_CO values obtained in γ rays, but, by construction, the inclusion of the DNM map cannot bias the X_CO estimate (Grenier et al. 2005; Planck Collaboration Int. XXVIII 2015). Furthermore, the recent analyses include the DNM in both the dust and γ-ray fits and we have verified in this work that the X_CO values change by less than 20% when the DNM phase is added or not in the fits.

Alternatively, the standard derivation of the X_CO factor as X_COτ = y_CO/ (2y_HI) is based on a uniform dust opacity, τ₃₅₃/N_H, across the gas phases. A change in dust emission properties, in particular an increase in τ₃₅₃/N_H at large gas column densities, can significantly bias the X_CO factor upward (Planck Collaboration Int. XXVIII 2015). We show in Sect. 7.2 evidence for such a departure from linearity in dense CO regions where τ₃₅₃/N_H can increase by a factor of up to six due to grain evolution compared to their properties in the atomic gas. Outside the densest molecular parts (brightest in CO), the molecular hydrogen is more diffuse and the increase in τ₃₅₃/N_H is more modest, by a factor ranging between 1.2 and about 2, so the resulting bias on the derivation of the average value of X_CO across a cloud is less than a factor of two.

Dust reddening and extinction maps have also been used to evaluate X_CO in nearby clouds. We find values in close agreement with the extinction study of Chen et al. (2015) in the main Taurus and Perseus clouds (their Taurus and Perseus E1 values). The present Taurus South complex overlaps their Taurus E1, E2, and E3 structures with respective X_CO factors of (0.84 ± 0.01), (1.41 ± 0.02), and (1.69 ± 0.04)×10²⁰ cm^-2 K^-1 km^-1 s bracketing our result in Taurus South.

Using the 2MASS extinction map of Dobashi et al. (2013), Paradis et al. (2012) found X_CO= (2.27 ± 0.9)×10²⁰ cm^-2 K^-1 km^-1 s in the Taurus-Perseus region, without decomposing the different clouds. Their A_V/N_H estimate of 5.33 ×10^-22 mag cm² in this region is 20% lower than the all-sky | b | > 10° average, but this difference is not sufficient to reconcile their measurement based on dust reddening with our values based on dust emission. The determination of A_V/N_H in the H iphase directly impacts the derivation of X_CO in the embedded molecular phase. The different methods used to infer A_V from stellar data yield systematic differences as large as 1 mag at low A_V, toward the atomic and diffuse ISM (Dobashi et al. 2013; Green et al. 2015; Juvela & Montillaud 2016), so further analyses are necessary to evaluate their impact on X_CO measurements.

Pineda et al. (2010) have used 2MASS data and the NICER method to map the dust reddening and to derive X_CO near 2.1 ×10²⁰ cm^-2 K^-1 km^-1 s in the Taurus Main complex. In this case, the difference with our result can be explained because they did not separate the extinction related to H ialong the CO sight lines, but they attributed all the dust to the molecular phase. In our analyses, we separate the dust or γ rays associated with each phase along all sight lines.

7. Dust evolution

7.1. Average dust opacities per gas phase

Our analyses yield average dust properties per H atom in each gas phase. The best-fit y_HI coefficients directly give the average dust opacities, , in the H istructures. We have used the γ-ray estimates of the X_CO factors to derive the average dust opacities in the CO-bright phase of each complex, = y_CO/ (2X_COγ). The results are listed in Table 2 together with similar measurements in the Chamaeleon clouds (Planck Collaboration Int. XXVIII 2015). To ease the comparison with other dust reddening measurements, we have converted the dust opacities to N_H/E(B − V) ratios using the E(B − V) /τ₃₅₃ slope of (1.49 ± 0.03) × 10⁴ found in a large correlation study between τ₃₅₃ and reddening measurements toward 53 399 quasars (Planck Collaboration XI 2014). We did not attempt to estimate dust opacities in the ionised gas because of the large temperature gradients which impact the τ₃₅₃ derivation toward the H iiregions, with uncertainties in τ₃₅₃ doubling near the hotspots.

The correlation existing between the γ-ray intensity and dust optical depth in the DNM yields an important measure of the average dust opacity in this phase. The value = q_HI/q_DNM assumes a uniform CR flux across the H iand DNM phases, an assumption that is corroborated by the fact that both phases exhibit comparable γ-ray emissivity spectra and moderate gas volume densities. According to Fig. 14, we find no opacity changes at parsec scales in the extended DNM structures. The average, = (16.5 ± 0.3) × 10^-27 cm², compares very well with the value of (17.2 ± 0.5) × 10^-27 cm² obtained in the DNM surrounding the Chamaeleon clouds (Planck Collaboration Int. XXVIII 2015).

Using optically thin H i, Planck Collaboration XI (2014) found marked gradients in dust opacity across the sky, with values ranging from 6.6 to about 11 ×10^-27 cm²in the atomic gas, and a high-latitude value of (7.0 ± 2.0)×10^-27 cm²in good agreement with the opacity of (7.1 ± 0.6)×10^-27 cm²derived in high-latitude cirruses with a different analysis method (Planck Collaboration Int. XVII 2014). The opacities continued to increase in the CO-bright phase, up to 18 ×10^-27 cm²for an X_CO factor of 10²⁰ cm^-2 K^-1 km^-1 s. Our results confirm that dust opacities vary from cloud to cloud and they can reach values well in excess of the high-latitude value, both in the atomic and molecular phase.

We find 50% to 135% larger opacities in the anticentre H iclouds than in the more diffuse H icirruses, even though we have corrected the N_HI column densities for a spin temperature of 400 K. The systematic uncertainties on τ₃₅₃/N_H due to the H ioptical depth correction are listed in Table 2 as the second set of errors. Given the likelihood variations shown in Fig. C.1, we have taken the standard deviation in the set of measurements as a function of spin temperature for each parameter to represent such an uncertainty. Considerations of those uncertainties cannot reconcile the individual cloud measurements and the cirrus value, even for the rather tenuous H iclouds in Cetus and Taurus North. The cloud-to-cloud dispersion indicates environmental effects beyond the latitude dependence noted by Liszt (2014) in the atomic gas. The opacities listed in Table 2 do not relate to the H imass of the corresponding clouds. The 15% dispersion in this sample rather relates to the mean visual extinction of the H icloud and to the surface fraction of large and small N_HI column densities across the cloud. We present in the next section evidence for gradual opacity changes as the gas density increases that explain why the average opacity in an H icloud depends on its structure. The observation of opacity variations within the atomic phase at column densities below 10²¹ cm^-2 challenges the current dust evolution models discussed in the next section.

We find systematically larger mean opacities in the CO-bright phase than in the H iphase, thereby confirming the continued increase in τ₃₅₃/N_H at large H₂ column densities. The latter has been evaluated independently in each cloud (independent X_CO factors), so the cloud-to-cloud changes in that Planck Collaboration XI (2014) have noted across the sky are not due to their use of a uniform X_CO factor. In Table 2, the mean opacity increases by 30% to 100% between the H iand CO-bright phases, independent of the H₂ mass locked in the cloud. The dispersion in rather relates to the surface fraction of dense regions with large W_CO intensities, SF_dense, in the averaging. To study environmental changes in τ₃₅₃/N_H within molecular clouds, we must take into account the added complication of spatial gradients in X_CO.

7.2. Environmental changes in dust opacity

In order to follow changes in dust opacity across the clouds, we have built two maps of the total N_H column density in the region. The first one, N_{H γ}, takes advantage of the CR interactions with gas in all chemical forms and thermodynamical states. The γ-ray intensity from the gas has been obtained from the LAT data in the overall energy band after subtraction of the γ-ray counts unrelated to gas in the best-fit model. In order to reduce the Poisson noise and get more robust photon statistics, we have degraded the resolution to 0.̊375. We have converted this γ-ray intensity into N_H using the average emission rate, , obtained in the atomic phase in the region. The second map, N_{H mλ}, uses the higher-resolution information from the H i, CO, and dust data, and the mass scaling provided by the best-fit γ-ray model (uniform H ispin temperature of 400 K, X_COγ factors in the clouds, γ-ray emissivity in the DNM, ionised, and CO_sat gas as in the H i). The opacity distribution at the smallest angular scales should be considered with care because the scaling factors used to construct N_{H mλ} have been derived at the scale of a whole cloud.

Figure 15 shows the 2D histogram of the correlation between the sightline integrals in dust optical depth and in N_{H γ} column density. The distribution deviates from linearity above about 4 × 10²¹ cm^-2. The opacity maps, shown in Fig. 16, indicate that τ₃₅₃/N_H evolves progressively within the clouds, especially when the medium becomes molecular, in the diffuse DNM and even more so in the CO-bright parts. The increase in opacity does not relate to an increase in the specific power, 4πR/N_H, radiated by the grains, so the rise in opacity cannot be attributed to a larger heating rate (if we consider that the large grains are in thermal equilibrium to equate the radiated and absorbed powers). The largest opacities rather correspond to a ~ 30% lower heating rate in the well-shielded CO regions.

Thanks to the total gas tracing capability of the CRs and to the complementarity of the N_{H γ} and N_{H mλ} maps, it is very unlikely that opacity variations as large as those seen in Fig. 16 are caused by large deficits in the gas column densities. Our results give strong support to earlier indications that the big grains in the CO filaments of the Taurus cloud have larger sub-mm emissivities than in the more diffuse media. The previous studies noted opacity changes by a factor ranging from ~ 2 (Flagey et al. 2009; Planck Collaboration XXV 2011; Ysard et al. 2013) to (Stepnik et al. 2003) above the diffuse-ISM value. In the present analysis, we find opacity enhancements exceeding a factor of three and reaching a factor of six toward the CO clouds, even though we provide a measure of the additional H₂ gas that is not linearly traced by W_CO in the CO_sat directions. The largest opacities do not rise as high in the Taurus clouds as they do in the Chamaeleon ones (Planck Collaboration Int. XXVIII 2015), but a quantitative comparison would require the separation of the DNM and CO_sat components in a new analysis of the Chamaeleon region. We also find that the opacity starts to increase in the DNM, therefore at lower gas densities than the few thousand per cm³ sampled in CO. The grain opacity at the H i–H₂ transition compares very well with that in the DNM surrounding the Chamaeleon clouds (Planck Collaboration Int. XXVIII 2015).

In order to determine how the opacity gradient depends steeply on the ambient ISM, Fig. 17 shows how the opacity varies with the molecular fraction in the gas column, f_H2, with the total N_{H mλ} column density as traced by the multi-wavelength data, and with the dust colour temperature. We observe a gradual increase as the gas becomes molecular and a steeper change at N_H above 5 × 10²¹ cm^-2 due to two factors: change in molecular fraction and in grain properties in the molecular gas. The amplitude of the gradient is best traced by the drop in grain temperature. Figure 16 further shows that the ambient molecular gas density is not the only factor of evolution because we detect notable cloud-to-cloud differences in the opacity gradients. Compare for instance the large opacities that often exceed 25 ×10^-27 cm² along the Main Taurus and Perseus CO clouds with the 50% lower values obtained toward comparable W_CO intensities near both ends (in longitude) of the California cloud.

Several studies have noted that the opacity enhancement occurred together with a reduction of the IR emission from polycyclic aromatic hydrocarbons (PAHs) and from very small grains stochastically heated by the ambient ISRF (e.g., Stepnik et al. 2003; Flagey et al. 2009; Ysard et al. 2013). Grain evolution in dense media was invoked to explain the radiation changes (see the previous references and Martin et al. 2012; Köhler et al. 2015). It involves accretion of aliphatic-rich carbonaceous mantles from the gas phase, grain coagulation into fluffy aggregates, and the formation of an ice mantle on the surface of the aggregates in the coldest regions. Recent theoretical studies such as Köhler et al. (2015) suggest that the spectral changes in β relate to gas density-dependent processes and that the variations in opacity, spectral index β, and colour temperature are closely connected. Figure 18 shows the quantitative relation between those three parameters in the observations. The data corroborate the prediction that the opacity should increase for β above 1.5 and temperature below 17 K. According to this model, the directions with temperatures near 16–18 K, β ≥ 1.5, and τ₃₅₃/N_Hmλ values about twice as large as the diffuse-ISM value of ~ 7×10^-27 cm²would be rich in grains with carbonaceous mantles. The mantle accretion would therefore occur primarily in the DNM. Aggregates would dominate the directions below 16 K, at 1.6 ≲ β ≲ 1.8, and where τ₃₅₃/N_H is at least four times above the diffuse-ISM value. Such conditions are found inside the sampled CO clouds. Ice coating of the aggregates would explain the largest opacity values in Fig. 18, at the lowest temperatures and largest spectral indices, β.

8. Conclusions

We have analysed the gas, dust and CR content of several nearby anti-centre clouds including the Cetus, Taurus, Auriga, Perseus, and California clouds. We have performed an iterative fit of the total gas column density as traced by γ rays and the dust optical depth at 353 GHz. We have modelled both tracers as linear combinations of the gas column densities traced by H i, ¹²CO, and free-free emissions, and we have extracted the morphology and column density of the DNM from the joint γ-ray and dust data. We have verified the robustness of our set of parameters through the use of jackknife tests; these parameters have enabled us to constrain the properties of the ISM in these clouds.

The main results are summarised as follows.

• On the cosmic rays. The γ-ray emissivity of the gas in the analysedclouds has the same energy spectrum as in other clouds of the localISM and we find no dependence of their average γ-ray emissivitywith Galactocentric radius (i.e., distance from the local spiralarm), nor with height above the Galactic plane. In the0.4–100 GeV energy band and at the precisionlevel of the current LAT observations, we find no evidence of CRexclusion or CR concentration in the clouds, up to the ¹²CO-bright molecular regions.

• On the dark gas. The γ-ray and dust data jointly reveal significant amounts of gas in addition to that seen in H i, free-free, and ¹²CO emissions. The diffuse large-scale structures are associated with the DNM at the transition between the atomic and molecular phases. They gather dense atomic hydrogen and diffuse H₂ in unknown proportions, but with column densities equivalent to those found in the H iand CO emitting parts. In the molecular phase, the γ rays and dust reveal filaments of dense gas in addition to that proportionality traced by W_CO where the ¹²CO line intensities saturate.

• On the H iiregions. The H iiregions, NGC 1499 and G156.6-18.5, are jointly detected in dust and γ-ray emission. The corresponding γ-ray flux cannot be attributed to the IC up-scattering of the bright stellar radiation by the local CR electrons. It is more likely caused by the hadronic interactions of the local CR nuclei in the ionised gas. Hence for NGC 1499, we have used the average γ-ray emissivity of the atomic gas in the region to infer a mean electron density of 4.3 ± 0.6 cm^-3 for an electron temperature of 8000 K, in good agreement with the modelling of the free-free and Hα emissions of the H iiregion.

• On the X_CO factors. We provide independent measurements of the X_CO factor from the dust and γ-ray analyses in six different clouds, against the same H iand CO data. As in the Chamaeleon complex, we find that the dust-derived values are systematically larger, by 30% to 130%, than the γ-ray estimates. The difference is likely due to chemical and structural evolution of the dust grains with increasing N_H. The X_CO factors measured in γ rays range from (1.04 ± 0.05)×10²⁰ cm^-2 K^-1 km^-1 s in Taurus South to (0.44 ± 0.04)×10²⁰ cm^-2 K^-1 km^-1 s in Perseus. Together with the estimate found in the Chamaeleon clouds (Planck Collaboration Int. XXVIII 2015), these measurements indicate that X_CO values tend to decrease with the average W_CO intensity of a cloud, with the average CO line width, and with the surface fraction subtended by the brightest CO clumps. The more diffuse CO clouds therefore tend to have larger average X_CO factors. Models of the formation and photodissociation of H₂ and CO molecules predict a marked decline in X_CO from the diffuse envelopes of molecular clouds to their dense cores. Hence the X_CO average should qualitatively vary with the surface fraction of dense regions in the cloud structure and with the CO line width, as in our data. We find a modest change by typically a factor of two in X_CO with cloud state, but our sample lacks a giant molecular cloud to extend the state range. The amplitude of these variations already limit the precision of molecular mass estimates based on the use of CO intensities and of a mean X_CO conversion factor. We further note that the low X_CO values measured here are at variance with the theoretical predictions from PDR chemistry for the moderate gas densities filling most of the volume of the observed clouds. The present measurements confirm a recurrent discrepancy by a factor of two to three between the X_CO averages obtained in nearby clouds and at large scale in the Galaxy. Simulations predict X_CO variations of this amplitude in the Galactic disc owing to the variety of CO-line blending situations occurring along the lines of sight (Shetty et al. 2011b; Smith et al. 2014). A dominance of diffuse, unbound or sheared clouds, or of giant molecular clouds, can drive the X_CO factor of the ensemble to large values. Additionally, the increased level of cross-correlation at large distance between the CNM, DNM, and CO-bright phases of the ISM could bias the large-scale X_CO determination upward.

• On the dust opacities. The dust opacity at 353 GHz, τ₃₅₃/N_H, appears to rise by a factor of three from low column densities in the atomic gas (about 10^-26 cm²) to cold grains in molecular gas at N_H≳5 × 10²¹ cm^-2. The rise can reach a factor of six at very low dust temperatures below 14 K. The amplitude of the rise is comparable to the variations observed in the Chamaeleon clouds (Planck Collaboration Int. XXVIII 2015). As the observed specific power radiated by the grains decreases in the cold molecular regions the changes cannot be attributed to a larger heating rate, but they are more likely caused by a chemical or structural change in the grains. The magnitude of the rise severely limits the use of the thermal emission of the large dust grains to trace the total gas. The linear regime is limited to N_H<3 × 10²¹ cm^-2 in these clouds. The confirmation of large opacity variations across clouds directly impacts the gas mass estimates inferred from dust emission at sub-mm and mm wavelengths to derive star-forming efficiencies in the Galaxy and in external galaxies. We quantify the coupled changes in grain temperature, opacity, and β index of the thermal radiation to help model the grain evolution.

In a second paper we will study how the DNM and CO_sat gas relate to the H i-bright, ¹²CO-bright, and ¹³CO-bright phases. We will derive the masses of the DNM and of the CO_sat phase for a sample of substructures in the Cetus, Taurus, California and Chameleon clouds, while taking care to avoid regions of cloud confusion along sightlines. We will follow the evolution of the CO-dark H₂ fraction both within a cloud and from cloud to cloud and we will place our results in the context of other observational and theoretical studies.

In a subsequent study, we will confront the gas, γ-ray, and dust distributions in the nearby clouds using dust column densities inferred from stellar reddening. The joint comparison of dust emission and extinction data to other gas tracers will provide additional information on the relation between the total gas and dust, including the abundant dark neutral medium, and will help to better constrain the evolution of dust properties, in particular with regard to potential variations of the extinction curve (Schlafly et al. 2016).

Adding more data points to the trends shown in Fig. 12 is essential to confirm that the average X_CO factor in a cloud tends to decrease for increasingly dense CO molecular clouds, in response to different X_CO gradients spanning each cloud (Bell et al. 2006; Glover & Mac Low 2011). Particular attention should be paid to the use of the same method to derive a consistent set of X_CO values, and to the addition of clouds with little confusion between the different gas phases to ensure a precise separation of the total gas tracers in each phase. The current γ-ray analyses cannot yield X_CO maps across individual clouds. The variations in dust opacity with increasing gas column density prevent the use of dust emission maps to probe X_CO gradients in nearby clouds. In this context, it is essential to increase the number of clouds with measured average X_CO values and to confront the observed X_CO trends with PDR models in order to interpret the slope and dispersion of the trends presented here. Confirming that X_CO factors vary from cloud to cloud within the local ISM because of their different dynamical and chemical structures stresses that large-scale X_CO averages can vary because of a different mix of cloud states in different Galactic regions, in addition to the Galaxy-wide X_CO gradients expected from metallicity and ISRF gradients. For instance, it opens the way for observable differences in X_CO between ensembles of clouds inside and outside of the spiral arms since the inter-arm clouds that have been sheared after their passage through an arm are more susceptible to photo-dissociation, thus more prone to having large X_CO averages (Smith et al. 2014).

---

Acknowledgments

We thank the referee for an interesting discussion. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Énergie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. The authors acknowledge the support of the Agence Nationale de la Recherche (ANR) under award number STILISM ANR-12-BS05-0016.

References

Appendix A: Component separation

We have decomposed the H iand CO spectra into sets of individual lines as in Planck Collaboration Int. XXVIII (2015). The method starts with the detection and localization in velocity of isolated lines or significant peaks in the spectra. Each HI and CO spectrum is then fitted by a sum of pseudo-Voigt line profiles centred on the centroids of the detected lines or peaks, within ± 2.7 to 3 km s^-1. The fits are generally good, with less than 15% difference between the observed and fitted integrals over all velocities. In order to preserve the observed photometry, we calculate the residuals between the observed and fitted spectra in each channel and we redistribute them among the fitted lines, proportionally to their intensity in this channel. We finally construct N_HI column-density and W_CO intensity maps of a specific cloud by integrating the fitted lines that have their velocity centroid within a chosen velocity interval for each (l, b) direction. This method corrects for potential line spill-over from one velocity interval (i.e., cloud) to the next. To determine boundaries between clouds in the position-velocity space (l, b, v), we have used the velocity centroids of all the fitted lines and their density in the (l, b, v) space to search for line clusters (i.e., clouds). To do so, we have used the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm (Ester et al. 1996) to isolate arbitrarily-shaped clusters and to draw contours of minimum density around them. The six H iand CO regions used in this analysis have been delimited by following those curves of minimal number density in (l,b,v) space. The cuts follow piecewise continuous lines in (b,v) over specific longitude intervals. The latter and the points used to define the broken lines are given in Table A.1.

Appendix B: Free-free emission in the 70 GHz data

In order to separate the free-free intensity in the Planck LFI map at 70 GHz from non-gaseous emission components, we have proceeded as follows. We have first subtracted the CMB intensity as modelled at 70 GHz by Bobin et al. (2016) from the joint analysis of WMAP and Planck observations. We have then taken advantage of the sparse spatial correlation between the thermal dust intensity I₇₀ at 70 GHz and the dust optical depth τ₃₅₃ at 353 GHz in the wavelet domain to filter out the dust component. To do so, we have transformed both maps with eight scales of the B-spline wavelet with the mr-transform tool of Starck & Pierre (1998). For each scale, we have obtained the best (χ²) linear regression between the wavelet distributions of I₇₀ and τ₃₅₃ in the regions of both significant dust optical depth (τ₃₅₃ > 5 × 10^-6) and low free-free intensity ( < 0.3 mK_CMB at 22 GHz). We have used this regression to remove the part linearly correlated with τ₃₅₃ from the 70 GHz transform at each scale. After this correction, we have applied the inverse transform to reconstruct the filtered map. We have then used the mr-detect tool with the same eight wavelet scales to detect and remove point sources. Only the first scale was used for source detection. The final map has been translated into units of specific intensity (Jy/sr) using the conversion coefficients³ for the 70 GHz band. We have applied the colour correction interpolated for the − 0.14 spectral index expected for free-free emission at 70 GHz.

Appendix C: H I spin temperature

The evolution of the maximum log-likelihood of the γ-ray fit in the 10^2.6−5 MeV band can be used to constrain the choice of average H ispin temperature that best matches the N_HI column density maps of the different clouds to the structure of the γ rays produced by CR interactions in the atomic gas. We find that the fit quality peaks near 400 K (Fig. C.1). We have not attempted

to use different spin temperatures for the different clouds present in the analysis region.

Appendix D: Best fit coefficients of the γ-ray and dust models

The best fit corresponds to an H_I spin temperature of 400 K. For each parameter we first give the statistical uncertainty from the fit and jackknife tests, then the standard deviation obtained by varying the H ispin temperature.

Appendix E: Previous γ-ray measurements of the X_CO conversion factor

All Tables

All Figures