Cosmic-rays, gas, and dust in nearby anti-centre clouds : III – Dust extinction, emission, and grain properties

,


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, τ 353 . 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, τ 353 /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 3.4 +0.3 −0.7 (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 2 gas that is not linearly traced by CO intensity, W CO , in the directions where the 12 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 2 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 , 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 anticentre 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 τ 353 /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 τ 353 /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 τ 353 /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 are given in Sect. 4 and compared to the results of the former "γ + τ 353 " 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 τ 353 /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.

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 in the 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).

Dust emission and extinction
The dust optical depth is defined as the product of the dust emission cross section per unit mass, κ ν in cm 2 g −1 , and the dust mass column density, M dust . The emission cross section can be expressed as κ ν = κ 0 ν ν 0 β with κ 0 being the value at the reference frequency ν 0 and β the spectral index (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 µ H N H with µ H = 2.27 × 10 −27 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: and then we can define the dust emission opacity or hereafter in short, opacity (Planck Collaboration XI 2014), as: 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 σ ext λ (Spitzer 1998): Assuming the grains can be described by spheres of effective radius a and mean mass density ρ gr , we have σ ext λ = Q ext λ πa 2 with Q ext λ the efficiency factor for extinction. The extinction optical depth can therefore be expressed as a function of the hydrogen column density as: Substituting this relation in the expression for the visual extinction, A V , and using the definition of the total-to-selective extinction factor R V = A V E(B−V) yields an expression of the specific reddening: A relation between the optical depth and the reddening can be obtained by combining Eqs. (1) and (5): Provided we can reliably trace the total gas column densities in all gas forms, variations in τ 353 /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, τ 353 , 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 τ 353 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 τ 353 , reddening E(B − V), and of their E(B − V)/τ 353 ratio are displayed for the region of interest in Fig. 1 3-1 mag) are shown to have a ratio close to the value of 1.49 × 10 4 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).

γ -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 1 parameter file 54-LRYusifovXCO4z6R30-Ts150-mag2 has been used to generate an energy-dependent template of the Galactic IC emission across the analysis region . 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).

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 12 CO J = 1→0 observations at 115 GHz from the momentmasked 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 −1 in the local standard of rest. The EBHIS frequency sampling is coarser, with a velocity resolution of 1.44 km s −1 , 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).
Ionised gas 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).

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 τ 353 for E(B − V): 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. N DNM H (l, b) and N CO s at H (l, b) 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 χ 2 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 isotropic intensity 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 −2 s −1 sr −1 MeV −1 , can be modelled as : with E(B − V) DNM and E(B − V) CO sat 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 5 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 factors and 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 obtained through γ-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 τ 353 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 τ 353 , small differences between the two maps in Fig. 3

Mapping DNM and CO sat 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 2 transition and the CO sat H 2 gas missed by the saturation of the 12 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 of the 7 K km s −1 contour in W CO intensity.
The DNM and CO sat templates of additional gas 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 τ 353 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.

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)/τ 353 ratio result from uncertainties in the zero level in τ 353 , fluctuations from the cosmic infrared background in the τ 353 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 τ 353 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 2 conversion factor (X CO = N H 2 /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.

Gas column density
The coupled analysis 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 freefree 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, then in the total gas. Figure 4 shows that the hydrogen column densities derived in the DNM from the "γ + E(B − V)" and "γ + τ 353 " 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 20 cm −2 in the "γ + τ 353 " map, and larger column densities, probably from the Galactic disc background, above 1.5 × 10 21 cm −2 near the Galactic plane (b > −8 • ) in the "γ + E(B − V)" map. But the total mass in the "γ + τ 353 " map is only 9% lower than in the 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 21 cm −2 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 "γ + τ 353 " one. Conversely, the densest filaments and clumps emitting in the rarer 13 CO isotope are enhanced in the "γ + τ 353 " map compared to the "γ + E(B − V)" one (beyond 6 × 10 21 cm −2 ). This difference results from the combination of the sharp rise in τ 353 /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 "γ + τ 353 " 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 "γ + τ 353 " 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 model excess noted towards the brightest CO parts of Main Taurus, Perseus, and California clouds in the "γ + τ 353 " 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 "γ + τ 353 " was caused by the marked rise in τ 353 /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.

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 τ 353 /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 "γ + τ 353 " 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 is very 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 −22 mag cm 2 of Bohlin et al. (1978) that applies to lines of sight with E(B − V) < 0.5 mag Fig. 6. Hydrogen column density map of additional gas in the 12 COsaturated 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 13 CO intensities above 2 K km s −1 (Ridge et al. 2006;Narayanan et al. 2008). in the solar neighbourhood and in the Galactic disc. They also agree with the value of (1.71 ± 0.07) × 10 −22 mag cm 2 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 −22 mag cm 2 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 −22 mag cm 2 in the denser DNM at the HI -H 2 interface (given by the best-fit parameters of the γ-ray model q HI /q DNM ). FUSE observations probed the same type of medium in diffuse H 2 clouds up to 1 mag in reddening (Rachford et al. 2009). The specific reddening changed between (1.68 ± 0.10) × 10 −22 mag cm 2 and (2.38 ± 0.96) × 10 −22 mag cm 2 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 −22 mag cm 2 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 −22 mag cm 2 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 τ 353 /N H ratios towards the same clouds, at the same angular resolution, and with the same analysis method (Remy et al. 2017). The implications of this result on dust grain properties will be discussed in Sect. 5.

X CO factors
The best-fit coefficients of the γ-ray and dust models yield independent 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 "γ + τ 353 " 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, τ 353 , 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 2 and CO molecules. The models predict strong inward gradients in the relative abundances of CO and H 2 , 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 The solid line marks the best-fit slope found across the whole field (mask0); it corresponds to a E(B − V)/N H γ ratio of (2.02 ± 0.48) × 10 −22 mag cm 2 averaged over all gas phases. Lower panel: evolution of the E(B − V)/N H γ ratio in equally spaced bins in N H γ . 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. τ 353 value partially captures their diffuse (low τ 353 ) or compact (high τ 353 ) character. The systematic decline seen in X CO with all gas tracers clearly confirms the presence of CO-to-H 2 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 τ 353inferred 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 "γ + τ 353 " analysis because the masked pixels have large column densities in the τ 353 -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 20 cm −2 K −1 km −1 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 τ 353 estimates of X CO differed more at larger τ 353 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.

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 τ 353 /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 of 0. • 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)/τ 353 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)/τ 353 ratio by the mean slope of 1.49 × 10 4 mag found in the correlation between τ 353 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.

Dust evolution in the diffuse medium
In the "γ + τ 353 " analysis of the nearby anti-centre clouds (Remy et al. 2017), we showed that the emission opacity τ 353 /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 −27 cm 2 H −1 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 −27 cm 2 H −1 (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 22 cm 2 mag of Bohlin et al. (1978) towards the sight lines intercepting atomic gas. They increase to an average of (2.09 ± 0.15) × 10 −22 mag cm 2 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 22 cm 2 mag in Figs. 11 and 9 correspond to very slightly reddened regions (<0.08 mag) where the difficult reddening measurements should be taken A71, page 11 of 17 A&A 616, A71 (2018) 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. 2015Ysard et al. , 2016Kö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, aromaticrich 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 τ 353 /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 . Such constraints are compatible with the small dispersion in E(B − V)/N H found in the atomic and DNM phases of the anticentre 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 2 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)/τ 353 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)/τ 353 ratio is also expected to decrease as the ambient radiation field hardens . 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.

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 3 cm −3 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 .
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 2 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 both sub-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 −22 mag cm 2 in the HI phase, to (2.09 ± 0.15) × 10 −22 mag cm 2 in the DNM, and (2.53 ± 0.12) × 10 −22 mag cm 2 in the CO-bright phase. Furthermore, Fig. 9 shows that the changes in E(B − V)/N H do not anticorrelate 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-togas mass ratio. Nevertheless, such variations cannot explain the gradual decline in E(B − V)/τ 353 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 2 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 13 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 τ 353 /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 3 , will also clarify the robustness of the reddening measurements towards the densest molecular cores.

Appendix A: Choice of the dust extinction map for the analysis
This section presents 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 the E(B − V)/N H ratios in the different 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 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. 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 data from 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 . 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 τ 353 > 10 −4 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 /τ 353 ratio of 4.62 × 10 4 mag found by Planck Collaboration XI (2014) using colourexcess 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 21 cm −2 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 by blue 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 21 cm −2 , 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 −22 mag cm 2 of Bohlin et al. (1978). At larger column densities, N H > 4 × 10 21 cm −2 , 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.  Green et al. (2015) with the extinction law of Cardelli et al. (1989). The blue and red contours correspond to pixel selections at 0 < A V < 1.6 mag and at A V > 4 mag, respectively, in regions where N H > 4.5 × 10 21 cm −2 . These selections in A V and N H 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.  (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.