© M. Elyajouri and R. Lallement 2019

1 Introduction

One of the most challenging problems in terms of the composition and evolution of interstellar medium (ISM) is the ubiquitous presence of the numerous unidentified diffuse interstellar bands (DIBs) detected in ours and other galaxies. Since their first discoveries in the 1920s by Heger (1922), these absorption features seen in spectra of Milky-Way reddened stars or in galaxy spectra have puzzled astronomers and spectroscopists (see McCall & Griffin 2013 for a recent historical review). In spite of an apparent general correlation between the DIB strengths and the interstellar extinction or reddening, all efforts have failed up to confirm the hypothesis that dust particles are responsible for the absorption of DIBs. One can cite in particular the most stringent limitations on polarization reported by Cox et al. (2007) who analyzed six strong DIBs (at 5780, 5797, 6196, 6284, 6379 and 6613 Å) and revealed no linear polarization, a polarization expected for large grain-related carriers. The lack of correlation between DIB strengths and the far-UV (FUV) extinction disfavors the small grain-related carriers as well (Desert et al. 1995; Xiang et al. 2017). Diffuse interstellar band carriers are therefore believed to be large molecules in gas phase. The detection of substructures in the profile of some DIBs (e.g., Elyajouri et al. 2018, and references therein) supports the molecular nature of some DIB carriers.

The most promising carrier candidates are thought to be carbon chains, polycyclic aromatic hydrocarbons (PAHs), and fullerenes. Indeed, at least four of five bands predicted to be associated with the buckminsterfullerene C${}_{60}^{+}$ have now been detected, making this fullerene cation the first identified DIB carrier (Foing & Ehrenfreund 1994; Cordiner et al. 2017; Campbell et al. 2016; Lallement et al. 2018).

It remains that the global correlation between most DIBs and optical extinction, A_V, suggests that dust grains that cause this extinction are associated with DIB carriers in some way. However, the underlying mechanisms and role of DIB carriers in the multiple physico-chemical reaction chains that occur throughout the ISM life cycle are still unknown and many open questions remain to be answered: is the DIB carrier density related to the depletion onto the grains? Could grains be the formation site of DIB carriers?

Fortunately, the new massive stellar spectroscopic surveys can now support statistical studies allowing comparisons between absorption by gaseous species (e.g., CaII, NaI, KI, C₂, CH, CH⁺), DIB carriers, and dust, helping to clarify their links. Such statistical studies of relationships between DIBs, gas, and dust grains, as well as ratios between pairs of diffuse bands, contain precious information on the molecular content and the evolution of macro-molecules and grains in response to the physical properties, the radiation field, the abundances in the ISM, and the grain size distribution. A hierarchy of DIBs ordered by increasing sensitivity to the ionizing field is currently being elaborated (e.g., Ensor et al. 2017). The link with the ionizing field is particularly interesting as it is related to the depth within the clouds, and thus to the grain properties.

An important observational clue is that some of the strongest optical DIBs have shown to level off with respect to dust reddening for lines-of-sight that probe denser molecular cores rather than the diffuse edges of interstellar clouds. Using the Lick Observatoryscanning spectrophotometer, Wampler (1966) was the first to notice such an effect for the DIB at 4430 Å in the spectra of background stars to dark clouds. This was also seen by Snow & Cohen (1974) whilst studying DIB strengths at 4430, 5780, and 5797 Å in Sco-Oph and Per. These latter authors concluded that formation or survival of the carrier is inhibited in cloud interiors. Since then, the same effect has been documented in a number of studies (Strom et al. 1975; Meyer & Ulrich 1984; Wallerstein & Cardelli 1987; Adamson et al. 1991; Herbig 1995). This behavior was referred to as the skin effect in the Herbig (1995) review and is supported by the higher degree of correlation of the DIB strength with HI compared to molecular hydrogen (Friedman et al. 2011). Moreover, the discovery of unusually weak DIBs in the spectra of the star HD 62542 by Snow et al. (2002) also confirmed such an interpretation, if, as is likely, the diffuse outer layers of the dark cloud in front of the star have been stripped away. However, a particular class of DIBs, the so-called C2-DIBs (Thorburn et al. 2003), seems to be abundant in the dense molecular phase (Ádámkovics et al. 2005). For recent reviews including the skin effect see Snow & McCall (2006) and Snow (2014).

More recently, Lan et al. (2015) reported this behavior for many DIBs towards lines of sight intercepting high-latitude molecular clouds using the Sloan Digital Sky Survey (SDSS) star, galaxy, and quasar spectra. These authors concluded that the break seen at high E(B–V) values can be quantitatively characterized by the decrease of HI In dense clouds, and higher molecular fraction.

Finally, important additional findings on the skin effect have been brought by Fan et al. (2017). The authors of this latter publication performed a large study of eight different, strong optical DIBs and produced a statistically significant study of the DIB-to-extinction ratio as a function of various atomic and molecular tracers of the ISM phase. In particular, they found that the relationship between the DIB-to-extinction ratio and the molecular hydrogen fraction f(H₂) is characterized by a lambda-shaped curve with a peak at f(H₂) ≃ 0.3, unequivocally showing that the DIB carriers reside predominantly in this range of molecular fraction and tend to decrease in both fully molecular cores and atomic gas. As a result, this study reveals how the physical state of the gas controls the DIB-to-extinction ratio. Here, the term skin effect incorporates all of these factors, since there is indeed a maximum volume density of DIB carriers for a particular molecular fraction, and a decrease above and below this specific fraction value. In what follows, we only essentially address the depletion in the dense molecular phase, and for this reason we refer to a depletion effect instead of a skin effect.

At least one of the near-infrared (NIR) DIBs, which allow us to probe more heavily obscured regions, and thus denser (higher extinction) clouds, also displays a depletion in the dense phase. Using VLT/KMOS moderate-resolution (R = 4000) NIR (H-band) spectra of 85 stars located behind the Barnard 68 dark globule – from the edge to the center – a tomography of the 15 273 Å NIR DIB in an individual dark interstellar cloud was performed for the first time (Elyajouri et al. 2017a), demonstrating the decrease of the 15 273 DIB carrier in the denser parts of such dark clouds.

Previous efforts to observe the skin effect either cover small regions around individual clouds, or have insufficient sample sizes to statistically constrain its amplitude. Taking advantage of the increasing number of stellar surveys with highly multiplex instruments, we developed an entirely automated fitting method to extract as much information as possible from stellar spectra. Such methods have been applied to NIR spectra of hot stars, in particular the telluric calibrators of the SDSS/ Apache Point Observatory Galactic Evolution Experiment survey in the H band (Elyajouri et al. 2016, hereafter EMRL16).

In this work, we expand the studies to the latest data release DR14 using spectra of red giant stars and covering a significant fraction of the sky. Our approach to the problem of the DIB–extinction link involves the use of over 100 000 individual lines of sight, focusing exclusively on the strongest DIB in the H-band, at 15 273 Å. We begin the following section with a description of the public catalogs used in this study, including a brief discussion of the changes and novelties made in the various surveys. Section 3 describes the methods used for DIB extraction and the resulting catalog. This catalog will be discussed in detail in a separate publication, with special focus on its applications. In Sect. 4 we revisit the relationship between DIB strength and extinction, showing a break at increasing extinction values. Special attention is paid to slight differences in the stellar models of the various pipelines, since correction for this variation is crucial to achieve a good DIB extraction, and we discuss the extinction estimates. Using our results in combination with existing extinction estimates and an all-sky dust map, we investigate in Sect. 5 both the relationship between the DIB–extinction ratio and the dust optical depth measured by Planck. Our main conclusions are summarized in Sect. 7.

2 Description of the public catalogs used in this study

2.1 SDSS/APOGEE spectra

The high-resolution (R ~ 22,000) IR Apache Point Observatory Galactic Evolution Experiment (APOGEE, Majewski et al. 2017) is one of the programs in the Sloan Digital Sky Survey III (SDSS-III, Eisenstein et al. 2011). Released APOGEE spectra have a high signal-to-noise ratio (S/N) (generally above 100) and are remarkably well corrected for telluric emission and absorption with residuals generally below 1% as can be seen in figures from Casey et al. (2016). The latest APOGEE Data Release 14 (DR14; Abolfathi et al. 2018) includes results from APOGEE-1 (September 2011–July 2014) in addition to two years of SDSS-IVAPOGEE-2 data (July 2014–July 2016), and DR14 APOGEE products contain data from approximately 263 000 stars, all identified by a 2MASS catalog name. Among the APOGEE stars, there are 231 000 main scientific targets, mostly red giants intended for various scientific programs and forming part of the bulge, the bar, the disk, and the halo. There are also 27 000 young and hot stars referred to as telluric standards stars (TSSs) which were used to correct telluric absorptions (Zasowski et al. 2013). These hot stars are distributed in all the fields observed by APOGEE (35 stars/field) and were chosen among the bluest in each field. Of these 27 000 stars, 17 000 from APOGEE Data Release 12 (DR12; Holtzman et al. 2015) were the subject of a detailed DIB study by Elyajouri et al. (2016, 2017b). They are particularly interesting because one can easily identify any interstellar signature¹. It is more difficult to extract the DIB from the late-type stars since the interstellar features are blended with deep stellar lines. This explains why most of our reliable DIB-fitting results of earlier studies were obtained for hotter stars. However, if the stellar models reproduce the data well, the extraction can be performed without ambiguities due to stellar blends. This is the case for the synthetic spectra of the APOGEE DR14 which includes an innovative data-driven technique called Cannon.

2.2 Cannon model spectra

Highly detailed model spectra were fitted to the APOGEE data and made publicly available for a large number of targets. These fitted models greatly facilitate the extraction of the diffuse interstellar bands present in the APOGEE spectral intervals. The DR12 synthetic spectra used the APOGEE Stellar Parameter and Chemical Abundances Pipeline (ASPCAP; García Pérez et al. 2016) for the derivation of the stellar atmospheric parameters and abundances.In DR14 a new and important change is the introduction of an innovative data-driven machine-learning technique known as Cannon (Ness et al. 2015). This algorithm parameterized the spectral fluxes as a function of a set of independently determined stellar parameters and abundances. For DR14, Cannon-2 code by Casey et al. (2016) was trained on ASPCAP products for a subset of high-S/N giant stars, and the model was applied to all objects within the range of parameters covered by the training set with some modifications listed in Abolfathi et al. (2018).

Here we analyze all available red giant stars (160 000) from the public Cannon data products. To do so we use the CannonStar file which contains Cannon labels, normalized observed spectra, and associated uncertainties, as well as synthetic spectra for each single star. Cannon results such as these can be found at path² /cannon/LOCATION_ID/cannonStar*.fits. The difference between Cannon and ASPCAP has significant consequences on our DIB extraction, as is detailed in the following section.

2.3 Starhorse extinction and distance database

In order to study the link between the 15 273 Å DIB strength and the color excess (Sect. 4), we made use of the stellar extinctions in the visible range A_V calculated with the Bayesian Starhorse method developed by Santiago et al. (2016) and Queiroz et al. (2018). The authors applied the Starhorse code to APOGEE DR14 data and the associated stellar parameters derived from the Cannon method, and therefore there is a consistency in using the results of the StarHorse technique (here version 1) applied to APOGEE DR14 and the Cannon DR14 model spectra.

Starhorse also contains spectrophotometric distances or parallax distances from the Tycho-Gaia Astrometric Solution (TGAS, Gaia Collaboration 2016) when available, as described in Anders et al. (2018). We also made use of these distances (see Sect. 5). All catalogs are available to the community via the LineA web page³.

2.4 Planck Galactic dust data

To trace Galactic dust properties in regions with DIB measurements and allow better comparisons with the interstellar dust along the sightlines, we made use of the Planck Legacy archival⁴ data. The Planck Collaboration used Planck and WMAP data to model the spectral energy distributions (SEDs) of the dust thermal emission and derive an all-sky map of the dust optical depth at 353 GHz, τ₃₅₃, and of the dust temperature, T_D. Here we used the Planck 2015 results described in Planck Collaboration X (2016).

3 DIB extraction

3.1 Automated fitting method

We used the APOGEE DR14 normalized spectra and their associated best fit Cannon models to extract the 15 273 Å DIB through simple Gaussian profile fitting. The general principles are the same as in the case of the automated code applied to TSSs of APOGEE DR12 and are described in EMRL16. Here we kept the same fitting technique and the same method of classification as in EMRL16 and Capitanio et al. (2017), that is, we fitted the data to the product of an adjusted synthetic spectrum $S_{\lambda }^{\alpha }$ (here Cannon instead of ASPCAP), a Gaussian DIB model DIB[σ, λ, D], and a smooth continuum (1 + [A] × λ) as shown in Eq. (1). $$M_{\lambda }={[S_{\lambda }]}^{\alpha }\times \text{DIB}[\sigma ,\lambda ,D]\times (1+[A]\times \lambda ),$$(1)

where each DIB is modeled as a Gaussian function with three free parameters associated with its Gaussian RMS width (σ), central wavelength (λ_c), and depth (D) as follows: $$\text{DIB}[\sigma ,\lambda ,D]=1+D\times e^{\frac{-{(\lambda -{\lambda }_c)}^2}{2{\sigma }^2}}.$$(2)

Prior to the final global fit, the DIB extraction code performs preliminary series of automated tests that allow one to measure the residuals between the data and the adjusted stellar model in three spectral regions: the local region around the DIB A = [λ_c − 10–λ_c + 10] Å, a region free of DIBs and telluric residuals B = [15 911–15966] Å, and a narrower region free of DIBs and telluric residuals and closer to the DIB C = [15 200–15 239] Å. The associated standard deviations R_A, R_B, and R_C are used for selection criteria and estimates of the DIB uncertainties. During the fit itself the adjusted Gaussian curve representing the DIB is rejected or selected depending on its depth D, its width σ, and spectrum quality criteria derived from the preliminary tests.

Stellar model adjustment S${}_{\lambda }^{\alpha }$

The Cannon stellar models are optimized for red giant stars, and they perfectly reproduce the observed data. This was not the case for the TSSs that are bluer stars for which it was necessary to adjust the stellar models. Here we simplified our constraints on the adjustment of the stellar synthetic spectra. We still allow for a global increase or decrease of all stellar lines, however the prior value for the adjustment coefficient α is always unity. The fitted coefficient is always very close to one and the continuum is also always very close to a straight line of value 1.

Fig. 1 Illustration of DIB extraction towards three red giants stars. In each figure: APOGEE spectra are shown with solid black curves. The initial stellar model provided by APOGEE Cannon is shown in solid blue, the DIB absorption alone is represented in orange (with an offset), and the solid magenta curves represent the final stellar+DIB modeled spectrum.

DIB depth upper limit

TSSs are nearby stars whereas the APOGEE main targets can be very distant and DIBs can be much stronger. For this reason no upper limit was applied to the DIB depth.

DIB depth lower limit

The stellar model is so close to the data in all spectral regions that the standard deviation R_A of the data-model residuals in the region close to the DIB (region A) is a good enough estimate of the measurement error on the DIB depth. As a consequence the new acceptance criterion is no longer based on a comparison between the DIB depth and the standard deviations of residuals in regions A, B, and C, but is now solely: Depth D ≥ R_A.

DIB radial velocity range

TSSs being nearby stars, the absorbing interstellar matter for those targets is characterized by a narrow radial-velocity range. This is no longer the case here for the distant giant stars. As a consequence we widened the allowed interval on radial velocities with respect to the one used for TSSs. Specifically, we now restrict the spectral range for the DIB center to [15 260–15 286] Å.

DIB width upper limit

We no longer used HI 21cm data to limit the total width of the DIB, instead we used limits based on the FWHM histograms from our study based on TSSs and Zasowski et al. (2015). The latter correspond to a broad range of star distances, broad enough to be valid here.

DIB width lower limit

We kept the same lower limit of 0.7 Å on the DIB width, which is supposed to be close to the intrinsic width as determined by Zasowski et al. (2015) and EMRL16 and with some allowance for noise-related uncertainty. Figure 1 shows some examples of our fitting method.

3.2 DIB upper limit estimate in case of no convergence or no DIB detection

If the adjustment described previously failed by absence of convergence (we refer to it as the non-convergence case), or there was convergence but the fitted absorption was found to be smaller than the noise level (we refer to it here as the upper-limit case, although as we see below for some spectra at the final stage of analysis this term is not appropriate) then we re-initiated a second fitting procedure with the goal of extracting equivalent widths (EWs) or upper limits in a more efficient way with more appropriate initial conditions. The total number of targets falling in the two categories and re-analyzed is 45 393. First, we fitted the combination of a continuum and a depth-variable stellar model to the data (i.e., removing a DIB contribution from the global model), this time masking the DIB area. In case of convergence, this ensures that the fitted continuum and the synthetic spectrum accurately reproduce the data outside the DIB.

• Case 1: if the fit fails again – this is the case when the model is not adapted to the data or if there is a problem of continuum – a new standard deviation of residuals is calculated based on the difference between the data and the initial model provided by the APOGEE pipeline.

• Case 2: if the fit converges then the new residuals (difference between the data and the no-DIB adjusted model) are used for an updated estimate of the standard deviation.

The second step consists in fitting a simple Gaussian to the data-model residuals in the region of the DIB. To do so, we use the same initial guesses and constraints used for the automated fitting method as described in EMRL16.

• If the Gaussian fit does not converge, either because it has not found a solution or because there is a real absence of detectable DIB, then a linear fit is performed, that is, on the continuum over the wavelength range chosen for the DIB EW [15 260–15 286] Å. The selected value for the EW is then zero and an upper limit on the EW is estimated from the standard deviation of the residuals around the linear fit Std_dev and the average value of the width of the 15 273 Å DIB taken as 1.7 Å:

$$\text{EW}=0\pm 1.7\times \sqrt{2\pi }\times {\text{Std}}_{\text{dev}}.$$(3)

Some representative examples are shown in Fig. A.1.

• If the fit finds a solution, then the coefficients of the fitted Gaussian (σ_Gauss; D_Gauss) is used either as DIB parameters in case 1 or as initial values for a new automated fit as detailed in EMRL16 in case 2. Equivalent widths in both cases will be derived as EW= $\sqrt{2\pi }D\sigma$ where the widths are σ = {σ_Gauss;σ_fit} and the depths are D={D_Gauss;D_fit}, respectively. The errors are extracted as follows:

$${\text{Error}}_{\text{DIB}}=2\sqrt{\pi }\times \sigma \times {\text{Error}}_{\text{depth}}+2\sqrt{\pi }\times \sigma \times {\text{Std}}_{\text{dev}}.$$(4)

These cases are illustrated in Fig. A.2. Finally, for a better understanding of our algorithm, Fig. A.3 shows a graphical representation that maps out our step-by-step process to determine the upper limits. It is interesting to note that this additional treatment is quite productive: for approximately two out of three cases, the change of initial guesses for the DIB in the fitting procedure allows the Levenberg–Marquardt algorithm to converge to a good solution and gives a significantly smaller χ² than for the first fit. This is explained by the existence of large Doppler shifts between the star and the absorbing ISM. As a matter of fact, in our first adjustment the initial condition for the DIB center was a null radial velocity of the absorbing medium in the stellar frame. In case of large differences between the star radial velocity and the absorbing ISM radial velocity, the initial position of the DIB center was accordingly distant from the actual DIB center, and in some cases the adjustment could subsequently stop at a secondary χ² minimum. This no longer happens in the second adjustment with an initial guess for the DIB center now very close to the best fit at the primary χ² minimum. Similarly, peculiar or poorly modeled stellar features could favor a secondary χ² minimum during the first adjustment, when they were combined with an initial position of the DIB far from the actual one. Using instead an initial guess much closer to the solution leads to the optimal solution.

Fig. 2 Equivalent width of the 15 273 Å DIB as a function of stellar extinction AV. Black points are the individual sightline measurements. The red circles are the EW medians in each ΔA V = 1.6 mag, and the yellow line is the fitted EW/AV relationship.Superposed are the background stars of Barnard 68 (the larger blue circles).

3.3 The catalog

The extraction procedure allowed us to extract EW measurements of the 15 273 Å DIB for 124 064 different sight-lines. The catalog of EWs and associated uncertainties, radial velocities, and widths are available at the Strasbourg Data Center (CDS) after merging with the updated catalog of DIBs in the TSS spectra of DR14, when available.

We note that because the DIB is a weak absorption, despite the high quality of the spectra and stellar models, measurement uncertainties on the EW become on the order of the EW itself for EW ≲ 0.02–0.05 Å. This corresponds to A_V ≲ 0.5–0.8 mag according to the average relationship (see below). This will result in some limitations in our analyses and correlative studies.

4 Diffuseinterstellar band and dust extinction

Taking advantage of improvements to models and spectra, we have revisited the link between 15 273 Å DIB strength and the extinction.

To do so, we cross-matched our catalog of 15 273 Å DIBs (Sect. 3) and the Starhorse Cannon database using TOPCAT (Taylor 2005). We removed stars that are flagged with NUMMODELS_BAD, EXTINCTION_BAD, and BRIGHT2MASS_WARN in Starhorse catalog. We finally obtained 124 064 stars with DIB detection and photometric extinction and distance.

Figure 2 shows the EW of the DIB as a function of the extinction along the sightline to the target. As expected and firstly discovered by Zasowski et al. (2015) for this DIB, there is a strong correlation between these two quantities that extends over more than two orders of magnitude. We computed the median value of the EW in extinction bins of 1.6 magnitudes,and fitted this median value to a power-law relationship with the extinction bin centers. Using the new catalog we find the following relationship: $${\text{EW}}_{\text{DIB}}=0.056\times A_{\text{V}}^{0.99\pm 0.03}\pm 0.001(\AA ),$$(5)

for the range of A_V=[0.55–13.35] mag. As explained above, our lower threshold is associated with EW error intervals. On the other hand, the upper limit for the linear fit is constrained by the limited number of targets with high extinctions and chosen to avoid poor statistics. This relationship is significantly different from the previous power-law formula derived by Zasowski et al. (2015) found over three orders of magnitude in both Av and EW: $${\text{EW}}_{\text{DIB}}=0.102\times A_{\text{V}}^{1.01}(\AA ).$$(6)

We investigated the reasons for this difference by comparing the EW results for the same individual stars and by visual inspection of the corresponding spectra and models. Figure 3 shows one of these individual comparisons. The figure shows that there are changes in the adjusted model spectra in many regions, and that by chance one of the most significant differences happens at the location of the 15 273 Å DIB. The Cannon stellar lines are deeper than the corresponding ASPCAP lines, and as a result the fitted DIB has a smaller EW. We further investigated the impact of these differences by cross-matching the catalog of Zasowski et al. (2015) and our catalog for the same target stars. Figure 4 displays the EWs derived by Zasowski et al. (2015) and in this work based on the ASPCAP and Cannon stellar spectra, respectively. Both are displayed as a function of the same extinction estimate, here the Cannon value. The slope of the average linear relationship between EW and the extinction is found to be about 15% smaller in our case compared to the Zasowski et al. (2015) derivations. This may be due to the deeper stellar lines in the Cannon spectral model and the subsequent weaker EWs. This interpretation is strongly reinforced by the observed dependence of the discrepancy between the two determinations on the DIB strength; see Fig. 5. The figure displays the average of the ratio between our determination and that of Zasowski et al. (2015) in bins of increasing EW. It is clearly seen that, the stronger the DIB EW, the smaller the difference between the two determinations. We subsequently conclude that a non-negligible fraction of the difference between our two relationships can be attributed to the use of different synthetic stellar models.

A second reason for the discrepant slopes is the use of different extinction estimates. Figure 6 displays the RJCE and Cannon values of theextinction A_V for the target stars listed in Table 1 of Zasowski et al. (2015). The average linear relationship between the two estimates is shown in the figure. According to this, RJCE values are on the order of 25% above the Cannon estimates for the set of stars. The two effects, that due to the use of the Cannon extinction instead of the RJCE estimate on the one hand, and that due to the slightly different Cannon spectrum on the other hand, result in a global decrease of the EW-to-extinction ratio and explain why our DIB-to-extinction ratio is about 40% below that found by Zasowski et al. (2015).

Figure 2 also shows that the linear relationship is no longer valid at high A_V values, and that an increasing number of DIBs have EWs that are significantly lower than predicted by this law. This is the so-called skin effect first detected by Snow & Cohen (1974) and then often found in the spectra of highly reddened stars. This effect is illustrated here in Fig. 2 using averages of DIB EWs in bins of 1.6 mag of A_V extinction. This choice of A_V interval is a compromise between the need for a minimum number of targets per bin and the maximum discretization over the whole range.

This depletion is better seen when zooming at the break region for high A_V values (right part of Fig. 2). The median value of the DIB-to-extinction ratio (equivalent to the ratio between the column of DIB carrier and the column of grains) is decreasing while A_V is increasing. This is clearly due to an increasing number of sightlines with relatively weak DIBs. On the other hand, a large fraction of the sightlines have strong DIBs that remain quasi-proportional to the extinction.

It is interesting to compare the skin effect seen here in a large variety of sightlines with the specific case of a dark cloud, that is, a quasi-perfect mono-cloud situation in which all the extinction is generated in a unique dense cloud. To do so, we added to the figure the results obtained in Elyajouri et al. (2017a) for the Barnard 68 dark cloud. The corresponding DIB EW data points correspond to some of the largest departures from the values predicted by the linear relation. However, they are still comparable to tens of other similar cases, suggesting that the marked disappearance of DIBs in Barnard 68 is far from an exception.

Fig. 3 Spectrum of 2MASS J19420307+2317393 from ASPCAP (solid black) and Cannon (dashed line) data products. The associated synthetic stellar models are shown as dashed curves. We clearly see the difference between ASPCAP (blue) and Cannon (red) models around the 15 273 Å DIB location (black arrow). An offset was applied to ASPCAP data and model to facilitate the comparison.

Fig. 4 Comparison between the DIB EWs derived in this work and those from the Zasowski et al. (2015) catalog. Both values are plotted as a function of the Cannon extinction value AV. Linear fits to both datasets (dashed blue and pink lines) are superimposed.

Fig. 5 Average ratio between the EW derived in this work using the Cannon models and the EW derived by Zasowski et al. (2015) based on ASPCAP models, here in bins of EWs from our work. The ratiois increasing and reaches an asymptotic value on the order of 0.95 for large EWs above 300 mÅ.

Fig. 6 Cannon and RJCE extinction estimates for target stars from Table 1 from Zasowski et al. (2015). The RJCE estimate used by these authors and listed in their table is ≃ 25% higher than the Cannon estimate used in the present work.

5 The DIB carrier depletion illustrated with Planck

A fraction of the APOGEE sightlines is out of the Plane and crosses the nearby clouds of the Gould belt. For those targets that are located behind the main clouds, the observed dust emission is entirely or almost entirely generated by the nearby clouds. For this reason it is interesting to investigate the link between the DIB-to-extinction ratio and the dust emission. One of the goals is to detect variations among the cloud complexes or conversely some general laws for the skin effect.

Figure 7 is a map of the Planck optical depth τ₃₅₃ in the region of the nearby Orion clouds (data from Planck Collaboration X 2016). Superimposed on the map are the APOGEE targets, with markers color-coded according totheir EW/A_V ratio. To ensure that the stars are located behind the Orion clouds, targets closer than 600 pc were eliminated. In addition, we discarded very low reddening targets and/or very low DIB EWs in order to avoid the corresponding uncertain ratios. The corresponding thresholds are A_V ≥ 0.2 mag and EW ≥ 0.03 Å. The EW DIB decrease in the dense phase is conspicuously illustrated in this figure: towards all cloud cores (τ₃₅₃ ≥ 10⁻⁵) the DIB-to-extinction ratio decreases dramatically. Reciprocally, in the more diffuse regions the ratio reaches values that are above the average ratio of ~0.06 Å mag⁻¹. This is the first time such a clear visualization is possible, and is due the combination of high-quality APOGEE data and the high quality and spatial resolution of the Planck maps. Similar figures, this time for the Taurus and Cepheus clouds are presented in Figs. 8 and 9. The same restrictions were applied on A_V and EW, and the lower limits for the target distances were chosen as 1 kpc for Taurus and 1.3 kpc for Cepheus. The same striking association between cloud cores and low DIB–extinction ratios is immediately visible.

A visual inspection of the three figures does not reveal any marked difference between the three cloud associations. In order to better quantify the DIB level-off, Fig. 10 shows log(EW/A_V) as a function of log(τ₃₅₃) for the three regions and the same selected targets located beyond the clouds. A first attempt to quantify the level-off was made by means of a linear fit in the log–log scale for the merged datasets. The resulting best-fit relationship between EW/A_V and τ₃₅₃ was found to be $$\text{EW}/A_{\text{V}}=7.6\times {10}^{-5}\times {\tau }_{353}^{-0.66\pm 0.01}(\AA {\text{mag}}^{-1}).$$(7)

We also performed the same adjustments for the three regions separately, this time taking into account uncertainties on the DIB/EW ratio (the fractional error on the ratio is estimated as the square root of the quadratic sum of fractional errors on the two terms). The adjustments are displayed in Fig. 10 and the results were: $$\begin{array}{l}\text{Taurus}:\text{EW}/A_{\text{V}}=5.3\times {10}^{-5}\times {\tau }_{353}^{-0.70\pm 0.01}(\AA {\text{mag}}^{-1}),\\ \text{Orion}:\text{EW}/A_{\text{V}}=1.0\times {10}^{-4}\times {\tau }_{353}^{-0.65\pm 0.02}(\AA {\text{mag}}^{-1}),\\ \text{Cepheus}:\text{EW}/A_{\text{V}}=1.6\times {10}^{-4}\times {\tau }_{353}^{-0.59\pm 0.02}(\AA {\text{mag}}^{-1}).\end{array}$$

These quantitative results show that the average relationships present some differences from one region to the other, with a slope varying from −0.70 to −0.59 and withoutoverlap of the allowed intervals.

In a search for the origin of such differences, we made use of 3D maps of the interstellar dust. Such maps are based on extinction measurements for stars at distributed distances and assign locations to the main cloud complexes (Lallement et al. 2019). Figure 11 displays the dust distribution in vertical planes crossing the clouds. For each complex, two vertical cuts are shown, for two longitudes chosen to correspond to large concentrations of target star sightlines. In each vertical plane, the range of Galactic latitudes of the target stars selected in Fig. 8 are indicated by shaded sectors. A black line delimits the maximum distance from the Sun to the densest regions in the corresponding sector. The dust differential opacity threshold has been chosen here to be 1.5 mag pc⁻¹. The figures show that most sightlines cross several dense structures, that is, a sightline does not probe a unique gas phase; instead, series of dense and more diffuse regions are distributed along the sightlines. This probably explains the large dispersion of data points in Fig. 10, since DIB-to-extinction ratios for the various phases will be smeared outin different ways depending on the sightlines. However, interestingly, the distance from the Sun to the dense cloud limit increases from Taurus to Orion then to Cepheus (although the increase from Taurus to Orion is much stronger than from Orion to Cepheus), that is, this distance increases while simultaneously the slope of the relation of DIB-to-extinction ratio with the dust thickness gets shallower. An increase in the limit distances and the number of structures along the sightlines implies that the mixing of diffuse and dense regions must be stronger in Orion than in Taurus and even higher in Cepheus compared to Orion. For those reasons, we interpret the change of slope among the three clouds shown in Fig. 10 as an indicator of the strength of phase mixing. In the case of Taurus, with its smaller number of structures and smaller resulting phase mixing, the level-off of the DIB is more distinctly linked to the dust column and the slope of the decrease is stronger. The opposite is true for Cepheus. If our interpretation is correct, this implies that the actual DIB-to-extinction ratio decrease from the diffuse to the dense phase is even stronger than suggested by the slopes in Fig. 10.

Finally, we attempted to use the Planck data and all APOGEE sightlines to study the DIB level-off in all types of configurations. To do so we removed all conditions on the target distances and kept only the lower limits on A_V and EW imposed to avoid overly large uncertainties on the DIB-to-extinction ratio. The total number of sightlines becomes 75 840. Figure 12 shows the DIB-to-extinction ratio as a function of τ₃₅₃. The Orion, Taurus, and Cepheus targets discussed above are marked by black signs, while all other targets are color-coded according to the Galactic latitude. Interestingly, one can distinguish three different regimes: at high latitude, thus for low dust optical depth and mostly nearby clouds, the DIB-A_V ratio apparently stabilizes, that is, it no longer increases with decreasing τ₃₅₃. It reaches a limit on the order of ≃ 0.23 Å mag⁻¹. Such a value is probably representative of very diffuse regions located close to the boundary of the Local Cavity (see, e.g., maps from Capitanio et al. 2017). The second regime corresponds to the trend observed for the off-plane Taurus, Orion, and Cepheus regions, as already discussed; however this time it also applies to a large number of targets within the Plane. The third regime corresponds to large optical depths and sightlines within the mid-plane (as indicated by the low values of Galactic latitude). Here the ratio is no longer depending on τ₃₅₃ and reaches a limit on the order of 0.05 Å mag⁻¹, a value very close to the average found from the global linear fit to the data. This value must correspond to a strong mixing of diffuse and dense phases all along the sightlines.

Superimposed in Fig. 12 are the Barnard 68 corresponding data points for the five regions delineated in Elyajouri et al. (2017a) that span values of A_V between 2 and 20 mag. To convert A_V into τ₃₅₃ we used the average relationship A_V ≃ 4.5 × 10⁴ τ₃₅₃, based on the classical A_V -to-E(B–V) ratio of 3.1 and the average E(B–V) versus τ₃₅₃ relationship of Remy et al. (2018). It can be seen from the Fig. 12 that these data points fall in the middle of the ensemble of data points for Taurus, Orion, and Cepheus, that is, again, the dark globule does not correspond to an exceptional situation in terms of DIB carrier depletion. This is interesting since Barnard 68 is immersed in a large opaque region with only a few ionizing stars, in strong opposition to Orion.

Fig. 7 Two-dimensional map of the Galactic interstellar dust in the Orion region (Planck 353 GHz optical depth, color scale at bottom right). Superimposed are the SDSS/APOGEE target stars contained in the field of view. The black and white coding of the markers follows the ratio between the 15 273 Å DIB strength and the extinction of the target light in the visible, with white (resp. black) indicative of a high (resp. low) ratio. Towards the cores of the clouds (in red) there is a spectacular decrease of the DIB/extinction ratio. At the cloud periphery, and also outside the clouds that are visible in the Planck image, the ratio is enhanced, reflecting the so-called “skin effect”.

Fig. 8 As in Fig. 7 but for the Taurus region.

Fig. 9 As in Fig. 7 but for the Cepheus region.

Fig. 10 DIB EW/AV ratio towards stars behind Taurus (top), Orion (middle), and Cepheus (bottom) cloud complexes as a function of the dust optical thickness measured by Planck τ(353 GHz). Linear fits in the Log–Log scale are shown in red.

6 Summary

We have analyzed red giants NIR spectra from the last SDSS/APOGEE survey data release (DR14) with the goal of extracting the 15 273 Å DIB. To do so, we used the code developed by EMRL16 with some adaptations motivated by the change of target types, namely from early- or intermediate-type stars to red giants. The use of the stellar models computed by means of the data-driven Cannon algorithm greatly improves the analysis. In the region of the 15 273 Å DIB the adjustment has proven to be good enough to allow the extraction of the DIB EW, width, and Doppler velocity through Gaussian fitting of the interstellar absorption. Careful and severe examinations of the DIB parameters, the continuum shape, and the quality of the adjustment were done. This resulted in a conservative selection of reliable DIB parameters for a total of 124 064 lines of sight. We have compared our DIB measurements with their associated visual extinctions extracted from the Starhorse database. We revisited the DIB−A_V relationship first reported by Zasowski et al. (2015) and found the DIB strength average to be linearly correlated with dust extinction over at least two orders of magnitude, however with a smaller EW-to-extinction ratio than the one derived by these authors for the red giants and later by EMRL16 for the hotter TSS targets: EW${}_{\text{DIB}}=0.056\times A_{\text{V}}^{0.99\pm 0.03}(\AA )$. We showed how this variation is explained by the differences between the ASPCAP and Cannon models and by the use of Cannon A_V extinctions instead of A_V extinctions deduced from the RJCE technique.

It has also been found from the global DIB–extinction comparative study that the 15 273 Å DIB tends to weaken with respect to extinction in the spectra of stars situated behind dense regions. With the goal of going one step further in the analysis of the DIB carrier depletion in dense cores, we investigated the link between the EW/A_V ratio and the distribution of dust as traced by Planck. In particularly, the Taurus, Orion and Cepheus cloud complexes were studied using targets distant enough to be located behind the clouds. The depletion effect is spectacularly illustrated for those cloud complexes, thanks to the high quality of the APOGEE data and of the Planck maps. We have fitted an average relationship between the DIB/A_V ratio and the Planck optical depth τ₃₅₃ through a log–log linear formula: $EW/A_{\text{V}}=7.6\times {10}^{-5}\times {\tau }_{353}^{-0.66\pm 0.01}(\AA {\text{mag}}^{-1})$.

We performed a separate study of the three cloud complexes in a search for potential differences in the DIB depletion effect. For the Taurus clouds (resp. Orion, Cepheus), and taking into account uncertainties on EWs and extinction, the coefficients of the log–log linear relationship are: (5.3 × 10⁻⁵;−0.70 ± 0.01), (1.0 × 10⁻⁴;−0.65 ± 0.02) and, (1.6 × 10⁻⁴;−0.59 ± 0.02), respectively. Using 3D maps of the dust clouds obtained by inversion of stellar individual extinctions (Lallement et al. 2019), we have shown that the decrease of the depletion gradient from Taurus to Cepheus is very likely due to a corresponding increasing number of structures encountered along the paths to the target stars behind the three complexes, and the subsequent increase of phase mixing. Our results provide first empirical average relationships between the 15 273 Å DIB carrier depletion in molecular clouds and the dust optical depth.

Finally, we used the entire dataset in an attempt to identify different depletion regimes for the APOGEE sightlines. We found that, not surprisingly, for large dust optical depths exceeding ~ 10⁻⁴, that is, for low latitude and extended sightlines, the DIB/A_V ratio tends to the mean value of about 0.05–0.06 × A_V which is characteristic of a strong mixing of diffuse and dense phases distributed along the extended paths. For intermediate τ₃₅₃ values between ~ 10⁻⁵ and 10⁻⁴, we clearly see the same depletion effect that is affecting the nearby off-plane cloud complexes, now extended to midplane targets. Finally, for low values of τ₃₅₃ below 10⁻⁵, the DIB/A_V ratio seems to be constant at ~ 0.2 × A_V, a value significantly above the strong mixing average. This high value corresponds to the nearby diffuse clouds at the periphery of the local cavity.

Fig. 11 Three-dimensional distribution of dust in vertical planes containing the Sun and crossing the Taurus, Orion, and Cepheus clouds. The north Galactic pole is toward the top in each figure. The cloud distribution is based on Lallement et al. (2019). For each cloud group, two vertical planes are selected in such a way that the number of target sightlines falling very close to at least one of the two planes is maximized. The Galactic longitudes corresponding to each plane are indicated in the bottom-right corner of each figure. The dashed black line encloses the densest structures (see text).

Fig. 12 DIB to extinction ratio EW/Av as a functionof τ353 for all APOGEE sightlines with EW ≥ 0.03 Å and Av ≥ 0.2 mag. These limitations ensure a reliable determination of the EW/Av ratio. Orion, Taurus, and Cepheus targets are marked by black signs, while all other targets are color-coded according to the Galactic latitude. Here, one can distinguish three different regimes: high-latitude, low dust sightlines that correspond to nearby diffuse clouds, and a high EW/Av of ≃ 0.23 Å mag−1; partially mixed phases in nearby off-plane cloud complexes following the Taurus-Orion-Cepheus trend (red curve); and distant, high-dust-optical-depth sightlines within the midplane for which the EW/Av ratio reachesa constant value of ≃ 0.056 Å mag−1 (blue dashed curve).

7 Discussion

This break seen at high A_V values is well known as the skin effect (Herbig 1995) and has been observed several times for optical DIBs. In the case of the 15 273 Å DIB, it was also observed recently in the B68 nearby dark cloud (Elyajouri et al. 2017a) with a strength comparable to multiple sightlines from the present study. Very recently, it has been more specifically characterized and quantified for optical DIBs by Fan et al. (2017) who demonstrated the existence of a maximum carrier abundance with respect to dust grains for a molecular H₂ fraction on the order of 0.3.

The origin of the DIB carrier depletion in dense molecular cores is still unclear. In the context of PAH-like carriers, UV radiation would certainly directly influence the charge state of PAHs (e.g., Cami et al. 1997; Ruiterkamp et al. 2005) and the resulting disappearance of some DIBs could simply be due to the evolution of the charge state from the outer to inner regions of the clouds. On the other hand, photo-processing and hydrogenation or de-hydrogenation may come into play (e.g., Vuong & Foing 2000). More recently, Bertaux & Lallement (2017) suggested that a general carrier disappearance could additionally be due to coagulation of DIB carriers into aggregates and accretion onto dust grains, as predicted by Jones (2016) in the context of new models of dust. It is well known that grains grow in compact cores, generating the so-called coreshine effect (Lefèvre et al. 2014), and that grain coalescence to form aggregates may also involve macro-molecules associated with DIBs. Future detailed studies of the DIB carrier depletion in regions with differing radiation field, pressure, temperature, and grain properties will hopefully help to disentangle the different contributions. In this respect, empirical relationships involving dust emission such as those established here may complement observed links with the physical state of the gas, such as those found by Fan et al. (2017), and bring additional clues.

Acknowledgements

M.E. acknowledges funding from the “Region Ile-de-France” through the DIM-ACAV project. R.L. acknowledges support from the CNRS PCMI national program. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support andresources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatory of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group,Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Appendix A Derivation of DIB EWs in the case of very small absorptions (illustrative examples)

Fig. A.1 Illustration of EW upper limit derivation, in the case of nonconvergence of the Gaussian fit (see Sect. 3.2). A Gaussian DIB profile with the maximum EW compatible with the standard deviation in the data–stellar model residuals is shown at an imposed wavelength of 15 272.42 Å in solid black (with an offset).

Fig. A.2 Illustration of DIB extraction in the case of convergence after derivation of an optimal DIB prior (the initial guess based on the preliminary computation of data–stellar model residuals and its Gaussian fit). The fitted DIB is shown in dashed orange with an offset.

Fig. A.3 Flowchart compiling our decision criteria to determine the DIB EWs in the upper limit and nonconvergence cases. Part of the spectra analyzed here are those classified as upper limit in our automated method described in EMRL16. The width σfit and the depthfit are the Gaussian parameters derived from the same automated method described in EMRL16. Stddev is the standard deviation derived from the linear fit to residuals over the wavelength region defined to measure the EW of the DIB.

References

All Figures

	Fig. 1 Illustration of DIB extraction towards three red giants stars. In each figure: APOGEE spectra are shown with solid black curves. The initial stellar model provided by APOGEE Cannon is shown in solid blue, the DIB absorption alone is represented in orange (with an offset), and the solid magenta curves represent the final stellar+DIB modeled spectrum.
In the text

	Fig. 2 Equivalent width of the 15 273 Å DIB as a function of stellar extinction AV. Black points are the individual sightline measurements. The red circles are the EW medians in each ΔA V = 1.6 mag, and the yellow line is the fitted EW/AV relationship.Superposed are the background stars of Barnard 68 (the larger blue circles).
In the text

	Fig. 3 Spectrum of 2MASS J19420307+2317393 from ASPCAP (solid black) and Cannon (dashed line) data products. The associated synthetic stellar models are shown as dashed curves. We clearly see the difference between ASPCAP (blue) and Cannon (red) models around the 15 273 Å DIB location (black arrow). An offset was applied to ASPCAP data and model to facilitate the comparison.
In the text

	Fig. 4 Comparison between the DIB EWs derived in this work and those from the Zasowski et al. (2015) catalog. Both values are plotted as a function of the Cannon extinction value AV. Linear fits to both datasets (dashed blue and pink lines) are superimposed.
In the text

	Fig. 5 Average ratio between the EW derived in this work using the Cannon models and the EW derived by Zasowski et al. (2015) based on ASPCAP models, here in bins of EWs from our work. The ratiois increasing and reaches an asymptotic value on the order of 0.95 for large EWs above 300 mÅ.
In the text

	Fig. 6 Cannon and RJCE extinction estimates for target stars from Table 1 from Zasowski et al. (2015). The RJCE estimate used by these authors and listed in their table is ≃ 25% higher than the Cannon estimate used in the present work.
In the text

Fig. 7 Two-dimensional map of the Galactic interstellar dust in the Orion region (Planck 353 GHz optical depth, color scale at bottom right). Superimposed are the SDSS/APOGEE target stars contained in the field of view. The black and white coding of the markers follows the ratio between the 15 273 Å DIB strength and the extinction of the target light in the visible, with white (resp. black) indicative of a high (resp. low) ratio. Towards the cores of the clouds (in red) there is a spectacular decrease of the DIB/extinction ratio. At the cloud periphery, and also outside the clouds that are visible in the Planck image, the ratio is enhanced, reflecting the so-called “skin effect”.

In the text

	Fig. 8 As in Fig. 7 but for the Taurus region.
In the text

	Fig. 9 As in Fig. 7 but for the Cepheus region.
In the text

	Fig. 10 DIB EW/AV ratio towards stars behind Taurus (top), Orion (middle), and Cepheus (bottom) cloud complexes as a function of the dust optical thickness measured by Planck τ(353 GHz). Linear fits in the Log–Log scale are shown in red.
In the text

Fig. 11 Three-dimensional distribution of dust in vertical planes containing the Sun and crossing the Taurus, Orion, and Cepheus clouds. The north Galactic pole is toward the top in each figure. The cloud distribution is based on Lallement et al. (2019). For each cloud group, two vertical planes are selected in such a way that the number of target sightlines falling very close to at least one of the two planes is maximized. The Galactic longitudes corresponding to each plane are indicated in the bottom-right corner of each figure. The dashed black line encloses the densest structures (see text).

In the text

Fig. 12 DIB to extinction ratio EW/Av as a functionof τ353 for all APOGEE sightlines with EW ≥ 0.03 Å and Av ≥ 0.2 mag. These limitations ensure a reliable determination of the EW/Av ratio. Orion, Taurus, and Cepheus targets are marked by black signs, while all other targets are color-coded according to the Galactic latitude. Here, one can distinguish three different regimes: high-latitude, low dust sightlines that correspond to nearby diffuse clouds, and a high EW/Av of ≃ 0.23 Å mag−1; partially mixed phases in nearby off-plane cloud complexes following the Taurus-Orion-Cepheus trend (red curve); and distant, high-dust-optical-depth sightlines within the midplane for which the EW/Av ratio reachesa constant value of ≃ 0.056 Å mag−1 (blue dashed curve).

In the text

	Fig. A.1 Illustration of EW upper limit derivation, in the case of nonconvergence of the Gaussian fit (see Sect. 3.2). A Gaussian DIB profile with the maximum EW compatible with the standard deviation in the data–stellar model residuals is shown at an imposed wavelength of 15 272.42 Å in solid black (with an offset).
In the text

	Fig. A.2 Illustration of DIB extraction in the case of convergence after derivation of an optimal DIB prior (the initial guess based on the preliminary computation of data–stellar model residuals and its Gaussian fit). The fitted DIB is shown in dashed orange with an offset.
In the text

Fig. A.3 Flowchart compiling our decision criteria to determine the DIB EWs in the upper limit and nonconvergence cases. Part of the spectra analyzed here are those classified as upper limit in our automated method described in EMRL16. The width σfit and the depthfit are the Gaussian parameters derived from the same automated method described in EMRL16. Stddev is the standard deviation derived from the linear fit to residuals over the wavelength region defined to measure the EW of the DIB.

In the text