© ESO 2018

1 Introduction

Large dust grains in the interstellar medium (ISM) absorb starlight from the UV to near infrared domain causing appreciable extinction and reddening in stellar observations (Draine 2003). These grains in thermal equilibrium with the ambient radiation field re-emit in the infrared to millimetre wavelengths. If the dust grains are well mixed with the interstellar gas, the reddening they cause and their thermal emission both provide indirect tracers of the gas column density.

The reddening caused by dust grains is generally estimated from the photometry of a set of stars in a given direction. A family of methods is based on the comparison of the observed stellar colours in various photometric bands with reference colours corresponding to zero extinction (NICE, NICER, NICEST algorithms: Lada et al. 1994; Lombardi & Alves 2001; Lombardi 2009). Such colour-excess methods have been applied to 2MASS photometry to derive two-dimensional (2D) maps of the dust reddening over the whole sky (Rowles & Froebrich 2009; Juvela & Montillaud 2016). Marshall et al. (2006) developed a method to estimate the dust reddening and distance simultaneously by comparing the observed colour from 2MASS data with the intrinsic colours and distances of stars provided by the Besancon Galaxy model (Robin et al. 2003). More recently, Green et al. (2015) used the combined information of Pan-STARRS 1 and 2MASS surveys to derive a three-dimensional (3D) reddening map in a fully probabilistic framework described in Green et al. (2014).

The thermal emission of the large grains has been recorded at 100 and 240 μm with IRAS and DIRBE and has been used to estimate dust optical depths over the whole sky (Schlegel et al. 1998). This map has been scaled to E(B − V) colour excesses using a set of quasar reddening measurements, but its structure reflects the distribution of optical depths. In a more recent study (Planck Collaboration XI 2014), the emission measured with Planck and IRAS between 353 and 3000 GHz (100 to 850 μm) has been used to model the spectral energy distributions (SEDs) of the thermal emission and to derive an all-sky map of the dust optical depth in emission at 353 GHz, τ₃₅₃. The latter relates to dust column density and grain emissivity. Assuming uniform dust properties (size distribution, chemical composition, and grain structure), the dust optical depth can be converted into extinction or reddening.

In addition to the information provided by dust tracers, the γ rays with energies above a few hundred million electron volts that are produced by the interaction of cosmic rays (CRs) with gas nucleons provide a measure of the total gas column density under the assumption of a uniform CR flux through a given cloud (Cesarsky & Volk 1978; Lebrun et al. 1982; Strong et al. 1988). As shown in Remy et al. (2017), the interstellar γ-ray spectra measured in the nearby clouds of the Galactic anti-centre region studied here are consistent with a uniform CR flux across the region and with a uniform penetration of the cosmic rays with energies above a few billion electron volts through the different gas phases of the clouds, from the atomic envelopes to the molecular cores. Therefore we can use γ rays as a robust tracer of the total gas column density, N_H, in this region.

Several studies have noted an increase in dust opacity (optical depth per gas nucleon, τ₃₅₃ ∕N_H) with increasing gas column density from the atomic to molecular phases. In the Taurus clouds, the opacities increase in the dense CO filaments 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 a study of the Galactic anti-centre region including the Taurus clouds (Remy et al. 2017), we have found opacity enhancements reaching a factor of six toward the CO clouds, even though we provided a measure of the additional H₂ gas that is not linearly traced by CO intensity, W_CO, in the directions where the ¹²CO line emission saturates (CO_sat). We also found that the opacity starts to increase in the dark neutral medium (DNM) at the HI –H₂ transition, that is, at lower gas densities than the few thousand per cubic centimetre sampled in the bulk of the CO-bright volume. The grain opacity found in the DNM compares well in the anti-centre clouds (Remy et al. 2017) and in the Chamaeleon complex (Planck and Fermi Collaborations Int. XXVIII 2015).

The opacity increase could be caused by an evolution of the grain properties changing their mass emission coefficient, as suggested by recent theoretical works (Köhler et al. 2015), or by spatial variations in the dust-to-gas mass ratio. The latter possibility also impacts the relation between dust reddening and gas column density.

This paper is a follow up study of our previous work (Paper I: Remy et al. 2017), which focused on the analysis of interstellar γ rays and dust thermal emission in the nearby anti-centre clouds. We aim here to study variations in the distribution of E(B − V)∕N_H ratios, across the different gas phases and in different types of clouds, using γ rays to trace the total gas column density. In order to directly compare environmental changes in E(B − V)∕N_H and in τ₃₅₃∕N_H, we use the same region, sampling, analysis method, and gas data as before.

The various ISM tracers used in this analysis and the expressions that relate E(B − V)∕N_H and τ₃₅₃ ∕N_H to dust properties are introduced in Sect. 2. The joint analysis of dust and γ-ray tracers is detailed in Remy et al. (2017) and summarised in Sects. 3.1 and 3.2. The E(B − V)∕N_H and τ₃₅₃ ∕N_H ratios are derived from fits of dust and γ-ray models to the data, as described in Sect. 3. The results of the joint “γ + E(B − V)” analysis aregiven in Sect. 4 and compared to the results of the former “γ + τ₃₅₃” analysis. In Sect. 5, we compare our results to the expectations of recent theoretical models of dust grain evolution (Jones et al. 2013; Köhler et al. 2015) and discuss how the variations in τ₃₅₃ ∕N_H and E(B − V)∕N_H relate to changes in the other physical quantities characterising large grains. The capability of dust tracers to probe the gas column density are discussed in Sect. 5.

Additionally, we compare the reddening maps that are currently available and comment on the choice made for this analysis in Appendix A. We present the best-fit coeffcients of the dust and γ-ray models in Appendix B.

2 Data

For our analyses, we have selected the same anti-centre region as that studied in Remy et al. (2017), extending from 139° to 191° in Galactic longitude and from − 56° to −3° in latitude. Regions with large gas column densities from the Galactic disc in the background and regions with a coarser angular resolution inthe HI map have been masked. The analysis region encompasses a variety of clouds, from diffuse and mostly atomic entities to dense molecular clouds forming stars at different rates, and spanning distances between 140 and 410 pc from the Sun (see Table 1 of Paper II: Remy et al. 2018).

2.1 Dust emission and extinction

The dust optical depth is defined as the product of the dust emission cross section per unit mass, κ_ν in cm² g⁻¹, and the dust mass column density, M_dust. The emission cross section can be expressed as κ_ν = with κ₀ being the value at the reference frequency ν₀ and β the spectralindex (Hildebrand 1983; Compiègne et al. 2011; Planck Collaboration XI 2014). The dust mass column density is expressed as a function of the hydrogen column density as M_dust = R_DGμ_HN_H with μ_H = 2.27 × 10⁻²⁷ kg the mean gas mass per hydrogen atom, and R_DG the dust-to-gas mass ratio. So the dust emission optical depth can be expressed as: (1)

and then we can define the dust emission opacity or hereafter in short, opacity (Planck Collaboration XI 2014), as: (2)

The dust extinction, A_λ, at a wavelength λ is proportional to the extinction optical depth along a line of sight, thus to the integral of the dust volume density, n_dust, and of the extinction cross section (Spitzer 1998): (3)

Assuming the grains can be described by spheres of effective radius a and mean mass density ρ_gr, we have with the efficiency factor for extinction. The extinction optical depth can therefore be expressed as a function of the hydrogen column density as: (4)

Substituting this relation in the expression for the visual extinction, A_V, and using the definition of the total-to-selective extinction factor yields an expression of the specific reddening: (5)

A relation between the optical depth and the reddening can be obtained by combining Eqs. (1) and (5): (6)

Provided we can reliably trace the total gas column densities in all gas forms, variations in τ₃₅₃ ∕N_H and in E(B − V)∕N_H can probe intrinsic variations in dust properties and in the dust-to-gas mass ratio. The knowledge of the gas column density and cloud type can inform us on environmental effects on dust evolution in the different gas phases.

We have used the dust optical depth at 353 GHz, τ₃₅₃, from Planck Collaboration XI (2014), which was derived by fitting the SED of the dust emission between 353 and 3000 GHz (100 to 850 μm) on a 5′ scale, as measured with Planck and IRAS.

The stellar reddening map used in the following comes from the 3D E(B − V) map of Green et al. (2015), integrated up to its maximal distance. This map was derived using Pan-STARRS 1 and 2MASS stellar data. We have tested other reddening maps to quantify the differences induced by the use of other stellar data and different methods. As the reddening map of Green et al. (2015) presents a better correlation with both the τ₃₅₃ optical depth and the gas column densities inferred in γ rays, we use this map to study the relation between E(B − V) and N_H (see the discussion on the choice of extinction map in Appendix A).

The distributions of optical depth τ₃₅₃, reddening E(B − V), and of their E(B − V)∕τ₃₅₃ ratio are displayed for the region of interest in Fig. 1. Regions of intermediate E(B − V) (~0.3–1 mag) are shown tohave a ratio close to the value of 1.49 × 10⁴ mag found by Planck Collaboration XI (2014), yet the ratio drops by a factor up to three at large extinctions in the direction of compact molecular clumps and filaments (further discussed in Sect. 5).

Fig. 1 Maps of the dust optical depth at 353 GHz (top left panel), of the stellar reddening by dust (top right panel), and of their ratio (bottom left panel). The map on the lower right shows in green the pixels with an angular resolution of 0. °11 in E(B − V), which have been retained for mask2, and in black the pixels with a lower resolution (see Sect. 3.3).

2.2 γ-ray data

We used the same γ-ray data as in Remy et al. (2017), namely six years of Pass 8 photon data, recorded between 0.4 and 100 GeV by the Fermi Large Area Telescope (LAT; Atwood et al. 2009), and the associated instrument response functions (P8R2_CLEAN_V6). In order to reduce the contamination of residual CR tracks in the photon data and by Earth atmospheric γ rays, we 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°; for details, see Nolan et al. 2012). To account for the spill-over of emission produced outside the analysis region, but reconstructed inside it, we modelled point sources and interstellar contributions in a region 4° wider than the analysis region.

The positions and 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 interstellar radiation field (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. 2012). Details on the isotropic spectrum for the extragalactic and residual charged-particle backgrounds are given in Atwood et al. (2013) and Acero et al. (2015).

2.3 Gaseous components

We have used the same data as in Remy et al. (2017) to describe the gaseous components. We will briefly summarize their properties in the following, and more details can be found in the aforementioned paper.

The bright atomic gas is traced by the 21 cm line emission of hydrogen. 74% of the analysis region is covered by the GALFA-HI survey at 4′-resolution (Peek et al. 2011). Elsewhere we have used the EBHIS survey at a resolution of 10.8′ (Winkel et al. 2016). Both surveys were re-sampled into the 0. °125-spaced Cartesian grid used for the analyses. We have verified the tight correlation in column densities between the two surveys where both are available in our analysis region (Remy et al. 2017).

In order to trace the bright 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.

HI and CO lines have been used to kinematically separate cloud complexes along the lines of sight and to separate the nearby clouds from the Galactic backgrounds. We have used the narrow-band GALFA data-cubes with their original velocity resolution of 0.18 km s⁻¹ in the local standard of rest. The EBHIS frequency sampling is coarser, with a velocity resolution of 1.44 km s⁻¹, but still adequate to separate clouds in velocity. Details of the pseudo-Voigt line decompositions and selection of longitude, latitude, and velocity boundaries of the different complexes are given in Remy et al. (2017). The separation resulted in six different complexes for both HI and CO components: Cetus, Main Taurus, South Taurus, North Taurus, Perseus, and California. They are shown in CO intensity, W_CO, in Fig. 2 and in N_HI column density in Fig. 1 of Remy et al. (2017).

Ionisedgas in the California nebula (around l =160. °1 and b = −12. °4), in the G159.6-18.5 HII region, and along the Eridanus loop is visible in Hα emission (Finkbeiner 2003) and in free-free emission at millimetre wavelengths (Bennett et al. 2013). We used the Planck LFI data at 70 GHz (release 2.01 data, Planck Collaboration I 2016) with an angular resolution of 14′ to measure the free-free emission in the analysis region. To do so, from the LFI data we successively removed the contributions from the Cosmic Microwave Background, from dust emission, and from point sources (see Remy et al. 2017, for details).

Fig. 2 CO intensity maps of the Cetus (Cet), Main Taurus (TauM), South Taurus (TauS), North Taurus (TauN), Perseus (Per), and California (Cal) clouds, and of the Galactic disc background (Gal). The solid curve and dashed line show the low-latitude borders of two analysis regions, mask 0 and mask 1, respectively.

3 Method and analyses

3.1 Dust and γ-ray models

The large dust grains responsible for the stellar reddening at near-infrared (NIR) and optical wavelengths are taken to be well mixed with the interstellar gas, so we have modelled the spatial distributions of the E(B − V) reddening as a linear combination of gas column densities in the various phases and from the different cloud complexes. The model is expressed similarly to Eq. (6) of Remy et al. (2017), substituting τ₃₅₃ for E(B − V): (7)

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 in the seven complexes. and stand for the column densities of the dark gas in the DNM and CO_sat components. They are deduced from the coupled analyses of the γ-ray and dust data (see Sect. 3.2). The y model parameters have been estimated using a χ² minimisation.

At the energies relevant for the LAT observations, the particle diffusion lengths in the interstellar medium exceed the cloud dimensions and there is no spectral indication of variations in CR density inside the clouds (Remy et al. 2017). The emissivity spectrum of the gas follows the average one obtained in the local ISM, q_LIS, given by Casandjian (2015). The interstellar γ radiation can therefore be modelled, to first order, as a linear combination of gas column densities in the various phases and different regions seen along the lines of sight. The model also includes a contribution from the large-scale Galactic inverse-Compton emission, I_IC (l, b, E), an isotropicintensity to account for the charged-particle and extragalactic backgrounds, I_iso (l, b, E), point sources, S_j(E) and contamination of sources present in the peripheral region outside the analysis region merged in a single map, S_ext(l, b, E). The γ-ray intensity I(l, b, E), expressed in cm⁻² s⁻¹ sr⁻¹ MeV⁻¹, can be modelled as : (8)

with and being the reddening maps in the DNM and CO_sat phases, respectively, extracted from the coupled dust and γ-ray analyses (see Sect. 3.2).

To ensure photon statistics robust enough to follow details in the spatial distributions of the different interstellar components, we analysed the data in five broad and independent energy bands, bounded by 10^2.6, 10^2.8, 10^3.2, 10^3.6, and 10⁵ MeV and in the total energy band. We left free normalizations (the q factors) to account for possible deviations in CR density and spectrum for the interstellar components and to account for possible changes in source spectra since we have used a longer time interval of LAT data than for the 3FGL catalogue. Cloud-to-cloud variations in CR flux are monitored by the q_HI,i scale factorsand will affect the normalizations equally in all energy bands, whereas a change in CR penetration in a specific cloud will be revealed 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(E) spectrum and the best-fit q_HI,i(E) normalizations. 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.

Figure 3 shows the distribution of γ rays of interstellar origin in the region, together with the stellar reddening displayed at the angular resolution of the Fermi LAT. The γ rays of interstellar origin are obtainedthrough γ-ray data analysis after subtraction of the inverse Compton, isotropic, and discrete source contributions (see Eq. (8)). Figure 3 illustrates striking similarities between the spatial distributions of the stellar reddening by dust and of the gas traced by CR interactions. The same likeness had been found between the dust optical depth τ₃₅₃ and the γ rays (see Fig. 3 of Remy et al. 2017). The spatial correlation supports the use of dust tracers to complement the HI and CO data in modelling the γ rays of gaseous origin. Yet, as for the optical depth τ₃₅₃, small differences between the two maps in Fig. 3 readily point to local variations in specific reddening E(B − V)∕N_H.

Fig. 3 Top: γ-ray counts of gaseous origin recorded in the 0.4–100 GeV energy band in a 0. °125 pixel grid. γ-ray emissions other than due to cosmic-ray interactions in the gas have been subtracted. The map has been smoothed with a Gaussian kernel of 0. °14 dispersion for display. Bottom: stellar reddening by dust displayed at the Fermi-LAT angular resolution for comparison.

3.2 Mapping DNM and COsat components

Both the stellar reddening and interstellar γ-ray emission trace the total gas to first order. They are modelled as linear combinations of several gaseous components including bright (i.e. detected in radio or millimeter emission lines) neutral gas in the atomic and molecular phases (HI, CO emissions), ionised gas (free-free emission), and dark gas. The analysis couples the dust and γ-ray models to constrain the structure of the additional gas not seen in HI, CO and free-free emissions, namely the gas in the DNM at the HI–H₂ transition and the CO_sat H₂ gas missed by the saturation of the ¹²CO lines at large optical thickness in the dense cores of molecular clouds. Both components show up as positive residuals over the expectations from the HI-bright, CO-bright, and free-free-bright gas components. To extract them, we 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 kept only the positive residuals above the noise and separated the DNM and CO_sat components in directions outside and inside ofthe 7 K km s⁻¹ contour in W_CO intensity.

The DNM and CO_sat templates of additionalgas used in the γ-ray model are inferred from the dust model and conversely. We cycled between the dust and γ-ray fits iteratively until the quality (log-likelihood) of the fit to the γ-ray data saturated. The spatial correlation between the γ-ray and dust excesses ensured by this iteration is essential to extract genuine gas quantities and discard local variations in CR density or in dust properties. The hydrogen column densities in the DNM and CO_sat gas phases were estimated using the γ-ray emissivity per hydrogen atom in the HI phase and assuming that the DNM and CO_sat phases are pervaded by the same average CR flux as in the atomic envelope. This assumption is validated by the uniformity of the γ-ray emissivity spectrum (i.e. CR spectrum) found in the different gas phases, including the dense CO and CO_sat phases. By comparing the results of the γ-ray analyses using τ₃₅₃ and E(B − V) as dust input, we can check for any small bias in the determination of hydrogen column densities in γ rays, which would be indirectly induced by the dust data used to infer the DNM and CO_sat distributions.

3.3 Choice of analysis masks

The E(B − V) map was constructed with a resolution of 0. °11 or 0. °23 depending on the number of stars available in each pixel. Directions with a resolution below 0. °11 are shown in black in Fig. 1. They correspond to low latitudes far from the densely populated regions near the Galactic plane (b < −30°) and to directions toward highly obscured molecular clouds. In the very diffuse environments below –30° in latitude where E(B − V) < 0.3 mag, the fluctuations seen in the E(B − V)∕τ₃₅₃ ratio result from uncertainties in the zero level in τ₃₅₃, fluctuations from the cosmic infrared background in the τ₃₅₃ map (Planck Collaboration XI 2014), and from an increased level of noise in E(B − V) at low stellar surface densities (Green et al. 2015). In order to compare dust emission and reddening on the same scale and to test the impact of the change in angular resolution across the E(B − V) map, we performed the analysis with three different masks:

• mask0, including the full region displayed in Fig. 3 (as for the τ₃₅₃ and γ-ray analysis of Remy et al. 2017);

• mask1, excluding latitudes below −30° to keep regions with high signal-to-noise in E(B − V) (dotted cut in Fig. 2);

• mask2, excluding all the pixels with an angular resolution in E(B − V) worse than 0. °11 (displayed in black in Fig. 1).

The California, Perseus, and Taurus molecular clouds are fully kept in mask1, but Cetus is masked out in both mask1 and mask2. Those masks discard 30% and 18% of the pixels with CO emission in the South and North Taurus clouds, respectively. Fifteen percent of the pixels toward the Perseus cloud are masked out with mask2. Since the CO-to-H₂ conversion factor (X_CO = /W_CO) tends to decrease from the edge to the core of molecular clouds (Bell et al. 2006; Bertram et al. 2016; Remy et al. 2017) masking out parts of a cloud changes the mean X_CO value of the cloud compared to the result obtained for the complete cloud (mask0). This effect is discussed in Sect. 4.3.

4 Results

4.1 Gas column density

The coupledanalysis of γ rays and dust, with the latter in emission or reddening, takes advantage of the common underlying gas distribution in the two tracers to constrain the amount and the structure of the gas ill-traced by the radio line surveys and free-free emission (i.e. DNM and CO_sat). This gas is jointly revealed by its dust and CR content. We first investigate how the choice of dust tracer affects the derivation of gas column densities in the DNM and CO_sat phases, thenin the total gas.

Figure 4 shows that the hydrogen column densities derived in the DNM from the “γ + E(B − V)” and “γ + τ₃₅₃” analyses exhibit very similar spatial structures and Fig. 5 illustrates that the column density values compare well across the whole range. The maps present only small differences; for instance, slightly more extended diffuse cloud edges below 4 × 10²⁰ cm⁻² in the “γ + τ₃₅₃” map, and larger column densities, probably from the Galactic disc background, above 1.5 × 10²¹ cm⁻² near the Galactic plane (b > −8°) in the “γ + E(B − V)” map. But the total mass in the “γ + τ₃₅₃” map is only 9% lower than in the “γ + E(B − V)” one.

Similarly, Fig. 6 illustrates that the two derivations of the CO_sat component present the same overall spatial structure. Yet, local differences in column densities are more accentuated in the CO_sat component than in the DNM (see Fig. 5). The main differences arise at intermediate column densities of (1− 3) × 10²¹ cm⁻² where the “γ + E(B − V)” map yields larger values in several structures. Hence, the total CO_sat mass in the “γ + E(B − V)” map is 22% larger than in the “γ + τ₃₅₃” one. Conversely, the densest filaments and clumps emitting in the rarer ¹³CO isotope are enhanced in the “γ + τ₃₅₃” map compared to the “γ + E(B − V)” one (beyond 6 × 10²¹ cm⁻²). This difference results from the combination of the sharp rise in τ₃₅₃∕N_H at the large column densities characterising these cores (see Fig. 16 of Remy et al. 2017) and from the loss of contrast in reddening because of the lower angular resolution of the E(B − V) map due to the loss of stars in those heavily obscured directions (see the grey contours in Fig. 6).

The γ-ray counts of gaseous origin are derived from the total data after subtraction of non-gaseous components estimated from the best-fit model (IC emission, isotropic backgrounds, and point sources). Those counts are displayed in Fig. 3. They are converted to total gas column density, N_{H γ}, using the average emissivity spectrum per gas nucleon measured in the atomic phase in the region. Figure 7 compares the gas column densities obtained in the “γ + E(B − V)” and “γ + τ₃₅₃” analyses by showing their ratio. The subtraction of point sources results in very localized differences, by less than ± 30%, between the two gas maps. Only 10% of the point sources present in the region have seen their flux slightly change between the two analyses, but the differences remain compatible with the individual normalization uncertainties given in the 3FGL catalogue (which used a different interstellar background model).

Away from those few point sources, the two column density maps agree within ± 3%, except toward the most diffuse regions, at b < −45° or at l <155° and b < −10°, where differences of about ± 10% result from the noise level in the reddening map. Hence, despite the differences discussed in Figs. 4 and 6, the N_{H γ} column densities are not significantly affected by the choice of dust tracer used to infer the DNM and CO_sat components that complement the HI and CO data in the γ-ray model in order to constrain the gaseous components. The robustness of the N_{H γ} map is due to the fact that most of the gas mass resides in the HI -bright and CO-bright parts (which are not inferred from fits to the dust) and that the large DNM mass corresponds to intermediate column densities where the two dust tracers yield comparable structures.

Figure 7 also compares the quality of the γ-ray fits obtained in the “γ + E(B − V)” and “γ + τ₃₅₃” analyses by showing the residuals (photon data minus model, in sigma units) in the overall energy band for each analysis. In both cases, the linear model provides an excellent description of the observed data in the 0.4–100 GeV band, as well as in the individual energy bands that are not shown here. The residuals remain largely consistent with noise fluctuations below ± 1σ. The small modelexcess noted towards the brightest CO parts of Main Taurus, Perseus, and California clouds in the “γ + τ₃₅₃” analysis (Remy et al. 2017) disappears when using dust reddening because of the reduced contrast in the E(B − V) map mentioned above. This confirms that the model excess in “γ + τ₃₅₃” was caused by the marked rise in τ₃₅₃∕N_H opacity that characterises the cold population of grains present in the densest CO cores (see Fig. 17 of Remy et al. 2017), rather than by a loss of CRs.

Fig. 4 Hydrogen column density maps obtained for the DNM with the “γ + E(B − V)” (top panel) and “γ + τ353” (bottom panel) analyses of the full region (mask0).

Fig. 5 Number distributions of hydrogen column densities in the DNM (black) and COsat (grey) components, obtained from the “γ + τ353” (dashed) and “γ + E(B − V)” (solid) analyses of the full region (mask0). The numbers of pixels for each component have been normalized by the maximum value found in the “γ + E(B − V)” analysis (solid curves).

Fig. 6 Hydrogen column density map of additional gas in the 12CO-saturated parts of the Taurus and Perseus clouds, derived from the “γ + τ353” (top panel) and “γ + E(B − V)” (bottom panel) analyses of the full region (mask0). The pixels within the grey contours have a lower angular resolution in E(B − V) (mask2). The magenta contours enclose regions with 13CO intensities above 2 K km s−1 (Ridge et al. 2006; Narayanan et al. 2008).

Fig. 7 Left panel: map of the ratios of hydrogen column densities derived from interstellar γ rays in the “γ + E(B − V)” and “γ + τ353” analyses, thewhite contours marking the positions of γ-ray point sources. Middle and right panels: residual γ-ray maps (photon counts minus best-fit model, in sigma units) in the 0.4–100 GeV band for the “γ + E(B − V)” (middle) and “γ + τ353” (right) analyses, for an HI spin temperature of 400 K. The black contours outline the CO clouds at the 1 K km s−1 level.

4.2 Specific reddening and opacity variations

Our analyses yield average dust properties per H atom in each gas phase for each cloud complex. The best-fit coefficients of the N_HI maps in the dust model directly give the average specific reddening E(B − V)∕N_H (alternatively τ₃₅₃∕N_H opacity) in the atomic envelopes of each cloud complex. We also used the best-fit coefficients of the W_CO maps and the γ-ray estimates of the X_CO factors to derive the average specific reddening (or opacity) in the CO-bright phase of each complex. The results are listed in Table 1 together with the opacities found in the “γ + τ₃₅₃” analysis (Remy et al. 2017).

We selected different masks to investigate a possible bias due to the variable resolution of E(B − V) map and we find only up to 10% variation in the average E(B − V)∕N_H when excluding the pixels of lower resolution. On average this effect is not significant, but it could be more important locally within the unresolved compact filaments of molecular clouds.

The average cloud opacities were found to systematically increase by 30% to 100% from the HI to the CO phase. The opacity maps showed a gradual and coherent increase which was related to grain evolution from the atomic to the molecular phase (see Figs. 16 and 17 of Remy et al. 2017, and the references discussed in the following section). The present results also show a systematic increase by 30% to 60% in average specific reddening between the HI and CO phases of the clouds, except for Perseus. A large fraction of this cloud is more coarsely mapped in E(B − V), but the results found with mask2 also yield consistent specific reddenings in the HI and CO phases of Perseus. However, this cloud isvery compact, both in HI and CO (see Fig. 1 of Remy et al. 2017), and it lies behind North Taurus, meaning that the separation between the HI and CO contributions to the total reddening is less efficient in this cloud than in the others.

The specific reddening E(B − V)∕N_H values found in the atomic parts of the North Taurus and Cetus clouds are close to the reference value of 1.7 × 10⁻²² mag cm² of Bohlin et al. (1978) that applies to lines of sight with E(B − V) < 0.5 mag in the solar neighbourhood and in the Galactic disc. They also agree with the value of (1.71 ± 0.07) × 10⁻²² mag cm² found in the atomic gas of the Taurus region from the visual extinction map of Dobashi et al. (2013; Paradis et al. 2012). We find, however, larger values in the range of (1.9−2.1) × 10⁻²² mag cm² in the atomic gas of the Main Taurus, South Taurus, and California clouds. Paradis et al. (2012) had noted variations of the same magnitude between local clouds in their extinction data based on 2MASS stellar photometry.

The average specific reddening reaches (2.09 ± 0.15) × 10⁻²² mag cm² in the denserDNM at the HI -H₂ interface (given by the best-fit parameters of the γ-ray model ). FUSE observations probed the same type of medium in diffuse H₂ clouds up to 1 mag in reddening (Rachford et al. 2009). The specific reddening changed between (1.68 ± 0.10) × 10⁻²² mag cm² and (2.38 ± 0.96) × 10⁻²² mag cm² depending on the inclusion or not in the fit of the data points at small reddening. The possible change in slope with reddening hinted at possible changes in dust properties with higher extinction. The average value that we find in the DNM is consistent with the FUSE observations and the higher values we get in the range of (2.2−2.7) × 10⁻²² mag cm² in the CO-bright clouds further support an evolution of the dust extinction properties.

We further note that the present E(B − V)∕N_H ratios are systematically 1.2 to 2.2 times lower than the values obtained in parts of the South Taurus, Main Taurus, and Perseus clouds by Chen et al. (2015) from their visual extinction map based on 2MASS and XSTPS-GAC stellar data. Further analyses of common regions using the same N_{H γ} data and direct comparisons of the E(B − V) and A_V maps are necessary to find the origin of the discrepancy.

The good correlation between the hydrogen column density derived from γ rays and the dust reddening is presented in Fig. 8 together with the evolution of the E(B − V) /N_{H γ} ratio as a function of the gas column density. The reddening map has been convolved with the appropriate LAT point-spread function and the angular resolution of both maps has been degraded to 0. °375 to improve the photon statistics in the histogram bins. At this resolution, the dust reddening and gas column density linearly correlate across the whole gas range. We detect no upward inflexion in the E(B − V)∕N_H ratio with increasing gas column density, contrary to the systematic rise seen in dust opacity across the same gas range (see Fig. 15 of Remy et al. 2017).

The spatial distribution of the specific reddening is shown in Fig. 9. We also map the variations relative to the average value of (2.02 ± 0.48) × 10⁻²² mag cm² found for the whole field (mask0) across all gas phases. The E(B − V) /N_{H γ} values increase in the DNM and molecular phases compared to the diffuse regions. The 20% drop towards Perseus is due to a relative lack of reddening, possibly due to the presence of young stellar objects rather than to an excess of gas column density. Large ratios are found in the region at 155° < l < 163° and −33° < b < −30° where the (independent) thermal emission of the grains has a spectral index greater than 1.8 (see Fig. 2 of Remy et al. 2017).

Across the whole region, the E(B − V) /N_{H γ} ratios vary typically by ±30% around the average value. This is far lower than the three-to-six fold variations found in the τ₃₅₃ ∕N_H ratios towards the same clouds, at the same angular resolution, and with the same analysis method (Remy et al. 2017). The implicationsof this result on dust grain properties will be discussed in Sect. 5.

Fig. 8 Upper panel: 2D histogram of the correlation between the total gas column density, NH γ, measured with the 0.4–100 GeV interstellar γ rays, and thedust reddening, , convolved with the LAT point-spread-function for an interstellar spectrum. The maps are sampled at a 0. °375 resolution. The solid line marks the best-fit slope found across the whole field (mask0); it corresponds to a /NH γ ratio of (2.02 ± 0.48) ×10−22 mag cm2 averaged over all gas phases. Lower panel: evolution of the /NH γ ratio in equally spaced bins in NH γ. The black dashed, white dashed, and black solid curves give the means obtained with mask0, mask1, and mask2, respectively.The shaded area gives the standard deviation in each bin using mask2 and the error bars give the standard error on the means.

4.3 XCO factors

The best-fit coefficients of the γ-ray and dust models yieldindependent estimates of the X_CO factors under the respective assumptions of a uniform CR flux or uniform dust properties per gas nucleon in the HI and CO phases of a given cloud complex. The results obtained from the “γ + E(B − V)” and “γ + τ₃₅₃” analyses for different clouds and with different spatial masks are presented in Table 2. The results obtained for the whole region (mask0) are displayed in Fig. 10 as a function of the mean dust optical depth, , in the cloud complex.

Remy et al. (2017) reported significant variations in X_CO when averaged over clouds that are structured differently. The X_CO factors tend to decrease from diffuse to more compact clouds, as expected from the models of the formation and photodissociation of H₂ and CO molecules. The models predict strong inward gradients in the relative abundances of CO and H₂, which result in a marked decline in X_CO from the diffuse molecular envelopes to the dense and well-shielded cores (Bell et al. 2006; Bertram et al. 2016). Therefore, the X_CO average over a cloud should vary with the surface fraction of dense versus diffuse regions and with the degree of saturation of the CO lines in the densest parts. Sorting the clouds in Fig. 10 by their mean value partially captures their diffuse (low ) or compact (high ) character. The systematic decline seen in X_CO with all gas tracers clearly confirms the presence of CO-to-H₂ abundance gradients inside the clouds. Table 2 further shows that masking out large fractions of the South Taurus and North Taurus complexes induces significant changes in their average X_CO factors in all the analyses, thereby confirming the existence of spatial variations in X_CO within a cloud complex because of the different states of its constituents. The Main Taurus and California cloud complexes are only marginally affected by the change in mask. This is why their average X_CO factors remain stable for the different masks.

The γ-ray values of the X_CO factor in the Perseus cloud should be taken with care because of the significant level of cross-correlation between the compact HI and CO phases of the cloud when convolved at the angular resolution of the LAT, especially at low energies (see Sect. 6.1 of Remy et al. 2017). This is even more true when further lowering the dynamic ranges by removing pixels in mask2. The X_CO factors in Perseus vary significantly when masking out the densest gas pixels: the τ₃₅₃ -inferred value decreases by a factor of 2.3 because the masked pixels correspond to regions where the grains have particularly large emission opacities; the change is still pronounced in γ rays (although reduced to a factor of 1.6) in the “γ + τ₃₅₃” analysis because the masked pixels have large column densities in the τ₃₅₃ -inferred CO_sat component; the change is reduced to 20% in the γ-ray fit coupled to the dust reddening. Based on the γ-ray and reddening results, the average X_CO factor in Perseus should be close to (0.5–0.6) × 10²⁰ cm⁻² K⁻¹ km⁻¹ s.

By construction of the models, the DNM gas component has negligible influence on the X_CO estimates (Planck and Fermi Collaborations Int. XXVIII 2015). Figure 10 further shows that coupling the γ-ray analysis to different dust tracers to infer the CO_sat components does not change the γ-ray values of X_CO significantly; they remain stable in the Cetus, Main Taurus, and California clouds and decrease by 20% in the North and South Taurus clouds when replacing dust emission by reddening. The figure also illustrates that the reddening-based measurements are systematically closer to the γ-ray ones than the measurements based on dust emission, but they still differ by 15% to 60%. Remy et al. (2017) noted that the γ-ray and τ₃₅₃ estimates of X_CO differed more at larger because compact clouds exhibit large gas column densities where the dust grains have particularly large opacities, but we do not find a similar divergence between the γ-ray and reddening estimates of X_CO. The reddening-based factors are systematically smaller than those obtained with dust in emission (see Fig. 10), meaning that the results corroborate the existence of a strong opacity bias in X_CO measurements when they are based on dust emission only.

5 Discussion and conclusions

Several studies have reported variations in dust opacity (e.g. Stepnik et al. 2003; Flagey et al. 2009; Martin et al. 2012; Ysard et al. 2013; Roy et al. 2013; Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014; Köhler et al. 2015; Remy et al. 2017). The extinction curves also respond sensitively to local conditions and contain unique information about the grains along the lines of sight (Fitzpatrick & Massa 2007). Changes in the extinction law have also been reported (Cardelli et al. 1989; Foster et al. 2013; Schlafly et al. 2016). In both cases, the changes have been attributed to the evolution of grain properties under the local conditions. Irradiation can modify their chemical composition and structure while mantle accretion, coagulation, and ice coating change their structure. In Sect. 2.1, we reiterated how the τ₃₅₃ ∕N_H opacity and the E(B − V)∕N_H specific reddening respond to several grain properties in emission and in extinction. The results presented here and by Remy et al. (2017) show that both the emission and extinction characteristics vary significantly across the anti-centre clouds, but not in strict spatial correlation and not with the same amplitude. In the following, we discuss how the variations relate with the spectral characteristics of the dust emission and how they compare with predictions from recent theoretical models (Jones et al. 2013; Köhler et al. 2015).

Figure 11 displays how the dust characteristics vary with the colour temperature, T_dust, and spectral index, β, of the thermal emission of the large grains in the 353–3000 GHz band (Planck Collaboration XI 2014). The clouds have been sampled at an angular resolution of0. °375 which corresponds to a linear scale of 1 to 3 pc in the different complexes. The upper plots map the changes in specific reddening and in opacity in absolute values, while the third plot shows the relative changes in the E(B − V)∕τ₃₅₃ ratio to illustrate how the extinction and emission properties evolve relative to each other. To show the amplitude of the relative variations, we divided the E(B − V)∕τ₃₅₃ ratio by the mean slope of 1.49 × 10⁴ mag found in the correlation between τ₃₅₃ optical depths and quasar colour excesses across the sky, mostly through atomic clouds (Planck Collaboration XI 2014). The bottom plot shows the mean gas column density inferred in γ rays across the range of dust spectra in order to link the dust properties with their gaseous environment, from the diffuse clouds where the grains are immersed in low-density and mostly atomic gas and where they are well exposed to the ISRF (at low β index and large T_dust temperature), up to the dense, molecular, and obscured media in our sample.

Figure 11 highlights that the dust properties per gas nucleon in emission and in extinction vary coherently across the ISM phases, but that the gradients in specific reddening and in opacity are not commensurate. The rise in opacity, is more than twice as large as the coincident rise in specific reddening, as the gas becomes denser. We discuss these evolutions in the diffuse and dense media in turn.

5.1 Dust evolution in the diffuse medium

In the “γ + τ₃₅₃” analysis of the nearby anti-centre clouds (Remy et al. 2017), we showed that the emission opacity τ₃₅₃ ∕N_H increases by a factor of two in the HI and DNM phases above the diffuse-ISM value of (7.0 ± 2.0) ×10⁻²⁷ cm² H⁻¹ found in tenuous high-latitude cirruses (Planck Collaboration Int. XVII 2014; Planck Collaboration XI 2014). Comparable variations have been noted in the atomic gas across the sky, with values ranging from 6.6 to about 11 ×10⁻²⁷ cm² H⁻¹ (Planck Collaboration XI 2014). Using the Planck optical depth at 353 GHz and N_HI column densities from GALFA, Reach et al. (2015) also reported up to three-fold variations in the diffuse ISM.

In the anti-centre region, we find that the specific reddening increases by less than 50% from the very diffuse gas to the denser HI and DNM parts of the clouds. The values remain close to the reference of 1.7 × 10²² cm² mag of Bohlin et al. (1978) towards the sight lines intercepting atomic gas. They increase to an average of (2.09 ± 0.15) ×10⁻²² mag cm² in the DNM, but we find no clear trend with N_H across the HI and DNM phases (see Table 1 and Fig. 9), probably because they span the same range of column densities even though the UV field and the gas chemistry differ between the two media. The E(B − V)∕N_H values seen below 1.5 × 10²² cm² mag in Figs. 11 and 9 correspond to very slightly reddened regions (< 0.08 mag) where the difficult reddening measurements should be taken with care. Hence, the amplitude of the opacity rise in the HI and DNM phases significantly contrasts with the small variations in E(B − V)∕N_H ratio.

Figure 11 shows that the increase in emission opacity and in specific reddening both relate to an increase in β index and a decrease in colour temperature. Recent theoretical studies (e.g. Ysard et al. 2015, 2016; Köhler et al. 2015) suggest that the spectral changes in β relate to gas density-dependent processes and that the variations in opacity, β index, and colour temperature are closely connected. According to Jones et al. (2017), the carbon and silicate grains in the diffuse ISM have H-poor, aromatic-rich carbonaceous mantles. Their optical properties correspond to the low β, large temperature (≳ 20 K), and low opacity part in Fig. 11b. By accretion of C and H atoms from the gas phase, an additional H-rich amorphous carbon mantle forms on the surface of grains. This mantle is aromatised in regions well exposed to UV radiation, and more aliphatic in attenuated regions. The change in thickness and composition of the mantle induces a large increase by typically 0.3 in β index, as well as a doubling of the emission opacity (emissivity), causing the grain equilibrium temperature to drop by ~ 2 K. The changes in opacities and temperatures found in the anti-centre clouds suggest that the accretion and evolution of the additional mantle would occur primarily in the dense HI and DNM phases, which correspond in Fig. 11d to the green-to-orange zones and to the red zone (at β < 1.65 and 17 < T_dust < 18.5 K). The amplitudes of the variations in τ₃₅₃∕N_H and T_dust are compatible with the model predictions, but for a smaller change in β index in the observations.

At UV and optical wavelengths, the extinction efficiency factor (defined as the sum of the absorption and scattering efficiency) tends to the limit Q^ext = Q^abs + Q^sca = 2 for grain sizes larger than 0.1 μm (Bohren & Huffman 1983; Desert et al. 1990). This condition is verified in the model of Köhler et al. (2015) as they find Q^abs ~ Q^sca ~ 1 up to 0.6 μm in wavelength for all the grains larger than 0.11 μm and the aggregates larger than 0.17 μm. Mantle accretion does not significantly change the size of the grains; it induces only a 10% increase in R_V, from 3.5 ± 0.2 to 3.9 ± 0.2, and a less than 20% increase in extinction efficiency at B and V wavelengths (Köhler et al. 2015). Such constraints are compatible with the small dispersion in E(B − V)∕N_H found in the atomic and DNM phases of the anti-centre clouds. Deviations by ± 30% suggest that the population of large grains is rather stable in size distribution, absorption cross section, and in proportion to the gas mass across those phases (see Eq. (5)). Thus the three-to-six fold increase in emission opacity (Eq. (2)) would primarily be due to a change in grain emissivity κ_ν caused by mantle accretion in the dense HI and diffuse H₂ media, rather than to an increase in dust-to-gas mass ratio, R_DG.

Schlafly et al. (2016) selected 37 000 stars with stellar parameters estimated from the APOGEE spectroscopic survey and used Pan-STARRS 1 photometry to constrain the extinction curve from optical to infrared wavelengths and to derive the R_V factor. Up to ~2 mag in reddening, their values of R_V vary in the 3.1–3.5 range around a mean of about 3.3 which is lower than the model prediction. The observed factors do not systematically increase with reddening in this range. The amplitude of these R_V variations is consistent with the dispersion we measure in E(B − V)∕N_H in the diffuse anti-centre clouds. Schlafly et al. (2016) further note that larger R_V values relate to lower β indices, at variance with the model of Köhler et al. (2015), but compatible with the decrease in E(B − V)∕N_H and E(B − V)∕τ₃₅₃ we see in Fig. 11 towards the lowest β indices and largest temperatures. In the proximity of UV sources or at low densities where there is less screening against UV radiation, the grains are more exposed to the ISRF which causes their darkening and an increase in R_V with time (Cecchi-Pestellini et al. 2010). The E(B − V)∕τ₃₅₃ ratio is also expected to decrease as the ambient radiation field hardens (Ysard et al. 2015). Figure 16 of Remy et al. (2017) shows how the dust specific power (radiance per gas nucleon) varies across the anti-centre region. For grains in thermal equilibrium, the specific power traces variations in the heating rate. Away from localized hot spots associated with HII regions, the relation between the specific power and specific reddening is unclear. Their variations in the HI phase correlate at small spatial scales in places and not in others.

Fig. 9 Spatial variations of the amount of reddening per gas nucleon, /NH γ, with the total gas measured by NH γ at 0. °375 resolution (top panel) and its variations relative to the average value of (2.02 ± 0.48) × 10−22 mag cm2 (bottom panel). Quantities with a tilde are convolved with the LAT response for an interstellar spectrum. The black contours outline the shape of the CO clouds at the level of 1 K km s−1 in WCO intensity.

5.2 Dust evolution in the molecular regions

Dust evolution continues in the denser ISM phases. The very small carbonaceous grains rapidly coagulate onto the large ones at gas volume densities above a few 10³ cm⁻³ that are typical of CO-bright media (Jones et al. 2017). Köhler et al. (2015) predict only a small impact on the β index in the 353–3000 GHz band, which should increase by less than 0.1 whether or not the final grains become coated with ice. The structural change and ~50% increase in emission opacity before and after coagulation should lower the grain temperatures below 18 K if the aggregates are ice free. In the case of an ice mantle, the grains radiate much more efficiently (by a factor of 3.5) and they cool below 17 K. The overall change in opacity between the diffuse atomic and dense molecular media is therefore predicted to be of the order of 3 without ice coating and 7 with ice (Köhler et al. 2015).

Opacity values that are 2 to 3 times above the diffuse-ISM average have been found in the local ISM (Flagey et al. 2009; Planck Collaboration XXV 2011; Ysard et al. 2013; Stepnik et al. 2003). Steeper rises have been found in nearby clouds, with factors of 3 to 6 in the Taurus and Perseus clouds (Remy et al. 2017), and 2 to 4.6 in the Chamaeleon clouds (Planck and Fermi Collaborations Int. XXVIII 2015). Variations of this magnitude agree with the above predictions. Yet, the opacity increase in the more massive California molecular cloud is less than two-fold (Remy et al. 2017). Whether this disparity could be due to the loss of small-scale contrast at larger distance or to a different dust state remains to be verified. Large opacities are associated in Fig. 11 with spectral indices β > 1.65 and colour temperatures below 17 K in reasonable agreement with the model predictions for coagulation. The largest opacities in the compact filaments of Taurus (highest column densities in Fig. 11d) would then be due to ice-coated aggregates with a β index 0.1 lower than in the model of Köhler et al. (2015). Directions with large extinctions (A_V > 3.3 mag) in the Taurus clouds are indeed associated with the detection of water ice (Whittet et al. 1988).

Köhler et al. (2015) use aggregates that are about three times larger than the size of the constituting grains (at the peak of the mass-size distribution; see Fig. 1 of their paper). Similarly, Ysard et al. (2013) observed a two-fold increase in opacity at 250 μm towards a molecular filament in the Taurus cloud and modelled it as the result of grain growth by an average factor of five. When aggregates are formed, Köhler et al. (2015) predict a larger R_V value of 4.9 ± 0.2 independent of ice coating. The R_V factors measured by Schlafly et al. (2016) start to increase with increasing reddening for E(B − V) > 2 mag, but they hardly reach values beyond 4 (only 1% of the sample), possibly because the survey does not point specifically to dense molecular regions, but scans larger scales than the compact CO filaments. Foster et al. (2013) have used deep red-optical data from Megacam on the MMT² combined with near-infrared data from UKIDSS to probe the extinction law in two 1°–size sub-structures of the Perseus cloud, known as B5 in the eastern end (at l =160. °6 and b = −16. °8) and L14484-L1451-L1455 in the western end (at l =158. °3, b = −21. °5). They find that the extinction law changes from the diffuse value of R_V ~ 3 at A_V = 2 mag to a dense-cloud value of R_V ~ 5 at A_V = 10 mag in bothsub-structures. These estimates compare with the theoretical value for aggregates and suggest changes in R_V from 3 to 5 with increasing gas density in the CO-bright phase.

The combination of grain growth in size and of lower extinction in the optical band (increased R_V) should lead to a marked decline in specific reddening from the HI to CO media (see Eq. (5)). In the anti-centre clouds, we find a gradual increase rather than a decrease of the mean E(B − V)∕N_H ratio, from (1.84 ± 0.10) × 10⁻²² mag cm² in the HI phase, to (2.09 ± 0.15) × 10⁻²² mag cm² in the DNM, and (2.53 ± 0.12) × 10⁻²² mag cm² in the CO-bright phase. Furthermore, Fig. 9 shows that the changes in E(B − V)∕N_H do not anti-correlate with the N_H map towards CO clouds, but they do towards Perseus. Inside the latter, the specific reddening is 30% lower than the average in the western end and 10% to 20% larger in the eastern end, whereas Foster et al. (2013) find the same increase in R_V in both regions. In the other anti-centre CO clouds, the relative stability of the observed E(B − V)∕N_H ratio may hint at a significant increase in gas-to-dust mass ratio and/or in Q^ext to compensate for the factor 4.8 increase due to grain growth and R_V increase. In their study of the Magellanic Clouds, Roman-Duval et al. (2014) found that depletion measurements suggest that gas-phase metals may accrete onto dust grains in molecular clouds, thus changing the dust-to-gas mass ratio, R_DG. Additionally, R_DG variations may result from the turbulent clustering of dust grains in molecular clouds. Hopkins & Lee (2016) simulated the behaviour of aerodynamic grains in highly supersonic, magnetohydrodynamic turbulence typical of molecular clouds. According to their simulation, the dust can form highly filamentary structures with a characteristic width much thinner than that of gas filaments, and sometimes the dust filaments form in different locations from the gas filaments. These differences result in localised variations of large amplitude in the dust-to-gas mass ratio. Nevertheless, such variations cannot explain the gradual decline in E(B − V)∕τ₃₅₃ as the gas becomes denser (see Eq. (6) and Fig. 11c), so intrinsic changes in the absorption and scattering properties of the grain aggregates in the optical bands must play a role in maintaining high specific reddening values in the dense CO clouds.

Alternatively, the lack of a significant decline in specific reddening with increasing H₂ column density may be due to limitations in the determination of the reddening in obscured regions where few stars are available to sample the extinction. Extinction estimates are derived from magnitude-limited stellar catalogues and therefore suffer from a bias to lower values because the most heavily extinguished stars are not detected (Sale 2015). This is particularly true for lines of sight with higher column densities and it could partially explain the lack of a significant decline in the E(B − V)∕N_H ratio. Moreover, the compact molecular filaments associated with ice mantle formation are not fully resolved in γ rays, nor in reddening (see Fig 6 where ¹³CO emission coincides with the lower resolution pixels). So the actual variations in E(B − V)∕N_H may be steeper than presently observed.

Figure 11 highlights the complex variations of dust properties and dust spectra with environment. Figures 10 and 11 as well as previous studies (Planck Collaboration XI 2014; Planck and Fermi Collaborations Int. XXVIII 2015) indicate that grain evolution severely limits the linear gas-tracing capability of dust emission to the HI phase unless we can understand and model the grain emissivity changes. Here, the extinction properties of the large grains have been shown to be less variable, though still not uniform across the ISM phases. In this context, it is essential to take advantage of the total gas tracing capability of the γ rays produced by cosmic rays to quantify τ₃₅₃∕N_H and E(B − V)∕N_H gradients in other nearby clouds, with both larger and smaller masses than those present in the anti-centre sample, in order to confirm the trends and relative amplitudes of the environmental changes in dust properties. The comprehensive photometry data on stars and quasars from future deep synoptic surveys, such as LSST³, will also clarify the robustness of the reddening measurements towards the densest molecular cores.

Fig. 10 Comparison of the XCO factors derived from dust and γ-ray measurements in mask0 as a function of the average dust optical depth at 353 GHz in the different clouds. The colours mark the use of different gas tracers: E(B − V) reddening in black, τ353 optical depth in blue, γ rays coupled to τ353 in magenta, and γ rays coupled to reddening in red.

Fig. 11 Evolution of the /NH γ specific reddening (top panel), NH γ opacity (second panel), relative variation in the ratio (third panel), and in total gas column density NH γ derived from interstellar γ rays (bottom panel), in 0. °375-wide bins, as a function of the colour temperature, Tdust, and spectral index, β, of the thermal dust emission. The tilde quantities have been convolved with the LAT response for an interstellar spectrum. The variations in the ratio are given with respect to the mean value of 1.49 × 104 mag found from quasar colours through diffuse atomic clouds (Planck Collaboration XI 2014).

Acknowledgements

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. This work performed in part under DOE Contract DE-AC02-76SF00515. The authors acknowledge the support of the Agence Nationale de la Recherche (ANR) under award number STILISM ANR-12-BS05-0016.

Appendix A Choice of the dust extinction map for the analysis

This sectionpresents a comparison of several extinction maps and a study of their correlation with the dust optical depth at 353 GHz from Planck Collaboration XI (2014) and with the hydrogen column density derived from interstellar γ rays in our previous analysis (Remy et al. 2017). The indications given by this comparison have guided the choice of the extinction map to be used with the γ rays to study theE(B − V)∕N_H ratios in thedifferent gas phases.

The main purpose of this preliminary comparison was to quantify the differences between dust extinctions derived with several methods and survey data, with a particular consideration of how the distance-extinction degeneracy affects 2D maps compared to 3D versions, and on how the stellar statistics in the optical and NIR surveys affect the results.

We have selected the three following extinction maps:

• a visual extinction (A_V) map constructed from the 2MASS Stellar data, with the median colour excess of the 49 stars nearest to each direction (Rowles & Froebrich 2009);

• an extinction map in the J band, A_J, derived from the 2MASS data with the NICER M2a processing at 12′-resolution (Juvela & Montillaud 2016);

• a 3D E(B − V) map based on Pan-STARRS 1 and 2MASS stellar data integrated up to its maximal distance (Green et al. 2015).

According to the extinction law of Cardelli et al. (1989), we converted the map of Juvela & Montillaud (2016) from A_J to A_V with a multiplicative factor of 3.55 and applied the extinction factor R_V = 3.1 to convert the map of Green et al. (2015) from reddening to visual extinction. The resulting three A_V maps are presented in Fig. A.1 together with contours delimiting a selection of regions where N_H and A_V are not linearly correlated. The pixel selections in N_H and A_V that were used to generate these contours are illustrated in Fig. A.3 as coloured areas. Figure A.2 shows the 2D histograms of the correlations between each of the A_V maps. Figure A.3 shows the 2D histograms of the correlations between the A_V maps and the dust optical depth at 353 GHz. It also displays the correlations with the gas column density derived from γ rays of interstellar origin in the 0.4–5 GeV energy band from Remy et al. (2017). The datasets were convolved to a common angular resolution and sampled in 0. °25 × 0. °25 bins in order to construct the histograms.

Fig. A.1 From left to right: Visual extinction maps derived from Rowles & Froebrich (2009); Juvela & Montillaud (2016); Green et al. (2015) with the extinction law of Cardelli et al. (1989). The blue and red contours correspond to pixel selections at 0 < AV < 1.6 mag and at AV > 4 mag, respectively, in regions where NH > 4.5 × 1021 cm−2. These selections in AV and NH are shown as the blue and red areas in Fig. A.3. The (l,b) directions that do not match both conditions have been masked to generate the contours that encompass the selected pixels.

Fig. A.2 2D histograms of the correlations between the visual extinction maps derived from Rowles & Froebrich (2009), Juvela & Montillaud (2016), and Green et al. (2015) with the extinction law of Cardelli et al. (1989) (RV = 3.1, AV∕AJ = 3.55). All the maps are compared at a common angular resolution in 0. °25 × 0. °25 bins. The yellow lines correspond to a one-to-one correlation. In each panel the magenta line corresponds to the best-fit of a linear regression.

Fig. A.3 Top: 2D histograms of the correlations between the AV data and the dust optical depth at 353 GHz (Planck Collaboration XI 2014) or the total gas column density, NH γ, measured with 0.4–100 GeV interstellar γ rays (Remy et al. 2017) (bottom). All the maps are compared at a common angular resolution in 0. °25 × 0. °25 bins. The orange line in the top plots corresponds to a ratio AV∕τ353 = 4.62 × 104 mag (Planck Collaboration XI 2014). The green line in the bottom plots corresponds to a ratio AV ∕NH = 5.34 × 10−22 mag cm2 (Bohlin et al. 1978). In each panel the magenta line corresponds to the best fit of a linear regression. The blue and red areas correspond to the pixel selections used to generate the contours in Fig. A.1.

Figure A.2 shows a good correlation between the visual extinctions obtained by Rowles & Froebrich (2009) and Juvela & Montillaud (2016). The slope is close to one, but there is an offset of 0.5 mag in the latter. As the two maps are based on the same dataset, we suspect a method-dependent bias. We note the same offset in the visual extinctions of Juvela & Montillaud (2016) when compared to the map of Green et al. (2015).

The map of Green et al. (2015) presents less noise at low A_V, but most of the pixels with A_V ≳ 1 mag exhibit extinctions above the values of the two others. Those differences may be explained by the use of optical datafrom Pan-STARRS to complement the NIR data from 2MASS, especially in diffuse regions. The precision is improved by larger stellar statistics. Conversely, the NIR data from 2MASS are essential to scan deeper into dense molecular clouds. In the nearby clouds of Orion, Taurus, and Perseus, Green et al. (2015) obtain larger extinctions when combining 2MASS and Pan-STARRS photometry than with the sole use of Pan-STARRS (Schlafly et al. 2014). We also find a tighter correlation between the gas column densities derived in γ rays and the extinctions based on 2MASS and Pan-STARRS data than with Pan-STARRS-based extinctions (not illustrated here).

Figure A.3 clearly shows that the dust optical depth strongly deviates from the visual extinction above 2 or 3 mag, independent of the method employed to derive the visual extinctions. At τ₃₅₃ > 10⁻⁴ the optical depth non-linearly increases above the extinction expectation. At A_V < 2 mag, the distribution of Green et al. (2015) compares well with the A_V∕τ₃₅₃ ratio of 4.62 × 10⁴ mag found by Planck Collaboration XI (2014) using colour-excess measurements towards 53 399 quasars and converted from E(B − V) to A_V with R_V = 3.1.

We note a significant saturation of the extinction at N_H > 4 × 10²¹ cm⁻² and A_V ≈ 1.5 mag in the maps of Rowles & Froebrich (2009) and Juvela & Montillaud (2016). This corresponds to regions in the Galactic disc (enclosed byblue contours in Fig. A.1). As noted by Juvela & Montillaud (2016), near the Galactic plane the zero point of the extinction scale is not well defined because of the unknown mixture of stars and dust along the line of sight. This ambiguity could be resolved only with knowledge of the 3D distribution of the stars and the ISM. Indeed, the reddening map constructed in three dimensions by Green et al. (2015) shows less deviation of the A_V ∕N_H ratio toward the Galactic plane.

At low gas column densities, N_H < 2 × 10²¹ cm⁻², the data from Juvela & Montillaud (2016) and Green et al. (2015) show a better correlation with the gas column densities derived from the γ-ray emission. The slope is consistent with the standard ratio A_V ∕N_H = 5.34 × 10⁻²² mag cm² of Bohlin et al. (1978). At larger column densities, N_H > 4 × 10²¹ cm⁻², the dispersion in A_V increases. Large extinctions deviating from the N_H trend are found in the densest cores of molecular clouds (within red contours in Fig. A.1). They are possibly due to the lack of suitable stars in heavily obscured regions.

As the reddening map of Green et al. (2015) presents a better correlation with both the dust optical depth and gas column densities, we selected this map to study the relation between N_H and E(B − V). The extinction maps from Rowles & Froebrich (2009), Juvela & Montillaud (2016), and more generally all the 2D extinction maps derived from near-infrared surveys such as 2MASS are more suitable to study relative variations in extinction within a cloud or within a region of smaller angular size than the one studied here. Regions near the Galactic plane and very diffuse regions with extinctions lower than about 0.5 magnitude are more prone to extinction biases. Three-dimensional methods adjusting simultaneously the distance and extinction of clouds along the lines of sight appear to be more robust than 2D methods at low Galactic latitudes (Green et al. 2015). Additionally, the combination of optical and NIR surveys allows us to scan a wider range in A_V. On the one hand, optical surveys increase the stellar statistics in translucent regions, which help to reduce the noise level below 0.2 mag in A_V. On the other hand, NIR surveys can probe deeper into the clouds and retain enough precision to A_V ~ 5 mag.

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

References

All Tables

All Figures

	Fig. 1 Maps of the dust optical depth at 353 GHz (top left panel), of the stellar reddening by dust (top right panel), and of their ratio (bottom left panel). The map on the lower right shows in green the pixels with an angular resolution of 0. °11 in E(B − V), which have been retained for mask2, and in black the pixels with a lower resolution (see Sect. 3.3).
In the text

	Fig. 2 CO intensity maps of the Cetus (Cet), Main Taurus (TauM), South Taurus (TauS), North Taurus (TauN), Perseus (Per), and California (Cal) clouds, and of the Galactic disc background (Gal). The solid curve and dashed line show the low-latitude borders of two analysis regions, mask 0 and mask 1, respectively.
In the text

	Fig. 3 Top: γ-ray counts of gaseous origin recorded in the 0.4–100 GeV energy band in a 0. °125 pixel grid. γ-ray emissions other than due to cosmic-ray interactions in the gas have been subtracted. The map has been smoothed with a Gaussian kernel of 0. °14 dispersion for display. Bottom: stellar reddening by dust displayed at the Fermi-LAT angular resolution for comparison.
In the text

	Fig. 4 Hydrogen column density maps obtained for the DNM with the “γ + E(B − V)” (top panel) and “γ + τ353” (bottom panel) analyses of the full region (mask0).
In the text

	Fig. 5 Number distributions of hydrogen column densities in the DNM (black) and COsat (grey) components, obtained from the “γ + τ353” (dashed) and “γ + E(B − V)” (solid) analyses of the full region (mask0). The numbers of pixels for each component have been normalized by the maximum value found in the “γ + E(B − V)” analysis (solid curves).
In the text

Fig. 6 Hydrogen column density map of additional gas in the 12CO-saturated parts of the Taurus and Perseus clouds, derived from the “γ + τ353” (top panel) and “γ + E(B − V)” (bottom panel) analyses of the full region (mask0). The pixels within the grey contours have a lower angular resolution in E(B − V) (mask2). The magenta contours enclose regions with 13CO intensities above 2 K km s−1 (Ridge et al. 2006; Narayanan et al. 2008).

In the text

Fig. 7 Left panel: map of the ratios of hydrogen column densities derived from interstellar γ rays in the “γ + E(B − V)” and “γ + τ353” analyses, thewhite contours marking the positions of γ-ray point sources. Middle and right panels: residual γ-ray maps (photon counts minus best-fit model, in sigma units) in the 0.4–100 GeV band for the “γ + E(B − V)” (middle) and “γ + τ353” (right) analyses, for an HI spin temperature of 400 K. The black contours outline the CO clouds at the 1 K km s−1 level.

In the text

Fig. 8 Upper panel: 2D histogram of the correlation between the total gas column density, NH γ, measured with the 0.4–100 GeV interstellar γ rays, and thedust reddening, , convolved with the LAT point-spread-function for an interstellar spectrum. The maps are sampled at a 0. °375 resolution. The solid line marks the best-fit slope found across the whole field (mask0); it corresponds to a /NH γ ratio of (2.02 ± 0.48) ×10−22 mag cm2 averaged over all gas phases. Lower panel: evolution of the /NH γ ratio in equally spaced bins in NH γ. The black dashed, white dashed, and black solid curves give the means obtained with mask0, mask1, and mask2, respectively.The shaded area gives the standard deviation in each bin using mask2 and the error bars give the standard error on the means.

In the text

Fig. 9 Spatial variations of the amount of reddening per gas nucleon, /NH γ, with the total gas measured by NH γ at 0. °375 resolution (top panel) and its variations relative to the average value of (2.02 ± 0.48) × 10−22 mag cm2 (bottom panel). Quantities with a tilde are convolved with the LAT response for an interstellar spectrum. The black contours outline the shape of the CO clouds at the level of 1 K km s−1 in WCO intensity.

In the text

	Fig. 10 Comparison of the XCO factors derived from dust and γ-ray measurements in mask0 as a function of the average dust optical depth at 353 GHz in the different clouds. The colours mark the use of different gas tracers: E(B − V) reddening in black, τ353 optical depth in blue, γ rays coupled to τ353 in magenta, and γ rays coupled to reddening in red.
In the text

Fig. 11 Evolution of the /NH γ specific reddening (top panel), NH γ opacity (second panel), relative variation in the ratio (third panel), and in total gas column density NH γ derived from interstellar γ rays (bottom panel), in 0. °375-wide bins, as a function of the colour temperature, Tdust, and spectral index, β, of the thermal dust emission. The tilde quantities have been convolved with the LAT response for an interstellar spectrum. The variations in the ratio are given with respect to the mean value of 1.49 × 104 mag found from quasar colours through diffuse atomic clouds (Planck Collaboration XI 2014).

In the text

Fig. A.1 From left to right: Visual extinction maps derived from Rowles & Froebrich (2009); Juvela & Montillaud (2016); Green et al. (2015) with the extinction law of Cardelli et al. (1989). The blue and red contours correspond to pixel selections at 0 < AV < 1.6 mag and at AV > 4 mag, respectively, in regions where NH > 4.5 × 1021 cm−2. These selections in AV and NH are shown as the blue and red areas in Fig. A.3. The (l,b) directions that do not match both conditions have been masked to generate the contours that encompass the selected pixels.

In the text

Fig. A.2 2D histograms of the correlations between the visual extinction maps derived from Rowles & Froebrich (2009), Juvela & Montillaud (2016), and Green et al. (2015) with the extinction law of Cardelli et al. (1989) (RV = 3.1, AV∕AJ = 3.55). All the maps are compared at a common angular resolution in 0. °25 × 0. °25 bins. The yellow lines correspond to a one-to-one correlation. In each panel the magenta line corresponds to the best-fit of a linear regression.

In the text

Fig. A.3 Top: 2D histograms of the correlations between the AV data and the dust optical depth at 353 GHz (Planck Collaboration XI 2014) or the total gas column density, NH γ, measured with 0.4–100 GeV interstellar γ rays (Remy et al. 2017) (bottom). All the maps are compared at a common angular resolution in 0. °25 × 0. °25 bins. The orange line in the top plots corresponds to a ratio AV∕τ353 = 4.62 × 104 mag (Planck Collaboration XI 2014). The green line in the bottom plots corresponds to a ratio AV ∕NH = 5.34 × 10−22 mag cm2 (Bohlin et al. 1978). In each panel the magenta line corresponds to the best fit of a linear regression. The blue and red areas correspond to the pixel selections used to generate the contours in Fig. A.1.

In the text