© ESO, 2014

1. Introduction

The study of low-mass stars and of planet formation starts by understanding the place where they form and evolve, that is, inside dense molecular clouds. There, the gas and the dust are in constant interaction through collisions that can lead to heat exchange and, in suitable conditions, to the freezing of gas molecules onto dust grains. Inside the molecular cloud, the dust content is known to evolve mainly via grain growth: by accretion of heavy gas particles onto the dust grains that increases the total dust mass (Hirashita 2012), and by the presence of sticky ice mantles (volatile species frozen onto the grains, Walmsley et al. 2004), which favors coagulation (Ossenkopf 1993; Ormel et al. 2009). Ice mantle formation beyond A_v ~ 3 mag (Whittet et al. 2001, 2013) implies a change in the grain properties throughout molecular clouds. In addition, interstellar grains evolve with time during the prestellar phase. They continue to grow or possibly reach a stationary state in the cloud envelope (Ormel et al. 2009), while in the densest region, the dust evolution becomes complex in the presence of an embedded object such as Class 0 or Class I protostars (André et al. 2000). The thermo-mechanical action of the protostellar object will affect the grains by shattering them in the outflows (Anderl et al. 2013) and by evaporating the grain mantles, releasing water (and other species) into the surrounding gas (Fraser et al. 2001), as seen by Herschel (Kristensen et al. 2012; van der Tak et al. 2013). In this context, to infer the molecular cloud stage from the dust properties is a complex problem that begins by understanding the dust grain content.

Grain properties can be investigated in different manners, for example, via a characterization of the extinction curve. In the optical and UV, this curve changes depending on the dominant grain growth mechanism: accretion or coagulation (Hirashita & Voshchinnikov 2014). In the near-infrared (NIR) and mid-infrared (MIR) ranges, this change is a clue to coagulation (Chapman et al. 2009; Ascenso et al. 2013) and to the presence of ice mantles (McClure 2009). Nevertheless, grain growth deduced indirectly from the extinction curve is sensitive to the wavelength normalization (Fitzpatrick & Massa 2007). In this context, a very efficient way to learn about the properties of larger grains (micrometer size) is through the recently discovered effect of MIR dust scattering, or “coreshine effect” (Pagani et al. 2010; Steinacker et al. 2010). The coreshine effect is widespread and detected in at least half of the molecular clouds investigated by Pagani et al. (2010) and Paladini et al. (in prep.), and thus can be used as a tool to explore the properties of the dust responsible for this phenomenon.

Coreshine is observed at those MIR wavelengths, 3–5 μm, where the scattering by large grains is strong enough to be seen in emission. The best examples are seen in the 3.6 and 4.5 μm Spitzer Infrared Array Camera (IRAC) filters (Pagani et al. 2010). When the absence of emission in the 5.8 and 8 μm IRAC filters excludes the presence of polycyclic aromatic hydrocarbons (PAHs), there is no need to consider an active emission process in the modeling, and only scattering and absorption have to be treated, as presented in Steinacker et al. (2010). The same restriction to absorption and scattering modeling is also pertinent for shorter wavelengths (optical and NIR), as demonstrated by Lehtinen & Mattila (1996) for a standard interstellar radiation field (ISRF).

The investigation potential of coreshine, combined with the modeling of other wavelengths, provides an opportunity to better constrain both cloud structures and dust properties. Indeed, thanks to the low opacities at MIR wavelengths and to the anisotropic scattering, coreshine provides access deep inside the clouds and brings information on their 3D structure. The feasibility of the method has been shown (Steinacker et al. 2010; Andersen et al. 2013), but to quantify our capacity of building a 3D model of a real cloud and deduce grain properties from a combination of wavelengths, we need to evaluate the impact of the free parameters and of the boundary conditions on the modeling. Here, we compare our models with observations and deduce general trends on the grain properties for the regions where most of the clouds present coreshine. The set of dust models will constitute a future data base to start modeling molecular clouds in absorption and scattering from visible to MIR. In a forthcoming paper, we will model a specific cloud while also including the far-infrared (FIR) emission to further constrain the dust properties and cloud structure.

In Sect. 2, we present a summary of the observations and the strategy adopted to analyze them. In Sect. 3, we describe the Monte Carlo radiative transfer code used for the simulations, focus on the importance of constraining the radiation field, and present the cloud model and dust content, which are the free parameters to be explored. In Sect. 4, we describe the results obtained from the observational data and compare them with our grid of models. We discuss the coreshine phenomenon and what a multiwavelength approach can add to investigating the grain properties. We present our conclusions in Sect. 5.

2. Observations and analysis

2.1. NIR data

NIR data have been obtained with WIRCAM, the CFHT Wide IR CAMera. Its field of view is 20′ × 20′ with a pixel size of 0.3′′. It is large enough to cover each of the selected targets, including a large area around them. The observations were obtained (and continue to be taken) for a total of eight sources (L183 and a set of sources in the Taurus region) in the standard J, H, and K_s (hereafter K) spectral bands. The observations, data processing, and data themselves will be presented in detail in a forthcoming paper. We used some J- and K-band observations to have a first look at the comparison with models, without performing any exact fits.

2.2. Coreshine data

To investigate the presence of coreshine phenomenon inside a sample of molecular clouds, we used Spitzer data taken by IRAC (Fazio et al. 2004) at 3.6, 4.5, 5.8 and 8 μm (Cryogenic mission) complemented by that from the Warm mission (3.6 and 4.5 μm only). When possible, the 3.6 and 4.5 μm data were compared with Spitzer 8 μm and/or WISE 12 μm maps¹ (Wright et al. 2010) to define the size of the core in absorption, and optical data (in R or B) to delimit the size of the surrounding cloud. These Spitzer observations were collected at different epochs and show a median frame-integration time per pixel from 50 s (from molecular cores to planet-forming disks (c2d) – Evans et al. 2009) to 1800 seconds (Hunting Coreshine Survey, HCS, Spitzer cycles 8 and 9 – Paladini et al., in prep.). The sensitivity for a unit exposure-time in the warm mode (HCS) is similar to that in the cold mode (P94 program: Search for Low-Luminosity YSOs – Lawrence & Keene 2004, c2d)², which allowed us to merge the observations. The aim of the HCS proposal was to obtain an unbiased sample, and our analysis confirms the ~50% detection rate of coreshine as in Pagani et al. (2010,Paladini et al., in prep.).

The coverage is entirely dependent on the program and on the target. The Spitzer field of view size is 5.12′× 5.12′ with a native pixel size of 1.2′′. The full width at half maximum (FWHM) of the point spread function (PSF) varies from 1.66′′ at 3.6 μm to 1.98′′ at 8 μm. For the nearby molecular clouds, the maps are large enough to include both the cores and their environment. The sensitivity for an extended source at 4.5 μm is about 90% of the 3.6 μm sensitivity for both warm and cold campaigns². However, the data must be corrected for column pulldown and zodiacal emission³ (Paladini et al., in prep.). It must be noted that even if Spitzer data are absolutely calibrated using point sources, the extended emission calibration can be as uncertain as 10%. Indeed, while the drift of the electronics zero level is small, the zodiacal light estimate remains slightly uncertain because of its time dependence and the spatial resolution of the model used to derive it. Although this is of no importance for differential measurements, such as the coreshine intensity defined as the net signal above the background, it becomes a problem for evaluating the absolute background intensity behind the clouds. Indeed, the observational value I_back,Sp can vary by a factor of two from one campaign to another because the uncertainty on the zodiacal light estimate is on the same order as the background value itself (and consequently cannot be used as a reference for I_back in the modeling, Sect. 3.2.2).

The sample of the 215 sources presented in Table C.1 and in Fig. 1 is essentially a compilation of previous surveys including the c2d and P94 surveys, which were the main archives for the identification of coreshine sources in Pagani et al. (2010), the new HCS survey (Paladini et al., in prep.), plus a few more targets of interest detected thanks to WISE (Wright et al. 2010) or recently identified in the Spitzer Archive (including the Gum/Vela region, Pagani et al. 2012). The total number of clouds that show coreshine is 108 with some preferential directions on the sky. Indeed, the Taurus, Perseus, Aquila, and Aries complex (hereafter Taurus-Perseus) reveals almost 100% coreshine detection (Fig. 2, Table C.1). Cepheus and Chameleon regions also display a large portion of coreshine cases. Finally, because of its high Galactic latitude and its complexity, L183 will be also one of the regions of interest for this study. Among these 108 positive coreshine detections, we chose to ignore the cores located in the Orion/Monoceros and ρ Oph regions where the coreshine occurrence drops to ~50% (Table C.1). In these regions the coreshine occurrence is likely dominated by the effect of local sources, which make the radiation field difficult to constrain, adding another degree of freedom. Therefore, they are beyond the goal of our global study of coreshine. With this selection criterion, we are left with 72 coreshine sources.

2.3. Analysis

From the observations, we quantified the scattering excess above the surrounding background. To proceed with the measurement, two steps are mandatory: i) removal of the point sources to keep only the extended emission; ii) subtraction of the background contribution (I_back,Sp) under the core to deduce the remaining intensity. We used the SExtractor software (Bertin & Arnouts 1996) and tuned its parameters to remove point sources and keep only the extended emission.

First, SExtractor computes the background from a large scale mesh. This coarse background is resampled to the original resolution, the fluctuations are dampened by applying a median filter to the mesh points. The result is subtracted internally from the original image. Next, the object detection step works on a smoothed image and requires to set up a detection threshold. The objects are detected from an image smoothed with a Gaussian kernel that has a FWHM equal to 1.5 pixels. Each detection is required to contain at least 3 pixels above 1.5σ. The difficulty here is to adjust the background mesh size so that the coreshine emission is removed as part of the background while the extended saturated stars remain so that they can be detected as objects by SExtractor. If the mesh size is too small, they will appear in the final image. On the other hand, if the mesh size is too large, the extended coreshine emission would not be included in the background estimate, but be detected instead as a source. The object identification file was checked, and for most coreshine cases, we chose a mesh size of 32 × 32 pixels, but sometimes we adapted the mesh size to 64 × 64 instead. Then, the object list can be used to either remove the sources from the image for display or to produce a mask of the point sources for measurements. This source subtraction method with the optimized parameters works on almost all Spitzer maps (even in crowded fields) and allows us to mask or to subtract compact objects and to retrieve only extended emission (Fig. 3).

Because coreshine is a phenomenon defined as a signal in excess with respect to the surrounding background emission, the coreshine intensity has to be background subtracted. To perform this subtraction, it is mandatory to interpolate the background value (I_back,Sp) at the core position. Because of the extended nature of the coreshine effect, the process cannot be automatized. After subtracting the point sources, we masked the coreshine region by hand and subtracted either a plane by least-squares fitting or repeated the background subtraction process with SExtractor, interpolating across the coreshine-masked region. This interpolation is usually safe when the masked region remains small, but can become less accurate for large sources. However, the interpolation of the background is normally close to a plane. We compared both subtraction methods for the very extended coreshine emission of L183 and found good agreement between the two methods, with a difference smaller than 10%.

The 3.6 and 4.5 μm cleaned images were convolved with a Gaussian kernel of 10′′ FWHM to reduce the noise by a factor of ~6. Then we measured the peak flux for the 3.6 μm image and built the 4.5 μm / 3.6 μm ratio image, eliminating all pixels with a 4.5 μm flux below 3σ (after smoothing). This is a good compromise to avoid the strong fluctuations due to noise. We computed the histogram of the ratio map and took the maximum and the FWHM to characterize its dispersion (Fig. 3). We selected four regions of interest (Fig. 2) where most of the clouds present coreshine. We applied this method to the 72 cores from these regions that show a positive coreshine signal (Table 1).

3. 3D radiative transfer modeling

Many 3D dust radiative transfer codes are available today (see review by Steinacker et al. 2013). In this section, we first present one of these 3D codes (Continuum Radiative Transfer – CRT)⁴. We discuss the main input parameters: the radiation field, the cloud model, and the dust properties. While the radiation-field parameter is important to constrain, we focus on the dust property variations to see the influence of grain growth (especially coagulation) in reproducing the observational trends region by region.

3.1. Monte Carlo radiative transfer code

The radiative transfer calculations were made with the CRT program (Juvela & Padoan 2003; Juvela 2005). This implements the basic Monte Carlo scheme where a number of photon packages are sent out from each of the radiation sources, the propagation and scattering of the photons is followed, and the intensity of the scattered radiation exiting the medium is registered. The density field can be discretized using spherical or cylinder geometries (see Ysard et al. 2012) or, as in the case of this paper, using full 3D Cartesian grids. CRT allows the use of multiple dust populations and spatial variations of their abundance. In addition to the density field, only the dust optical opacity (e.g., relative to hydrogen), the albedo, and the scattering phase function need to be specified to calculate the scattering. For the scattering phase function, that is, the probability distribution of the scattering angles, CRT allows the use of the Henyey-Greenstein approximation with the asymmetry dust parameter g = ⟨ cos θ ⟩, or one can use any scattering function tabulated as a function of the scattering angle θ. We used the usual Henyey-Greenstein approximation. Its validity is discussed in Sect. 3.3.2.

In our calculations, the main source of radiation is the interstellar radiation field that is described with all-sky DIRBE maps (Hauser et al. 1998) in HEALPix format (Górski et al. 2005), using a separate map to describe the sky brightness at each of the simulated wavelengths (Sect. 3.2.1). The influence of using a single wavelength instead of several averaged wavelengths taking into account the filter response has been tested and shows no statistical differences for the coreshine wavelengths. The simulation runs were made with 100 million photon packets, and we estimated the numerical uncertainty on the modeling result to be 1 kJy sr^-1. The original positions of the emitted photon packages are weighted with the sky intensity, so that more packages (with correspondingly lower weight or smaller photon number) are generated from the Galactic plane and, in particular, from the direction of the Galactic center. CRT uses the forced-first scattering method (Mattila 1970) to ensure adequate sampling of scattered flux in regions of low opacity. To improve the quality of the scattered-light images, the peel-off technique (Yusef-Zadeh et al. 1984) was used, where after a photon package is scattered, CRT always explicitly calculates the fraction of photons that scatter toward the observer and escape the cloud without further interactions (I_sca, Sect. 3.2.2).

The images of scattered light were built using the peeled photons, and they represent the surface brightness visible for an observer far outside the cloud. CRT has an option for calculating peel-off images for several directions during the same run. In the present case, images were only calculated in one direction as determined by the relative locations of the selected clouds and the observer. In practice, the cloud model was viewed along one coordinate axis and the background DIRBE maps were rotated so that the illumination geometry was correct.

3.2. Incident and background interstellar radiation fields

To model dust extinction and emission, the ISRF has to be determined. Two different quantities are needed: the sky brightness in all directions, which determines the illumination of the cloud, and the surface brightness behind the cloud, which determines the net effect of absorption along the line of sight of our observations. The background-field determination requires a precise treatment, which we discuss in more detail below.

3.2.1. Incident ISRF

The all-sky illumination has to take into account the contribution of stellar sources; O and B stars dominate the UV field (Habing 1968) and mostly K stars/red giants for the longer wavelengths. The diffuse part due to ambient stellar light scattered from small grains, UV light reprocessed in PDRs, and PAHs emission must also be considered. Different Galactocentric distances can lead to different intensity estimates (Mathis et al. 1983), but taking into account the fact that the molecular clouds we study are close to us (with a distance range from 100 pc – L183 – to 325 pc – Cepheus) we made the approximation that the illumination seen by the objects is the same as the one observed from Earth. Only an anisotropic ISRF was considered in this study, since its presence is essential to be able to see scattered light in excess of the background field (Appendix A).

The Galaxy, and especially the Galactic center, is the main source of the anisotropic ISRF that illuminates the clouds. The zodiacal-subtracted mission average (ZSMA) DIRBE survey directly provides the 3D all-sky maps⁵⁶ with an accurate estimate of the sky flux at J, K, 3.5 μm, 4.9 μm and 12 μm wavelengths (and up to 240 μm, which is beyond our scope here, Hauser et al. 1998). We used HEALPix maps with parameter N_SIDE = 256, giving sky pixels with a size of 13.7′, which is smaller than the DIRBE resolution of ~40′ (Hauser et al. 1998). At NIR wavelengths the radiation field was directly obtained thanks to filter conversion from DIRBE to 2MASS (Levenson et al. 2007): DIRBE 1 (J) was divided by 0.97 and DIRBE 2 (K) by 0.88 to obtain the NIR radiation map inputs. The ISRF at the MIR wavelengths was obtained by rescaling the DIRBE 3 (3.5 μm), 4 (4.9 μm) and 5 (12 μm) map fluxes to Spitzer fluxes (at 3.6, 4.5 and 8 μm, respectively) thanks to filter-color corrections and wavelength-scaling deduced from the Galactic spectrum by Flagey et al. (2006, their Table 2). To deduce the illumination at 5.8 μm, we took the 0.3 observed color value R_{5.8 / 8.0} for GLIMPSE (Benjamin et al. 2003; Flagey et al. 2006), consistent with Li & Draine (2001) and not dependent on the line of sight, which assumes that the stellar contribution at this wavelength is completely negligible compared with the diffuse contribution. This presumes that we are dominated by PAH emission and there is no offset between bands (Fig. 9 in Flagey et al. 2006).

3.2.2. Cloud background field

To illuminate the cloud, it is mandatory to consider all the contributions (both stellar I^∗(λ) and diffuse I_diff(λ)). On the other hand, the observed signal is the combination of radiation coming from behind the cloud, attenuated on its way through the cloud, and of the fraction of the radiation field that is scattered toward Earth. Therefore we have to be careful about how we evaluate this background value, I_back. The DIRBE resolution implies that the stellar contribution in a particular beam is always present and wavelength dependent: $\begin{array}{ccc}& & {\mathit{I}}_{\mathrm{back}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{diff}\mathit{,}\mathrm{\left(}\mathit{l,b}\mathrm{\right)}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{\ast }\mathit{bg}\mathrm{+}\mathrm{CIRB}\\ & & {\mathit{I}}_{\mathrm{diff}\mathit{,}\mathrm{\left(}\mathit{l,b}\mathrm{\right)}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathrm{DIRBE}}_{\mathrm{\left(}\mathit{l,b}\mathrm{\right)}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{-}{\mathit{I}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{,}\end{array}$where bg is the fraction of the diffuse light on the line of sight coming from behind the cloud, and CIRB is the cosmic infrared background due to unresolved galaxies from the early Universe (e.g. Gorjian et al. 2000; Levenson et al. 2007). To evaluate this stellar contribution we considered two different approaches. In the first method, we made a sky-direction independent guess on the proportion of the stellar contribution relative to the diffuse contribution in each DIRBE filter. We then used the J/K band ratio to deduce the extinction on the line of sight, and evaluated the stellar contribution in the 3.5 and 4.9 μm filters by using stellar color ratios and the extinction previously deduced (see Bernard et al. 1994 and Appendix B for details). This method is powerful because it gives an all-sky map of the stellar contribution at each wavelength. Nevertheless, it implies the major hypothesis that the signals in DIRBE 1 and DIRBE 2 bands are only due to the stellar contribution. Even though the relative contribution of the diffuse emission is weak in the J and K bands (10–20%), the previous hypothesis is too strong since we are also interested in this residual value. Despite this limitation, this method gives a good approximation at 3.6 and 4.5 μm of the diffuse cloud background field intensity.

The second method is based on the flux subtraction of the sources by summing them from 2MASS and WISE point source catalogs (Skrutskie et al. 2006; Cutri et al. 2012a,b). This has been done in several studies that tried to estimate the CIRB, and we basically applied the same method as described in Levenson et al. (2007). The DIRBE HEALPix maps are made of pixels of equal area but of different shape (healpixels, Górski et al. 2005). Each healpixel contains the average of all pointings that fall inside that pixel for all orientations, because DIRBE observed each position in the sky many times with different position angles. The probability that a star contributes to the measured flux of a healpixel is therefore not trivial to evaluate (Levenson et al. 2007). For each cloud direction, the healpixel shape is different, and we have to compute the probability map of the stellar contributions for each of them separately. We started from the probability that a star falls inside the DIRBE pixel for any position angle (this probability is 1 inside a 20′ radius and 0 outside a 28′ radius) and added the DIRBE pixel probability map for all positions inside the healpixel shape (discretized to the arcmin level). The resulting weight map was renormalized to unity in its center. This map is displayed in Fig. 4 for the Taurus-Perseus region. This weight mask was applied to a catalog of point sources to compute their contribution to the healpixel.

For each cloud of interest, we retrieved the DIRBE pixel value that covers the cloud and its eight neighbors, and the shape of that particular healpixel. We computed the weight map as described above and retrieved the 2MASS and WISE point-source catalogs of the corresponding region (from the Vizier database). The point-source fluxes were obtained by using the zero-flux reference of Jarrett et al. (2011) for WISE and of Cohen et al. (2003) for 2MASS. Then we converted for each band the individual stellar contribution to Spitzer fluxes or kept it in 2MASS fluxes, summed with the appropriate weight maps and subtracted from the DIRBE pixel value (scaled to Spitzer or 2MASS fluxes). The number of point sources considered varies from 5 000 in the 2MASS catalog at high Galactic latitude (36°, L183) up to 30 000 sources in the WISE catalog for lower latitudes (−14°, Chameleon). An example measurement using nine pixels along with a point-source contribution estimate is shown in Fig. 5.

The results for the four regions are summarized in Table 2. For each line of sight and wavelength, we give the value of the central pixel, the average value of the nine pixels, and the dispersion. The I_diff values obtained from individual lines of sight were chosen to be representative of the complete region for the Taurus-Perseus complex, the southern part of Chameleon, and the L183 direction. The Cepheus region is more heterogeneous. It shows a gradient in the direction of the Galactic plane and has to be considered more cautiously (Fig. 6).

The CIRB contribution (Eq. (1)) is negligible compared to I_diff in J, but is in the 50–100% range in K. It represents ~25–50% of the diffuse flux at 3.6 μm and 4.5 μm. At 5.8 and 8 μm, the star contribution is considered to be negligible and the background emission is directly measured from the interpolated DIRBE fluxes in these bands, as explained in Sect. 3.2.1.

The method is efficient but can become inaccurate when very bright stars enter the field because their flux is not always correctly estimated in the WISE or 2MASS catalogs, although it represents a major portion of the total flux in the DIRBE pixel. For a discussion about nonlinearities and saturation caused by bright stars in determining the photometry, we refer to the explanatory supplements of the WISE and 2MASS missions⁷. The adopted correction also assumes that the completeness of the input catalogs is similar to the depth to which point sources are removed in our analysis of the diffuse signal. Indeed, the faint end of the source distribution is only a few kJy sr^-1 in J and K and even less at coreshine wavelengths (Levenson et al. 2007). The uncertainty on the method can be estimated from the fluctuations between the nine adjacent pixels and from the fact that no pixel should show a flux lower than the CIRB flux (Table 2). The zodiacal light subtraction we used is also challenged by some authors, and the result slightly depends on the adopted zodiacal light correction (see Table 5 in Levenson et al. 2007).

The two methods agree well qualitatively at 3.5 and 4.9 μm, but not in J and K bands. While the first method assumes that all the J band flux is due to the stars, the second method finds a sizable fraction of the flux to be due to the diffuse light (10 to 25%). This is expected since the standard interstellar grains in the diffuse medium are more efficient at scattering light at 1.2 μm than at 2.2 μm. The diffuse contribution reaches a minimum in the K band. The diffuse emission fraction is also higher in the 3.5 and 4.9 μm bands since the stellar contribution decreases with increasing wavelength. When we obtained results below the CIRB intensity, we replaced I_diff by the CIRB intensity in our models. This occurred only in the Cepheus direction, for the K, 3.6, and 4.5 μm bands. The second method is more reliable for deriving the measurement of I_diff(λ) when we investigate a particular direction (modeling a cloud), while the first is more useful in exploring what occurs at Galactic scales in the MIR.

The measurement of I_diff(λ) is not the only parameter we need to know to reconstruct the absorption part of the final map; the portion of the diffuse emission that is in front of our cloud (foreground = fg) and the one that is behind (background = bg) is a key element as well (see Eq. (1)). These quantities are directly linked to the portion of dust in the diffuse medium that is located in front and behind the cloud. They are evaluated based on a model of dust in our Galaxy up to 300 pc by Lallement et al. (2014, bg1 and has to be taken as an indication rather than as a precise value. That is why, for certain lines of sight, and when possible, we also used another estimate based on a different method by Marshall et al. (2006) to evaluate the uncertainty on the fg and bg values. This method, known to be unreliable below one kpc, has been refined and gives better approximations at smaller heliocentric distances (Marshall et al., in prep., bg2).

At the end, it is the product of I_diff(λ) with the bg value that gives I_back (Eq. (1)). The contrast level of the emergent flux (I_final) in final model maps will depend on this I_back value as follows: ${\mathit{I}}_{\mathrm{final}}\mathrm{=}\mathrm{(}{\mathit{I}}_{\mathrm{sca}}\mathrm{+}{\mathit{I}}_{\mathrm{back}}\mathrm{\ast }\exp \mathrm{(}\mathrm{-}\mathit{\tau }\mathrm{)}\mathrm{)}\mathrm{-}{\mathit{I}}_{\mathrm{back}}\mathit{,}$(3)}with I_sca the scattered-light image, τ the integrated extinction opacity map, both obtained from the radiative transfer code, and (I_sca + I_back ∗ exp(−τ)) the transmitted signal from the cloud. The focus on four different regions allows us to explore different (I_diff, bg) pair values and see the impact on the emerging intensity in the modeling.

3.3. Molecular cloud content

3.3.1. Cloud model

Because modeling in three dimensions is mandatory (Appendix A), we chose to use a general 3D shape that represents most of the data clouds. An inclined ellipsoid was taken as the cloud model. This cloud model was built to correspond to 40 pixels for a cloud radius (core+envelope) equal to 16 000 AU (Fig. 7). This corresponds to a resolution ranging from 2′′ to 7′′ depending on the distance of the region, a few times the Spitzer resolution. Because the coreshine varies slowly throughout the cloud (red area Fig. 3) and because we are interested in the surface brightness (in MJy sr^-1), this resolution is not crucial. We made a compromise by limiting the total number of cells (with a size of 104 × 104 × 104 cells) to reduce the computational time. However, we kept enough cells to dedicate reasonable physical space to the cloud in our modeling cube while keeping some room for the external part. This also explains why we chose a jump at lower density to the envelope instead of a smooth variation. This has a minor impact on the modeling.

The minor axis has been designed to be half the major axis, and the inclination angle is 60° to see a possible gradient effect as a function of the direction of the Galactic center. The adopted density profile for the core is a Plummer profile adapted to an ellipsoid (from Doty et al. 2005 and Whitworth & Ward-Thompson 2001),$\mathit{n}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}{\mathrm{[}{\frac{\mathit{r}}{{\mathit{a}}_{\mathrm{0}}}}^{\mathrm{]}}}^{\mathit{\alpha }}}\mathit{,}$(4)with$\mathit{r}\mathrm{=}\sqrt{\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{y}}^{\mathrm{2}}}{{\mathit{b}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}}\mathit{,}$(5)where n is the number density, n₀ is the reference density, a₀ the radius associated with $\frac{{\mathit{n}}_{\mathrm{0}}}{\mathrm{2}}$, (x,y,z) defines the position, and (a,b,c) are parameters specifying the shape of the ellipsoid, with an index of α = 2.5. Generally used on a sphere, this type of profile (Fig. 7) is common and realistic for simple molecular clouds like L1544 in the Taurus-Perseus region (Doty et al. 2005).

While the detailed shape does not really matter for our toy cloud model, the column density has to be more representative of the range of our observations. We defined two models, one with a central density of n₀(H₂) = 2 × 10⁶ cm^-3 (M_high = 1.5 M_⊙) and another one with n₀(H₂) = 5 × 10⁵ cm^-3 (M_low = 0.4 M_⊙), which give a peak column density of 9.2 × 10²² cm^-2 and 2.3 × 10²² cm^-2, and we assumed a gas-to-dust ratio of 133 (Compiègne et al. 2011). The 3D cube model of more than one million cells is divided into two regions: an envelope and a core. The envelope, of constant density (1000 cm^-3), is always filled with a standard diffuse grain-size distribution since NIR studies seem to be able to reproduce scattering in the outer parts of the investigated clouds with such distributions (Juvela et al. 2012). The core itself is divided into two parts of approximately equal thickness (Fig. 7) in which different grain characteristics are tested. The threshold between the two core layers is situated at a density of 2 × 10⁵ cm^-3 (for M_high) and 5 × 10⁴ cm^-3 (for M_low). Grain property gradients can thus be investigated from the envelope standard distribution to the inner core evolved grains.

3.3.2. Dust properties

The main dust models used in this paper are based on an extrapolation of standard grains able to fit the observations in the diffuse medium (comprised of a mixture of silicates and carbonaceous grains, in the classical proportion of 3/4 to 1/4, respectively) from DustEm. DustEm is a software package implemented by Compiègne et al. (2011) that is able to compute the dust emission from a spherical grain-size distribution and grain optical properties (obtained internally by using the Mie theory and the IDL code DustProp). Since we focused our study on scattering and absorption and not on emission, DustEm was used as a tool to average the grain properties on their size distribution wavelength by wavelength. In each cloud layer (see Sect. 3.3.1), we mixed the two dust species and adjusted their size distributions independently. This hypothesis is justified because the two species are assumed to be able to coagulate in separate ways (deduced from the extinction curves, Hirashita & Voshchinnikov 2014).

Our extrapolation started from the biggest grain-size distributions of DustEm (hereafter aSil and CBx2 according to the notations in Compiègne et al. 2011, Table 3). The goal is to investigate to what extent grain growth inside molecular clouds, especially coagulation, is able to explain the coreshine observations. For both species, 1 μm grains are known to be sufficiently efficient in scattering and producing coreshine (Pagani et al. 2010 and Fig. 8) but some questions remain: is the suppression of the smallest grains required to observe coreshine? Is there a size limitation for the distribution? What is the influence of the change of the slope of the power-law or of the complete dust distribution itself? What is the impact of taking into account ice mantles and fluffy grains?

In this perspective, we tried different short cutoffs (smallest size for the grain-size distribution: a_min) and long exponential cutoffs (e^{− (a/a_cut − 1)2}, where a is the grain radius, a_cut the cutoff radius) which were applied to the original size-distribution power-law (Compiègne et al. 2011, Table 3). In particular, the index of the original power-law is −3.4 for the silicates and −2.8 for the carbonates. The a_min value evolves from 4 nm (aSil, CBx2) up to 50 nm (C50, S50) while the high cutoff value changes as well. Nevertheless, we also chose to compare grain-size distributions with the same high cutoffs and different short cutoffs to test the influence of small grains (e.g., from Cx1m to C50). All the grain-size distributions between CBx2 and Cx1m, as well as aSil and aS1m, were also used to reproduce a smooth increase in size throughout the cloud (GRAD model hereafter). The change of the whole distribution itself was investigated through different grain types from the literature (Weingartner & Draine 2001; Ormel et al. 2009), as was the porosity (Ysard et al. 2012) and the fractal dimension (Min et al., in prep.).

For all the grain models, the optical properties were averaged in each cell of the cloud. DustEm provides averaged parameters for each grain-size distribution, taking into account the grain-distribution law and the dust mass for each species. Therefore, we considered three averaged parameters: the mean scattering efficiency ⟨ Q_sca ⟩, the mean absorption efficiency ⟨ Q_abs ⟩, and the mean phase asymmetry factor ⟨ g ⟩. This phase function, the probability for a photon packet to be scattered in a given direction, is well-approximated by the Henyey-Greenstein expression ${\mathit{P}}_{\mathrm{HG}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{1}\mathrm{-}\mathrm{⟨}\mathit{g}{\mathrm{⟩}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\mathrm{⟨}\mathit{g}\mathrm{⟩}\mathit{\mu }\mathrm{+}\mathrm{⟨}\mathit{g}{\mathrm{⟩}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(6)with $\begin{array}{ccc}\mathit{\mu }\mathrm{=}\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}& & \end{array}$and $\begin{array}{ccc}{\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{P}}_{\mathrm{HG}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathrm{=}\mathrm{1}\mathit{,}& & \end{array}$which is accurate enough for spherical grains up to a form factor x = 2π ⟨ a ⟩ /λ equal to 10, with ⟨ a ⟩ the mean size of the grain-size distribution and λ the wavelength. In this paper, the form factors of the size-averaged dust models vary from 0.04 to 0.4 at 3.6 μm (Table 3). Only the axial backward scattering might not be taken sufficiently into account. The full radiative transfer equations are described in Steinacker et al. (in prep.). Nevertheless, the averaging process on the whole size distribution was tested by comparing a full distribution and the same distribution cut into three parts and put together in the same cell. No difference was found when the size-distribution discretization had been increased up to 600 steps in DustEm.

To understand the competition between scattering and absorption at NIR and MIR wavelengths, it is important to examine the individual properties for the different sizes and grain components (Fig. 8). Figure 8 shows that for a given total dust mass, ⟨ Q_sca ⟩ increases faster than ⟨ Q_abs ⟩ when ⟨ a ⟩ increases. This is undoubtedly linked to the emergence of coreshine in dark clouds (Steinacker et al. 2010). However, real clouds present a mixture of different grain types and grain sizes, which implies some degeneracy. They are partly removed by modeling the cloud at several discrete wavelengths (from J band to 8 μm).

Because the grains of Compiègne et al. (2011) are not the only ones able to explain the observations in the diffuse medium, and to explore beyond spherical grains, we added a sample of dust grain varieties. In our grid of models, silicate and graphite mixtures (Weingartner & Draine 2001, WD31, WD55, WD55B), porous grains without ices (Ysard et al. 2012, YSA 0%, YSA 10%, YSA 20%, YSA 40%), monomer fractal aggregrates with different sizes (Min et al., in prep.⁸, MIN0.2, MIN0.8, MIN1.2, MIN2.4, MIN4.0) or compact agglomerates with ices (Ormel et al. 2009, ORI2, ORI3) and their counterpart without ices (Ormel et al. 2009, ORNI2) were included with different maximum sizes (see Table 3 and references for more details). We did not attempt to provide a review of all the grain types available in the literature, among which a–C:H grains (Jones 2013), and iron inclusions in silicates (Jaeger et al., in prep.⁹) are other possibilities. The goal was to open the scope of what types of grains will be suitable to explain our observations.

4. Results

The coreshine strength is linked to the grain properties and to the environmental conditions. The modeling has to account for the absolute coreshine intensity, which is a contrast problem and a competition between scattering and absorption. These extinction processes are linked to the dust properties, and in particular to the grain-size distribution. However, these properties are degenerate when considering a single wavelength. Therefore, we adopted a multiwavelength approach in which the albedo variation between wavelengths, dependent on the dust population itself, is well-characterized by intensity ratio measurements.

In Sect. 4.1, we present the measurements of the two key quantities, the 3.6 μm intensity and the coreshine ratio (between 4.5 and 3.6 μm wavelengths) for the observational sample of the four regions. In Sect. 4.2, we expose the conditions necessary for coreshine to occur and retrieve the correct absolute intensity from modeling, while in Sect. 4.3, we detail the use of the different intensity ratios as a tool.

4.1. Surface brightness and coreshine ratio in observations

We investigated two thirds of the detected coreshine cases, located in four different directions of the Galaxy, to search for either spatial or physical regional variations. The two main regions we ignored (Orion/Monoceros and ρ Oph) present an equal number of detections and nondetections, which indicates that these regions might be more sensitive to specific local conditions and are beyond the scope of this work. Figure 9 gives the 4.5 μm/3.6 μm ratio as a function of the 3.6 μm intensity for all the sources (values are given in Table 1). The horizontal bars represent a ±10% uncertainty that is a conservative value of the error budget, which is dominated by the background removal, while the vertical bars represent the range of ratios around the peak ratio value (Sect. 2.3 and Fig. 3). Some of the weakest sources (21) have a signal-to-noise ratio at 4.5 μm that is too low to safely estimate the coreshine ratio. Their 3.6 μm intensity ranges from 10 to 40 kJy sr^-1. For the others, the 3.6 μm intensity ranges from 21 to 390 kJy sr^-1, and their ratio from 30 to 110%. Starless cores are in a narrower range, 10 to 100 kJy sr^-1 and ratios from 35 to 51% (the L1517C case – 53% – is weak and noisy and probably not different from its neighbors, L1517 A and B – 51%). This upper limit of ~50% is well-explained by the fact that the incoming DIRBE illumination ratio is only 70% at 4.5 μm with respect to 3.6 μm and that the scattering efficiency of grains up to 1 μm is always lower at 4.5 μm than at 3.6 μm (Fig. 8). In contrast, the stronger and redder coreshine flux of many sources with local or embedded young stellar objects (YSOs) is clearly linked to the YSOs themselves, which provide more photons than the ISRF locally, with a much redder color because of their dust cocoon.

However, some of the sources with nearby or embedded protostars may remain similar to the starless cores. There are three possible reasons: the source is i) too weak (like the VeLLO – very low luminosity object – in L1521F); ii) too deeply buried (like the strong jet driver L1157–mm/IRAS 20386+6751); or iii) too far outside the cloud, although this is difficult to estimate since the relative position of the cloud and the YSO are not known with enough precision along the line of sight (see Table 1 and references therein). This is the case for the L1157 region, which shows a low ratio (41%) near the embedded driving source of its strong outflow, and partly along the southern outflow itself, while it displays a high ratio (71%) just beyond the tip of this southern outflow where no YSO can be seen. A possible explanation is that the outflow provokes a shock region (B1 seen in H₂) that emits in the MIR with a 4.5 to 3.6 μm ratio of 2.5 that might change the coreshine ratio of that part of the cloud.

For several sources the embedded object has a limited range of influence, and if the core is extended, a multimodal ratio distribution appears (Fig. 10). In the cloud (G303.72-14.86) presented in Fig. 10, the strongest source is at the top end of the cloud, inside the blue polygon, and two other interacting but weaker stars are situated just below, inside the red polygon. The sources themselves are masked, hence the white spots. The situation for YSOs is quite complex, because other physical effects (e.g., dust alignment due to polarization – Chapman et al. 2013 – and/or dust destruction) could have a profound impact on the coreshine effect. Indeed, for L1152 and L1228 there is a suggestion from NIR (JHKs) and MIR (Spitzer) that there are regions where shocks have destroyed grains (Chapman et al. 2009), which would potentially affect the coreshine effect. This is probed by the fact that even though the radiation is enhanced by the local source, the 5.8 μm band remains in absorption, which would not be possible with too large grains.

Although the scattering effect does not seem to change much across the Galaxy, as shown in Pagani et al. (2010, their Fig. 4), we expected to find some difference between the four regions when taking into account the background and the contrast problem (see next section). Of course, Pagani et al. (2010) indicated scattering effects toward the Galactic center, while we show below that it is impossible to see any, but Fig. 9 clearly shows that away from the Galactic plane and bulge, there is no difference between the regions either in absolute flux or in ratio values, while the background values I_back vary by a factor of up to around four between Cepheus and Taurus-Perseus (Table 2). Moreover, the main incoming radiation field, from the Galactic center, is either behind the cloud, on its side, or behind the observer (in the Taurus-Perseus direction, which implies some efficient backward-scattering for anticenter cores, Steinacker et al. 2014a). The dispersion seems to be dominated by individual cloud properties instead, with little sensitivity on the coreshine ratio.

Finally, the very high proportion of coreshine cases in these regions is an interesting question. Either the grown grains are already present in the local diffuse medium of these regions and become apparent when the dust is concentrated enough, or the grains coagulate relatively fast after the clouds are formed and before the turbulence dissipates in the cores (Steinacker et al. 2014b). In the first case, the presence of the grown grains in the diffuse medium might explain the excess of emission at 3.6 and 4.5 μm seen by Flagey et al. (2006). In principle, this could be verified by examining the differential MIR spectrum on lines of sight that sample the diffuse medium inside and outside the Gum/Vela region for which the low number of cores with coreshine and the high presence of PAHs may indicate the erosion of large grains by the blast wave (Pagani et al. 2012). The second hypothesis would require the modeling of grain growth together with cloud contraction and turbulence dissipation.

4.2. Coreshine emergence

To study the coreshine emergence as a function of the galactic position, we used the all-sky background map obtained after stellar subtraction (following the first method, see Sect. 3.2.2), and we modeled this for several elevations toward the Galactic center direction, 90° from it (in longitude) and toward the Galactic anticenter direction (Fig. 11) for a given cloud (M_low) and grain-size distribution (aS1m/Cx1m) models. We easily explain that it is not possible to observe coreshine in the Galactic plane because of the strong background field, even with favorable grain properties (grain-size distribution up to 1 μm grains and ⟨ Q_sca ⟩ greater than ⟨ Q_abs ⟩). For a Galactic longitude (l) of 0°, the bulge dominates the background field so strongly that coreshine is not able to appear until a galactic latitude (b) around 10°. In contrast, in the anticenter direction coreshine is able to appear rapidly with elevation, as low as b ≥ 3°. When the main illumination field comes from behind, at small angles, coreshine appears for higher background values than for other illumination directions (~300 kJy sr^-1 instead of ~120 kJy sr^-1, Fig. 11). The values are compatible to what was found by Steinacker et al. (2014a); here we also compare these criteria directly with the detections. This modeling correctly explains the observations (Fig. 1), but the b values given in Fig. 11 have to be taken qualitatively since the values can change with the presence of local sources, the grain properties, and the background portion estimates. Finally, of all the clouds that contain enough micron-size grains to efficiently scatter the MIR light, only those outside the Galactic plane are detectable via their coreshine emission, and the 50% detection of positive coreshine cases (Pagani et al. 2010; Paladini et al., in prep.) is only a lower limit in terms of grown grains in clouds¹⁰.

In addition to the Galactic plane contrast problem, clouds will also appear in absorption if ⟨ Q_sca ⟩ is lower than ⟨ Q_abs ⟩. Nevertheless, there are not many cases in real absorption outside of the Galactic plane in our data collection (Table C.1). Scattering is always present in dark clouds, even when the cloud is seen in absorption. The coreshine phenomenon starts to appear when the scattering signal is able to exactly compensate for the extinction of the background field. Incidentally, clouds would totally disappear if this equilibrium was reached at MIR wavelengths. The coreshine phenomenon can be understood either by the presence of large grains (~0.5–1 μm) or a stronger local radiation field that would only enhance I_sca (Eq. (3)). The enhancement of the local radiation field has to be considered cautiously since it has to be consistent with FIR emission of the clouds. Here, we are not referring to embedded YSOs, but to the large-scale local radiation field that surrounds starless cores. Incidentally, Evans et al. (2001) showed that for a selection of prestellar cores in Taurus and ρ Oph regions, FIR observations seem to require a lower ISRF than standard to be fitted. In contrast, for the ρ Oph region, Roy et al. (2014) argued that the ISRF is one order of magnitude higher than standard.

At 5.8 and 8 μm, all the clouds of the cold Spitzer survey (Pagani et al. 2010) appear to be in absorption. While at 8 μm, this is partly due to the 9.7 μm silicate absorption feature wing that intercepts half of the 8 μm filter width, at 5.8 μm, it is only due to the background field strength. Indeed, the diffuse part I_diff increases globally more and more with the wavelength (Table 2) in the IRAC bands range, as can be seen in Flagey et al. (2006). In parallel to this increase of the diffuse background field, the scattering efficiencies drop. This provides additional constraints, in particular, our modeling can eliminate the grain models for which the 5.8 μm map shows emission. This is especially true for distributions that include grains in the range 2–5 μm for both silicates and carbonaceous grains (aS5m/Cx5m) and places an upper limit for the grain size based on the 5.8 μm diagnostic deep inside the core. Nevertheless, this limit in size has to be taken cautiously since it is dependent on the grain-size distribution law and on the grain composition. It has to be understood as un upper limit on the abundance of grains larger than 2 μm, and the WD55B distribution, for example, remains compatible with absorption at 5.8 μm. The 8 μm map always appears in absorption in our simulations and is not, as expected, sensitive to the grain properties, but could be used to constrain the silicate column density.

4.3. Intensity ratios as discriminants

We aim to separate the grain models according to their ability to explain the observed intensity ratios between bands. We start from the coreshine ratio (4.5/3.6) and extend the method to a comparison between NIR and MIR bands (K/3.6) and to NIR ratio (J/K ratio). The four modeled directions cover a range of local properties like the Galactic center direction and the differential local background field R_{λ1 /λ2}(I_back) = I_back(λ₁) /I_back(λ₂) (Table 2).

To measure the ratios, we took the median value inside a polygon on the simulated map ratios. We obtained the model uncertainty by comparing these ratios computed with different I_back estimates (obtained with the second method, Sect. 3.2.2), namely the center-pixel values versus the nine-pixel average values, and the two bg estimates (bg1,bg2) (Table 2). The suitable grain models are those that satisfy all the criteria, taking into account the uncertainties on the background variation. For the M_low cloud model (see Sect. 3.3.1), the chosen polygon is very central and corresponds to a mean value of A_V = 23 mag (assuming R_v = 5.5), while the mean polygon extinction for the M_high model is about A_V = 21 mag, chosen to be between the two layers and outside the internal depression.

4.3.1. Grain properties deduced from the coreshine ratio

We plot the coreshine ratio derived from the model calculations as a function of the 3.6 μm coreshine intensity for all the grain models and separately for the four modeled directions corresponding to the observed regions (Figs. 12 to 15). All the models that do not show coreshine at both 3.6 and 4.5 μm are not displayed on the ratio plots, especially the WD31, ORI2, ORNI2, and YSA models. Figure 12 confirms that the two different cloud models that investigate different zones with close averaged visual extinction give coherent results and that we correctly represent the observational range. Moreover, the plots are not really sensitive to the background portion itself concerning the coreshine ratio range (bg = 95 or 75% for the Taurus-Perseus region, Fig. 12, and bg = 75 or 50% for the L183 direction, Fig. 13). However, this implies some degeneracy in the grain-property solutions for some regions (Fig. 12). Conversely, regions where the I_sca term is dominant (Eq. (3), Fig. 13) give more reliable results. This regime is reached for the L183 and Cepheus regions. For the following, we keep the estimates by Lallement et al. (2014) as a reference (bg1, Table 2).

The first question is to determine whether small grains have any impact on the coreshine modeling. We found that no differences appear between models that contain small grains between 4 and 10 nm and those without grains below 10 nm (from aSil to S10 grains and from CBx2 to C10, Tables 3, 4). Removing the small grains at constant dust mass has been tested up to 50 nm (S50, C50) and the results are very similar to those with the same high cutoff (a_cut = 1μm), but with a starting size of 4 nm (aS1m, Cx1m), especially when taking into account the uncertainties (see, e.g., Fig. 13). The suppression of the smaller grains is therefore not mandatory for explaining the coreshine phenomenon because they neither contribute to nor attenuate the signal significantly.

On the other hand, we tested the large grain-size increase beyond the expected limits (Andersen et al. 2013). The more the mean size of the distribution ⟨ a ⟩ increases, the more the 3.6 μm flux increases as well. The grain models on the right side of the ratio plots (Figs. 12 to 15) have the highest number of big grains, with an exponential cutoff at 1 μm (aS1m, Cx1m), 2 μm (aS2m, Cx2m), or 5 μm (aS5m, Cx5m). However, when the high cutoff of the grain-size distribution increases, the coreshine ratio starts to saturate, and consequently, grain growth alone is not a pertinent answer to increase the ratio beyond 50% for the Taurus-Perseus and L183 regions (Figs. 13 and 14) and 65% for Cepheus and Chameleon regions (Fig. 15). This is well-explained by the combination of three effects. First, the integrated illumination field intensity at 4.5 μm deduced from the scaled DIRBE map is only 70% of the 3.6 μm intensity. Second, an I_back lower at 4.5 μm than at 3.6 μm helps to increase the ratio, as explained in Sect. 4.2. Since the ratio values of I_back at 4.5 and 3.6 μm  R_{4.5 / 3.6}(I_back) remain at around 1 for the Taurus-Perseus region and L183 direction, they have a minor impact on the coreshine ratio, while the lower I_back(4.5 μm) toward Chameleon and Cepheus provides an explanation of why ratios near 65% are reached (Fig. 15). In the third place, the ratio of the scattering and absorption efficiencies (Fig. 8) becomes flatter with wavelength for both grain types when the grain size is increased up to 5 μm.

We also confirm what has been found previously (Pagani et al. 2010; Steinacker et al. 2010). Classical diffuse medium grains (WD31 Weingartner & Draine 2001, aSil and CBx2) are not efficient enough to scatter in the MIR range. Indeed, it is well known that the ratio of visual extinction to reddening, R_V = 3.1, is not valid in dense, cold environments, and R_V = 5.5 has been advocated (Weingartner & Draine 2001). This change of slope has been explained by grain growth. While the WD55 model gives a 3.6 μm coreshine flux between 20 and 50 kJy sr^-1 and a coreshine ratio of about 20%, depending on the local conditions (direction and I_back values), the WD55B model, which includes grains up to 10 μm in size, has been found promising by other observations in the MIR (Ascenso et al. 2013), and allows the coreshine ratio to increase up to 40%.

The previous results were obtained for compact spherical grains, and we also investigated the consequences of the fluffiness and of the coagulation, which has to be understood here as an agglomeration of smaller grains. We tested the evolved grains of Ormel et al. (2009), with or without ices, with different evolution timescales. We aimed to compare three different populations: i) a reference population, without ices, and an evolution timescale of 1 × 10⁵ years (ORNI2), which happened to produce no coreshine in the modeling; ii) another one starting from the same distribution, and which evolved in the presence of ice mantles (ORI2); and iii) one that has grown during a longer timescale of 3 × 10⁵ years with ices (ORI3). Despite the presence of ices, which favor coagulation, the ORI2 grain-size distribution has not had enough time to reach a sufficient size (Fig. 1 in Ormel et al. 2009) to make coreshine appear, as our models confirm. In contrast, the ORI3 model traces some evolution like the WD55B model and shows a ratio compatible with the observations. For these two models, even though the ratio is compatible with the observations, the 3.6 μm flux is higher than observed, which is acceptable since we could always adjust the cloud model to obtain the correct flux with an identical ratio.

These comparisons show that the first key role of ices is to favor growth, but the 3.05 μm H₂O ice feature included in the edge of the IRAC filter might also play a role in increasing the coreshine ratio by changing the dust optical properties in the same way as the silicate feature enhances the absorption in the 8 μm filter. Since, without ices, the largest size reachable in Ormel bare-grain models is about 1 μm, similar to the size obtained with the ORI2 model, which produces no coreshine, it is not possible to compare the ORI3 model to the same size distribution without ices to separate the ice effect from the pure-growth effect.

Scattering optical properties are supposed to be more sensitive to grain surface variations than absorption efficiencies. Consequently, we expect to see a direct impact of the fluffiness on the coreshine. To calculate optical properties for fluffy grains, one can focus on the porosity with a simple approach where the portion of silicates, carbonaceous grains, and vacuum in the grains do not depend on their sizes (Ysard et al. 2012), or one can explore the fractal dimension, which means to focus on the asymmetry of the grains themselves, built by agglomeration of monomers (Min et al. 2006). First, we studied the influence of the degree of porosity. The grains from Ysard et al. (2012) that were able to explain the FIR emission in L1506C (Ysard et al. 2013), a condensation in a Taurus filament, are not efficient enough to produce coreshine at 4.5 μm, while the 3.6 μm coreshine flux, about 13 kJy sr^-1 (compared with 33 kJy sr^-1, Table 1), does not vary from 0% porosity to 40% porosity (YSA models). Because porosity does not seem to be an answer to change the 3.6 μm coreshine intensity, and correspondingly the coreshine ratio, we examined the influence of fractal structure. As preliminary results, we tested fractal aggregates that show optical properties quite different from the compact spherical grains (Min et al. 2006). These fluffy monomer aggregates behave like small spherical silicate particles for the 9.7 μm silicate feature (Min et al., in prep.). Their sizes vary from 0.2 μm to 4 μm¹¹.

We had to approximate the real phase function provided by Min et al. by an equivalent asymmetry coefficient (⟨ g ⟩, Eq. (6)). As we discuss below, the phase function is not expected to modify the coreshine ratio by much, especially because this approximation applies to both wavelengths. Indeed, the two coreshine wavelengths are close enough so that the respective variation in the phase function is weak. These grains are really promising since they do not necessarily imply a significant grain growth to produce a ratio about 40 to 50% and could be an interesting answer for clouds, which would not show much evolution from other tracers such as depletion (e.g., L1521E Tafalla & Santiago 2004; Pagani et al. 2010). Furthermore, we find the same saturation of the ratio for the agglomerates as for pure spherical grains with grain growth, which seems to confirm a maximum size efficiency, as can also be deduced from Fig. 8. Fractal aggregates behave like big spherical grains to produce a high coreshine ratio.

To limit the number of models, we restrained ourselves to two layers for the part dedicated to the core (Sect. 3.3.1), that is why we did not aim to fit the observations exactly. We introduced a finer slicing in ten layers inside the core, filled with grains with a high cutoff that increased from 0.2 μm to 1 μm (GRAD model). We obtained a ratio close to a two-layer core model with a_cut of 1 μm (aS1m and Cx1m mixed), but with less flux (see, e.g., Fig. 14), which is expected for the observations. This again emphasizes that the cloud model and the way we fill it with several layers would act only on the absolute coreshine fluxes and would not modify the ratio.

The phase function can slightly change the ratio, especially if backward-scattering becomes a non-negligible quantity. This is the case for the Taurus-Perseus region, for which the main anisotropic radiation source (the Galactic center) comes from behind us. Nevertheless, even if we consider an extreme case with ⟨ g ⟩ = −0.99, which corresponds to a probability of 75% to have axial backward-scattering, the ratio in the Taurus-Perseus region, with the same ⟨ Q_sca ⟩ and ⟨ Q_abs ⟩ absolute values, increases from 35 to 40% for grains up to 1 μm (aS1m, Cx1m). The variation is even smaller for less extreme probability changes. In the most extreme case, even if the ratio is not really sensitive to the phase function modification, the absolute intensity at 3.6 μm has been enhanced by almost a factor 30 to reach a few hundred kJy sr^-1. Incidentally, the emergence of coreshine for some dust models will be dependent on this parameter. We estimate that the Henyey-Greenstein approximation is satisfactory for most of the directions since the ratio is only slightly sensitive to phase function variations. We will investigate the impact on the absolute flux by considering the true phase function when we perform a full cloud model, especially for regions in the regime where the influence of I_back is dominant (Sect. 4.2).

4.3.2. Impact of a local source

Since the presence of an embedded object is often found to enhance the coreshine ratio and the coreshine intensity in its vicinity (Table 1, Figs. 9 and 10), we explored the influence of a Class 0 or a Class I object in our model. The full treatment of a protostar reddened by its compact and dense dust cocoon is beyond the scope of this work. To investigate its possible effects qualitatively, we inserted blackbodies (stars) in the model with temperatures and fluxes typical of Class 0 or I objects at coreshine wavelengths. We ran the test with three different temperatures, 400, 950, and 2150 K, and adjusted the flux accordingly. At 400 K, the 4.5/3.6 ratio is 3.8 and is representative of a Class 0 case, as was observed for IRAM 04191 (Chen et al. 2012). At 950 K, it is almost flat (1.2), and we obtain a flux and a coreshine ratio similar to what could be deduced for a Class 0/I object like IRAS 04016+2610 or IRAS 04361+2547 (Robitaille et al. 2007). At 2150 K, it is typical of a solar type Class I protostar and its 4.5/3.6 ratio is 0.8, like IRAS 04325+2402 (Robitaille et al. 2007), which is close to the DIRBE all-sky ratio of 0.7. The class distinction is not very important here since it depends very much on the geometry of the source and other factors that are beyond the scope of this paper. Several papers discuss these classification problems (e.g., Robitaille et al. 2007; Enoch et al. 2009; Kirk et al. 2009).

For the densest cloud model (M_high), we introduced these three different stars at the edge of the cloud.The test with the Class I object (2150 K) shows a small enhancement of the coreshine ratio in the vicinity of the object (Fig. 16). With the Class 0/I object at 950 K, the ratio is increased to ~0.8. Finally, with the Class 0 object (400 K), depending on its actual flux, the ratio can vary from below 0.5 for a negligible contribution (as in L1157 near the embedded object) to ≥2 if the local source flux dominates. We can therefore clearly reproduce the trend of the observations. However, we are not trying to make exact fits, and to proceed with the modeling, one has to constrain the local radiation field that emanates from the embedded source to avoid the degeneracy between the grain properties and the local source flux and color.

Indeed, the question arises whether an embedded protostar can be helpful in modeling the coreshine effect. If the protostar is weak and deeply embedded, its properties will be difficult to assess and its impact on the cloud scattered-light will only add another degree of freedom. Similarly, if the source is clearly outside of the core, its distance to the cloud along the line of sight will remain a free parameter. If the YSO is clearly embedded in the cloud, but not too deeply, or in contact with it, so that its contribution can be measured directly, little uncertainty is added. This has three more advantages: first, the ISRF becomes dominated by the local source and is therefore better constrained; second, as the scattered flux will be higher, the error due to the background uncertainty in the modeling becomes smaller; and third, the 5.8 μm criterion becomes even more constraining in enhanced ISRF conditions. Indeed, a stronger ISRF increases the capability to produce coreshine (either smaller grains can shine or big grains shine more). In a local stronger field, the nondetection of emission at 5.8 μm therefore places a more stringent constraint on the abundance of grains above 1 μm, as presented in Sect. 4.2. However, the presence of a Class 0 or Class I object indicates a more evolved cloud and, presumably, more evolved grains. Therefore studying starless clouds or protostellar clouds are both important and have caveats.

4.3.3. Extension to the NIR

The information obtained from the coreshine ratio, which is more sensitive to the grain models than to the other free parameters, can be extended to NIR wavelengths. Specific studies on NIR scattering have been done before (Malinen et al. 2014 and references therein). Our approach assumes that the different wavelengths are close enough to investigate the same volume of the cloud and far enough to see a variation in the slope for the different grain types (Fig. 8). Because our NIR observations are limited to a few cases, we only set a range for the expected ratio values, which we admit does not provide much constraint. Nevertheless, we attempted to see whether the multiwavelength approach could lead to a sample of suitable grain models.

First, we studied the K/ 3.6μm (= NIR/MIR) ratio as a function of the coreshine ratio toward L183 and the Taurus-Perseus complex. The observational range obtained for this ratio is between 3.5 and 6 for L183 and 3.2 and 7 for Taurus-Perseus region (Figs. 17 and 18). This eliminates models with a valid coreshine ratio but a NIR/MIR ratio that is too low (ellipse, Fig. 17). In both directions, the selected models correspond to relatively large grain-size distributions, which is what we expected because the K band is assumed to sample a region close to the core (see Fig. 19 and Andersen et al. 2013). Nevertheless, when the grain-size distributions contain too many large grains (aS2m/Cx2m, aS5m/Cx5m), the NIR/MIR ratio becomes too low. Grain growth on spherical grains decreases the NIR/MIR ratio because when ⟨ a ⟩ increases, the K intensity increases, but more slowly than the 3.6 μm coreshine intensity. This mainly applies to clumps with A_V> 10 mag, when the K band is already significantly saturated. In the same way, the fractal aggregates, even the 0.2 μm monomer aggregates, also seem to predict a NIR/MIR ratio that is too low. In contrast, ice mantles (ORI3) give a higher NIR/MIR ratio and are one of the more suitable grains for the Taurus-Perseus region. In both directions, grain size distributions that contain bigger silicates than carbonates are more likely to be able to explain the observations (purple squares – Figs. 17 and 18). Some spatially constant grain-size distributions are also able to explain them, especially those up to 1 μm with or without ices (aS1m/Cx1m, ORI3). Finally, the L183 direction is well explained by silicates up to 1 μm and with different grain-size distributions for the carbonates (from a high cutoff of 0.15 μm, CBx2, to 1 μm, C50) while the Taurus-Perseus direction can tolerate a wider variety of grains as a solution, for example, model A has standard grains in the external core layer and grains up to 1 μm (S50/C50) in the inner layer (Table 4). In this case, only the second layer is visible in coreshine (Fig. 19), and almost the same region emits in K band while the J signal is more extended.

The simple cloud model we used could be a limitation to investigating shorter wavelengths, which are assumed to be more sensitive to the smaller grain part of the distribution contained in the extended envelope. Nevertheless, we obtained interesting results about the J/K ratio (NIR ratio). The typical observational range obtained for the J/K ratio is from 0.3 to 3, and this ratio is highly position-dependent, which is well-reproduced in the modeling (model A, Fig. 19). Globally, it shows a bimodal distribution: one around 3 for the part dominated by the J region of emission, and one below 1 for the central part dominant at the K wavelength. In particular, the model-output median values for this ratio, which averaged this bimodal behavior, varies only from 0.3 to 0.9 (Table 4), for the remaining suitable grains. Therefore it cannot be used as an indicator to differentiate the models and is assumed to be more sensitive to the cloud structure than to the grain properties.

In conclusion, the preliminary constraints set by the multiwavelength approach give us suitable grain models as the solution for the Taurus-Perseus and L183 regions in standard ISRF illumination conditions. L183 ratios are well reproduced by grain models that contain larger silicates or a mix of larger carbonates in one layer and larger silicates in the other layer, while the Taurus-Perseus region also admits spatially constant grain-size distribution models as solutions (Table 4). In particular, grains that might explain the observations in the Taurus-Perseus region are in general larger than those able to explain the L183 observations. The addition of water ices decreases the coreshine ratio but increases the NIR/MIR ratio and could be an answer for some cores in the Taurus-Perseus region, but has to be considered carefully. Larger silicates also have the tendency to increase the NIR/MIR ratio while decreasing the coreshine ratio, and carbonates increase the latter and yield higher flux at 3.6 μm. In any case, grain composition has an important role to play, and the classical spherical silicates/carbonates approach has to be modulated, for example, the fractal dimension is a main actor even if it is not able to explain the NIR/MIR ratio alone for both regions.

5. Conclusion

We explored a radiative transfer model to study the scattering of NIR and MIR light in dark clouds by varying the grain properties and evaluating the impact of other parameters such as the ISRF and the background-field intensities. We compared our large sample of coreshine observations with the modeling and selected a handful of cases to explore the capability of coreshine to separate dust properties. Our results are the following:

• The emergence of coreshine is a contrast problem that has to betreated carefully. The cloud background-field estimation is a keypoint to the modeling. We adapted the Levenson method withcareful flux conversion to the pertinent wavelengths. Moreover,the choice of keeping the standard ISRF is the safest approach andalso a real challenge that will benefit from emission modeling infuture studies.

• Merging previous samples of coreshine observations to build sufficient statistics, some Galactic regions appear to be favored. This could be considered in a larger picture as the presence of previously grown grains in the initial diffuse medium of these individual regions. It opens the scope of modeling coreshine with regard to individual regions or for clouds in close environmental conditions (such as L1517A,B,C, and D, or L134 and L183).

• The use of 3.6 and 4.5 μm coreshine bands and especially their ratio, referred to as coreshine ratio, adds constraints on the grain properties in the core.

• The highest coreshine ratios and fluxes are obtained for clouds that contain an embedded object, probing the influence of local change in the radiation field, which has been tested qualitatively in our modeling.

The main conclusions for the grain properties are:

• Small grains have no influence in the modeling on the coreshineratio.

• For any given coreshine ratio, the absolute 3.6 μm flux value is somewhat adjustable from both the cloud model and the phase function.

• The size increase is mandatory but not sufficient to explain coreshine. The coreshine ratio quantity tends to saturate with the grain size increase, both for spherical compact grains and for agglomerates of monomers. While bare spherical grains show a saturation of this ratio above a high exponential cutoff of 1 μm, fluffiness helps to raise the saturation limit.

• Since dust grains inside molecular clouds are expected to be icy, the role of ice mantles needs to be further investigated. Indeed, we expect ices to favor fluffiness and growth, but in the meantime, the only model we tested shows that it decreases the coreshine ratio. New models of icy grains are urgently needed.

• Both the NIR/MIR ratio for the core outer layer and the absence of emission at 5.8 μm for any layers eliminate a mix of silicates and carbonates that both include grains above ~2 μm in meaningful quantity.

• Porosity has no impact at the studied wavelengths.

• The comparison with other wavelengths for a valid coreshine ratio can help to distinguish between the cloud models and the grain properties.

• When the NIR and MIR do not sample the same zone, we might be able to peel the cloud layer by layer thanks to the multiwavelength approach. If the same region is sampled, the grain-size distribution properties have to be customized to be compatible with the observations at the different wavelengths.

We gathered a database for the grain behavior that covers grain growth by coagulation (up to 5 μm in size). We extended our study to include aggregates and ice mantles and found several promising grain types that are able to reproduce the observations. We constrained the ISRF and the background value for different lines of sight and deduced the impact of the other free parameters on the modeling. The complementarity with NIR observations introduced to highlight the promising perspective of a 3D multiwavelength cloud modeling has to be investigated on a real cloud following all the key points above and extending the study to FIR emission.

Online material

Appendix A: Isotropy versus anisotropy for the incident radiation field

Coreshine (and any types of scattering processes) can only be seen in emission if it follows a number of conditions, as discussed in Sect. 4.2. Another compulsory condition is the anisotropy of the incident ISRF. This anisotropy can come from large scales (such as the Galactic structure) or small scales (such as stars, near or far). To illustrate the importance of anisotropy we considered two simple cases: only isotropic illumination, or the combination of one type of anisotropic source with the isotropic illumination.

If the ISRF is isotropic, there is no preferred direction for photons to travel. To see the cloud in emission, we need to introduce some anisotropy to concentrate photons toward a privileged direction, that is, toward the observer. The cloud would therefore have to produce this anisotropy. To increase the number of photons toward the observer, the cloud would have to act like a telescope mirror pointed at Earth to collect photons from many directions to redirect them in a single direction (this would of course lower the number of photons scattered in some other directions). Since the cloud has no such focusing capabilities and the observer no privileged position, an isotropic ISRF cannot make a cloud glow. In fact, for any direction across the cloud there are as many deviated photons away from that line than there are photons from other directions that are deviated into that propagation line and all the scatterings cancel out. This is in an ideal case without absorption. In presence of absorption, a part of the photons are lost and the cloud can only appear as a darker object against the background sky, never in emission. The field must therefore be anisotropic.

In the anisotropic case, there is one noticeable direction, which is the line going from the localized light source to the cloud (we consider a single source of photons, such as a nearby star, superimposed on an isotropic ISRF). This is the path with the highest number of photons traveling toward the cloud. As in the isotropic case, photons enter from all directions and are partly deviated, partly untouched (and partly absorbed, but we still ignore absorption here). The difference to the isotropic case is that across the cloud along this particular line, there will be

more deflected photons away from it than photons brought back onto it. The opposite can occur for some or all of the other directions (depending on details such as the phase function of the grain-scattering properties). Therefore, away from this path, the number of deviated photons that leave the cloud increases (while it decreases along the path). Eventually, there can be more photons deviated toward the observer from the anisotropic source than photons coming from behind the cloud and deviated away from the observer. The net effect is to show the cloud in emission.

To illustrate this effect we ran a model based on a cube of constant density, tilted at 30° angles on two axes to see the edges and three sides. The cube is either in an isotropic field or a composite of isotropic and anisotropic fields. It is filled with dust with only scattering capability (absorption coefficient is set to zero). Figure A.1 shows the results. The left column shows slightly different images of scattered photons, but taking into account the background illumination absorption due to the cloud opacity (displayed in the central figure), following Eq. (3), the cloud completely disappears in the isotropic field. This is explained by the fact that all the photons scattered toward the observer (as seen in the left panel) are exactly compensated for by the photons scattered away from the line of sight for the radiation field coming from behind the cloud. A close inspection of the image reveals the cube only by the numerical noise of the Monte Carlo radiative transfer code . If the cloud opacity had not been set to zero, the cloud would have been seen in absorption against the background, while for the isotropic+anisotropic case, the cloud would have appeared in emission or in absorption, depending on the balance (Sect. 4.2).

Appendix B: Global method for subtracting the stellar contribution

We assume that the intensity in the DIRBE1 (1.25 μm) and DIRBE2 (2.2 μm) bands is only due to the sum of individual stellar contributions (PSC): ${\mathit{I}}_{\mathit{K}}\mathrm{=}{\mathit{I}}_{\mathit{K}}^{\mathrm{\ast }}\mathrm{=}\begin{array}{c}\sum \end{array}{\mathit{I}}_{\mathrm{PSC}\mathit{,}\mathit{K}}\mathrm{=}\mathit{DIRBE}\mathrm{2}$(B.1)\newpage${\mathit{I}}_{\mathit{J}}\mathrm{=}{\mathit{I}}_{\mathit{J}}^{\mathrm{\ast }}\mathrm{=}\begin{array}{c}\sum \end{array}{\mathit{I}}_{\mathrm{PSC}\mathit{,}\mathit{J}}\mathrm{=}\mathit{DIRBE}\mathrm{1.}$(B.2)

The intrinsic color I_int of the stellar component was measured from high Galactic latitude and low dust emission regions and is also equal to ${\mathit{I}}_{\mathrm{int}}\mathrm{(}{\mathit{\lambda }}_{\mathrm{1}}\mathit{/}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{)}\mathrm{=}\frac{{\mathit{I}}_{\mathit{\lambda }\mathrm{2}}^{\mathrm{\ast }}\mathrm{\times }\exp \mathrm{\left(}{\mathit{\tau }}_{\mathit{\lambda }\mathrm{2}}\mathrm{\right)}}{{\mathit{I}}_{\mathit{\lambda }\mathrm{1}}^{\mathrm{\ast }}\mathrm{\times }\exp \mathrm{\left(}{\mathit{\tau }}_{\mathit{\lambda }\mathrm{1}}\mathrm{\right)}}\mathrm{·}$(B.3)

We deduce the extinction on each line of sight by using the measured value of the intrinsic color between J and K (I_int(J/K) = 1.), its definition (Eq. (B.3)) and the previous assumptions (Eqs. (B.1) and (B.2)), ${\mathit{I}}_{\mathrm{int}}\mathrm{(}\mathit{J}\mathit{/}\mathit{K}\mathrm{)}\mathrm{=}\mathrm{1.}\mathrm{=}\frac{{\mathit{I}}_{\mathit{K}}^{\mathrm{\ast }}\mathrm{\times }\exp \mathrm{\left(}{\mathit{\tau }}_{\mathit{K}}\mathrm{\right)}}{{\mathit{I}}_{\mathit{J}}^{\mathrm{\ast }}\mathrm{\times }\exp \mathrm{\left(}{\mathit{\tau }}_{\mathit{J}}\mathrm{\right)}}\mathrm{=}\frac{\mathit{DIRBE}\mathrm{2}}{\mathit{DIRBE}\mathrm{1}}\mathrm{\times }\exp \mathrm{(}{\mathit{\tau }}_{\mathit{K}}\mathrm{-}{\mathit{\tau }}_{\mathit{J}}\mathrm{)}\mathit{.}$(B.4)

Finally, we use Eq. (B.3) to yield the stellar contribution in each band, ${\mathit{I}}_{\mathrm{3.6}}^{\mathrm{\ast }}\mathrm{=}{\mathit{I}}_{\mathrm{int}}\mathrm{(}\mathit{J}\mathit{/}\mathrm{3.6}\mathrm{)}\mathrm{\times }{\mathit{I}}_{\mathit{J}}^{\mathrm{\ast }}\mathrm{\times }\exp \mathrm{\left(}{\mathit{\tau }}_{\mathit{J}}\mathrm{\right)}\mathit{/}\exp \mathrm{\left(}{\mathit{\tau }}_{\mathrm{3.6}}\mathrm{\right)}\mathit{.}$(B.5)

Taking into account our assumption (Eq. (B.2)) and the conversion coefficient deduced from the extinction curve of Rieke & Lebofsky (1985) with R_V = 3.1, we obtain $\begin{array}{ccc}& & {\mathit{I}}_{\mathrm{3.6}}^{\mathrm{\ast }}\mathrm{=}{\mathit{I}}_{\mathrm{int}}\mathrm{(}\mathit{J}\mathit{/}\mathrm{3.6}\mathrm{)}\mathrm{\times }\mathit{DIRBE}\mathrm{1}\mathrm{\times }\exp \mathrm{\left(}{\mathit{\tau }}_{\mathit{J}}\mathrm{\right)}\mathit{/}\exp \mathrm{\left(}\mathrm{0.5}{\mathit{\tau }}_{\mathit{K}}\mathrm{\right)}\\ & & {\mathit{I}}_{\mathrm{3.6}}^{\mathrm{\ast }}\mathrm{=}{\mathit{I}}_{\mathrm{int}}\mathrm{(}\mathit{J}\mathit{/}\mathrm{3.6}\mathrm{)}\mathrm{\times }\mathit{DIRBE}\mathrm{1}\mathrm{\times }\frac{\mathit{DIRBE}\mathrm{2}}{\mathit{DIRBE}\mathrm{1}}\mathrm{\times }\frac{\mathrm{1.}}{\exp \mathrm{\left(}\mathrm{0.5}\mathrm{\right)}}\\ & & {\mathit{I}}_{\mathrm{3.6}}^{\mathrm{\ast }}\mathrm{=}{\mathit{I}}_{\mathrm{int}}\mathrm{(}\mathit{J}\mathit{/}\mathit{K}\mathrm{)}\mathrm{\times }{\mathit{I}}_{\mathrm{int}}\mathrm{(}\mathit{K}\mathit{/}\mathrm{3.4}\mathrm{)}\mathrm{\times }\mathit{DIRBE}\mathrm{2}\mathrm{\times }\frac{\mathrm{1.}}{\exp \mathrm{\left(}\mathrm{0.5}\mathrm{\right)}}\mathit{,}\end{array}$with I_int(J/K) = 1. and I_int(3.4 /K) = 1.7 (Bernard et al. 1994).

We obtain the diffuse emission map, ${\mathit{I}}_{\mathrm{diff}}\mathrm{\left(}\mathrm{3.4}\mathrm{\right)}\mathrm{=}\mathit{DIRBE}\mathrm{3}\mathrm{-}{\mathit{I}}_{\mathrm{3.4}}^{\mathrm{\ast }}\mathit{.}$(B.9)

In the same way, using I_int(4.9 / 3.4) = 2.1, we compute the other diffuse emission map, ${\mathit{I}}_{\mathrm{diff}}\mathrm{\left(}\mathrm{4.9}\mathrm{\right)}\mathrm{=}\mathit{DIRBE}\mathrm{4}\mathrm{-}{\mathit{I}}_{\mathrm{4.9}}^{\mathrm{\ast }}\mathit{.}$(B.10)

Appendix C: Source sample

---

Acknowledgments

This research has made use of observations from Spitzer Space Telescope and data from the NASA/IPAC Infrared Science Archive, which are operated by the Jet Propulsion Laboratory (JPL) and the California Institute of Technology under contract with NASA. Based on observations obtained with WIRCam, a joint project of CFHT, Taiwan, Korea, Canada, France, and the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institute National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research makes a large use of the CDS (Strasbourg, France) services, especially Aladin (Bonnarel et al. 2000), Simbad and Vizier (Ochsenbein et al. 2000) and of the NASA Lambda data base at GSFC. We thank J.-Ph. Bernard, Nathalie Ysard and Vincent Guillet for fruitful discussions, Michiel Min for providing us his agglomerates of monomers in advance of publication and the DIM ACAV and “Région Ile de France” for financial support. M.J. and V.-M.P. acknowledge the support of Academy of Finland grant 250741. M.A., J.S. acknowledge support from the ANR (SEED ANR-11-CHEX-0007-01). We thank the referee for the careful reading of our paper and the valuable comments.

References

All Tables

All Figures