© ESO, 2015

1. Introduction

Dust grains are present in most cosmic objects, and they add a rich spectrum of investigation methods to astrophysical research (see, e.g., Draine 2003; Henning 2010). Through the absorption, emission, and scattering of radiation they allow us to trace the local density and temperature structure of objects (for molecular cloud cores see, e.g., Nielbock et al. 2012; Lippok et al. 2013; Launhardt et al. 2013), and to reveal the action of young stellar objects (YSOs), which are enshrouded at short wavelengths (e.g., Henning et al. 1990). In the interstellar medium (ISM), dense regions are predominantly cooled by dust radiation (Goldsmith 2001). Chemical reactions on their surface contribute to the chemical processes in the gas (Herbst et al. 2005) with the available surface depending on the grain size distribution and the fractal degree. Dust grains are sensitive to magnetic forces because of their charge and composition (Lazarian 2007). In the planet formation process, small dust grains are the seeds for forming larger bodies in the accretion disks around YSOs (Steinacker et al. 2013a).

Important questions therefore include the origin of the grains, their composition, and size distribution. Dust grains are known to be produced and destroyed in supernovae events (e.g., Wesson et al. 2010; Indebetouw et al. 2014) and in the atmospheres of stars in their late evolutionary phases (e.g., McDonald et al. 2009). Their distribution into the ISM and their complex processing including destruction and growth is the subject of ongoing research (Jones & Nuth 2011; Andersen et al. 2011; Zhukovska & Henning 2014). Grain growth processes have been discussed mostly in the framework of accretion disks and planet formation. Cold and dense prestellar molecular cloud cores may also exist long enough to allow the collisional growth of grains (e.g., Ossenkopf 1993; Ossenkopf & Henning 1994; Weidenschilling & Ruzmaikina 1994; Ormel et al. 2009, 2011). There is evidence that grains are ice-coated in the dense ISM (Whittet et al. 1983), supporting the growth process by efficient sticking. However, it remains to be investigated to what extent gas turbulence can provide the relative velocities that are needed to coagulate grains as assumed in these molecular cloud core models.

Concerning modeled grain size distributions, Mathis et al. (1977) were the first to propose a distribution for the diffuse ISM. They suggested a power-law distribution with a slope of −3.5 and with a size limit of 0.25 μm for non-graphite grains (abbreviated as “MRN distribution”) and emphasized that only weak constraints could be derived for the larger grains in the distribution. In the following we denote size distributions as “MRN-type” when they have a power-law index of −3.5. When we refer to grain sizes, we will give the grain radius value. Since then a large number of observations and models considering thermal emission and extinction of grains have suggested both varying size distributions and grains with radii beyond 0.25 μm in the denser parts of the ISM (Stepnik et al. 2003; Kiss et al. 2006; Ridderstad et al. 2006; Schnee et al. 2008; Chapman et al. 2009; Chapman & Mundy 2009; McClure 2009; Paradis et al. 2010; Veneziani et al. 2010; Ysard et al. 2013).

Scattering by dust grains is another process that is sensitive to the grain size. Lehtinen & Mattila (1996) modeled the 1.25, 1.6, and 2.2 μm band (J, H, and Kn) scattered light images of the Thumbnail nebula and concluded that the limiting grain size exceeds 0.25 μm. The strong decline of scattering efficiency of ISM grains beyond the Ks band caused people to assume that the cores become dark in scattered light in the mid-infrared (MIR). It would take even larger grains of micron-size scale to elevate scattered light intensities back to the surface brightness (SFB) measurable by MIR detectors. Such radiation was indeed found inside the cloud L183 and named coreshine (Steinacker et al. 2010). This term distinguished it from cloudshine, as Foster & Goodman (2006) named the near-infrared (NIR) scattered light from the core parts that are not shielded from NIR radiation. Beneficial for the detection of this radiation was the fact that extinction has a minimum in the wavelength range 3-5 μm, allowing interstellar radiation to penetrate deeper into the core.

Steinacker et al. (2010) present the most complex spatial structure modeling that has been performed so far for a molecular cloud core. Based on Spitzer/IRAC data at 3.6, 4.5, and 8 μm, the 3D radiative transfer (RT) modeling of L183 required the assumption of population of grains with sizes up to 1 μm. They argue that the excess of SFB in the 4.5 μm band could not be fit by emission of transiently-heated grains due to the absence of a strong feature in this band. Despite the strong central extinction (A_V above 150) and high Galactic latitude (b = 36.8°), the physical conditions in the L183 clump are not intrinsically different from those seen in nearby cores. Correspondingly, the investigation data of Spitzer PI and Legacy survey data of a sample of 110 cores led to the conclusion that about half of all cores show hints of coreshine (Pagani et al. 2010b). It was shown that scattered light from cores is observable in a variety of star formation stages, from prestellar cores to disks around forming stellar objects, as well as for a wide range of morphologies, i.e. single or binary systems and in filamentary structures (see also Stutz et al. 2009, for a small subsample of cores with 3.6 μm emission). In the same study, the RT analysis of model cores revealed that without background, the scattered light morphology does not vary strongly with the position in the Galaxy aside from weak enhancement towards the Galactic center (GC). For more massive cores with 10 M_⊙, the morphology showed a stronger crescent of enhanced emission due to the combined action of extinction and multiple scattering.

Since then further studies have investigated the occurrence and conditions of coreshine detection, and this work is also devoted to scattered MIR radiation in molecular clouds and the grain properties causing it. Pagani et al. (2012) find a drastically reduced occurrence of scattered light (3 detections plus 3 uncertain cases out of 24 objects) in the Gum/Vela region. It has been suggested that a connection to the action of a supernova remnant blast wave affects the grain size population in this region. In general the analysis of coreshine in regions with an enhanced interstellar radiation field (ISRF) like the Vela region is impeded owing to the confusion with emission of stochastically-heated grains like PAHs (polycyclic aromatic hydrocarbons).

Observations in the bands at 2.2, 3.6, and 8 μm of the core L260 were modeled by Andersen et al. (2013). Their simple core model was able to reproduce the observed SFB profiles that assume an MRN-type size distribution with grains up to about 1 μm in size and a locally enhanced incident radiation field (enhancement factor 1.7).

Steinacker et al. (2014b) have analyzed the coreshine seen from the low-density core L1506C and find that it requires the presence of grains with sizes exceeding the MRN distribution. Using grain growth models, they argued that L1506C must have passed through a period of higher density and stronger turbulence to grow the grains in situ.

The conditions of detecting scattered light from low-mass molecular cores at 3.6 μm have been re-investigated in Steinacker et al. (2014a) including background effects. The limits derived from the RT equation indicate that extinction by the core prohibits detection in bright parts of the Galactic plane and especially near the GC. They show that the scattering phase function favors the detection of scattered light above and below the GC and, to some extent, near the Galactic anti-center. It was discussed that an enhanced radiation field may give rise to a coreshine signal even in the presence of a strong background.

Lefèvre et al. (2014) have investigated the key parameters for reproducing the general trend of SFBes and intensity ratios of both coreshine and NIR observations. Based on a careful determination of the background field, they find that for a sample of 72 sources, starless cores show 3.6 μm/4.5 μm SFB ratios above 2, while cores with embedded sources can have lower values. To constrain the dust properties, they used a rich grid of size-averaged dust models based on an extrapolation of standard spherical grains able to fit the observations in the diffuse medium. Moreover, the effects of fluffiness, ices, and a handful of classical grain size distributions were also tested. As spatial density distribution, an inclined ellipsoid with a Plummer-like profile with two masses was assumed. Their density structure is well-suited to centrally-condensed cores like L1544, which show a central depression in their SFB pattern, while the vast majority of cores with coreshine have no depression. They find that normal interstellar radiation field conditions are sufficient to explain coreshine with suitable grain models at all wavelengths for starless cores. According to their multiwavelength approach, the standard interstellar grains are not able to reproduce observations and only a few grain types meet the criteria set by the data.

This paper investigates coreshine in a sample of ten cores. Coreshine has been detected for all cores except for Rho Oph 9 where the emission could also be explained by expected PAH emission. We use deep warm Spitzer IRAC 3.6 and 4.5 μm band observations (Cycle 9 Coreshine Follow-up: Program ID 90109, PI R. Paladini, Paladini et al., in prep.)¹. The observations were motivated by the fact that silicate grains of a size around 1 μm are expected to have their largest scattering-to-absorption efficiency in the wavelength range 3–5 μm and that the largest grains should have their strongest impact in the 4.5 μm band.

In Sect. 2 we briefly summarize the data and their acquisition and processing. The SFB ratio model is derived in Sect. 3, outlining the different RT transfer methods and optical depth effects. The modeling is described in Sect. 4, and the maximum grain sizes are analyzed for each core assuming an MRN-type size distribution. In Sect. 5, we discuss the results and consider the influence of local isotropic radiation fields, and summarize our findings in Sect. 6. The determination of the off-core SFB is described in Appendix A, and in Appendix B we describe the opacities of the dust model.

2. Data

The 3.6 and 4.5 μm IRAC data used for this work were taken during cycle 9 of the Spitzer warm mission. These observations targeted a set of ten Galactic cores partially selected from the “Hunting for Coreshine” program (Program ID 80053, PI Paladini). In earlier data, nine cores show a clear indication of coreshine emission at 3.6 μm. Extinction at 8 μm can be found for all cores in the cold Spitzer data except for PLCKECC G303.09-16.04 (we use the abbreviation ecc806 in this paper) and L1506C. It is a core from the Planck Early Release Cold Core Catalogue (ECC), which is part of the Planck Early Release Compact Source Catalog (ERCSC; Planck Collaboration VII 2011; Planck Collaboration XXIII 2011). We selected a varied sample of cores for this study, some with particularly strong coreshine (e.g., L260, ecc806, L1262, L1512), some also previously well-studied for coreshine emission (L260, L1506C) or otherwise (L1544, L1498). Some of these and other cores have particular features of interest: binary core L1262, L1517A having similar cores nearby, L1512 with a particularly simple shape, L1544 with a central coreshine depression, L1506C and L1498 with unusually low densities, L1439 near more evolved young stellar objects, and Rho Oph9 with an expected strong PAH contribution. The original survey did not reach a comparable signal-to-noise (S/N) at both 3.6 and 4.5 μm, with the consequence that coreshine emission at 4.5 μm was typically at the level or even below the instrumental noise. This situation was rectified in Cycle 9. These new observations were carried out in full array mode with 100 s integration. A detailed description of the source selection, observing strategy, and data processing will be provided in a forthcoming paper (Paladini et al., in prep.). In this work we have used the IRAC Level 2 mosaics (or pbcd, post basic calibration data), generated by the Spitzer Science Center (SSC) pipeline and downloaded from the Spitzer Heritage Archive (SHA). These images are affected by column-pulldown, an instrumental effect that consists in a depression of the intensity registered by pixels along a given array column, and which is caused by a bright source, i.e. a star or a cosmic ray. The IRAC images are also not zody-subtracted, meaning that the contribution from zodiacal light is not removed by individual frames or from the final mosaics at any stage of the automatic processing.

Basic core properties gathered from the literature are given in Table 1. The sources are ordered by decreasing coreshine peak SFB (see also Table 2). The estimates for the central column density have error bars (factors of a few) partly because the values are derived using a specific dust model. The ranges given in the Table are either specified in the cited work or based on assuming uncertainties by a factor of 2. Values without citation are based on simple estimates and just serve as an indication. The given masses and sizes (mean core diameter) are derived using different methods and also depend on the definitions of the core boundary. Moreover, it usually takes more than one parameter to characterize the complex core structure. In some cases we have estimated the size from the Spitzer map. To overcome the variation in the optical depth due to the uncertainties, we will rely in the modeling on a maximum optical depth at 3.6 μm (2.2 μm in the case of L260) that is given in the last column. Its choice and role in the modeling is described in Sect. 3.7. For modeling individual cores with more sophisticated density models and using multiwavelength 3D RT, it will be worth re-evaluating the mass estimates based on the original data and using a more sophisticated structure model with more than just a few size parameters.

The meaning of τ_max is explained in Sect. 3.7.

To give an impression of the overall SFB levels, the source geometry, sizes, and the influence of neighboring stellar sources, we show in Fig. 1 the 3.6 μm off-subtracted SFB images of all sources scaled to the same distance. For the example of L260, Appendix A describes how we have determined for all cores this “off-core surface brightness” that is subtracted from the core SFB. The linear scale in the plane-of-sky (PoSky) is indicated by the 10 kau bar in the lower lefthand corner. To indicate the approximate location of the central part of the core, we have marked the minimum of the corresponding 8 μm maps (for L1498, L1262, L260, L1439, L1544, L1517A, and L1512) where the cores are seen in extinction or the maximum of thermal emission maps in the FIR/mm (ecc806 and L1506C). Details are given in Sect. 4.

3. Modeling approach

Coreshine modeling involves RT calculations in a complex 3D core geometry. The scattering integral in the RT equation mixes radiation from all directions in all points. Moreover, to describe a 3D density structure, about ten to hundreds of free parameters are needed depending on the complexity of the structure. This makes coreshine modeling with an anisotropic illumination of the dust grains numerically difficult. A scan of the entire parameter space of size distributions, chemical compositions, grain shapes, and spatial density distributions is therefore too demanding for the current computer power. In the following sections, we describe our basic assumptions for reducing the parameter space and the complexity of the problem, the RT approximations, and the general modeling approach.

3.1. Basic assumptions

To reduce the complexity of the problem we make three basic assumptions.

Assumption 1: the grain size distribution s(a) da is normalized to have the same gas-to-dust mass ratio across the core, and it has no spatial change in shape or size limits.

This assumption is supported by the shallow gradients of the SFB pattern observed across the cores with coreshine. In turn, Steinacker et al. (2010) used a grain model with grain size increasing with gas density, and the SFB model showed stronger gradients than observed. Nevertheless, there are theoretical arguments for why the opacities might change across the core. For example, grain coagulation and the thickness of the ice mantles are expected to increase toward the center of the cores, depending, however, on how quickly turbulent mixing can smooth these gradients again (see also Andersen et al. 2013). Relaxation of this approximation will be the subject of further publications.

Assumption 2: the size distribution is a power law, and the sizes range from a₁ to a₂.

We have varied the power-law index and discuss its impact on the findings in Sect. 4.1. While there are studies suggesting that the spectral index of the distribution may be different from −3.5 as observed in the ISM (Weingartner & Draine 2001) or deviate from power laws (Weidenschilling & Ruzmaikina 1994; Ormel et al. 2009), we use this distribution as a typical case where the total grain surface area is dominated by the small grains while the mass predominantly rests in the large grains. It needs to be verified whether currently available data allow constraining the detailed size-distribution beyond the power-law picture.

Assumption 3: the optical properties of the grains are derived from Mie theory using spherical ice-coated silicate and carbonaceous grains.

Numerical calculations by (Ormel et al. 2009) that are based on aggregate collisions show that grains in cores go through a phase of compaction so that the spherical approximation might be better than initially expected from a growth picture that involves fluffy aggregates with large geometrical filling factors. Nevertheless, future modeling with a more complex dust model should include non-spherical grains (as discussed, e.g., in Lefèvre et al. 2014). Whittet et al. (1983) showed that grains in the Taurus cloud contain ice, and later studies have confirmed this finding for other clouds. Andersen et al. (2014) find for the Lupus IV molecular cloud complex that the limits for the occurrence of ice mantles and coreshine are similar and possibly related. We use a mass ratio of 1:4 for the carbonaceous and silicate component and assume the same size distribution. The details of the dust model are described in Appendix B.

3.2. Optical depth

The essential quantity for characterizing the impact of the dust on the transport of radiation is the size-integrated optical depth between two points on a line x = x₀ + ℓn with the starting point x₀ and the direction n. Based on the approximations listed in Sect. 3.1 it has the form ${⟨\mathit{\tau }⟩}_{\mathit{s}}\mathrm{\left(}\mathit{\lambda ,}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\underset{\mathrm{0}}{\overset{{\mathit{\ell }}_{\mathrm{1}}}{\mathrm{\int }}}\mathrm{d}\mathit{\ell }\mathit{t}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\underset{{\mathit{a}}_{\mathrm{1}}}{\overset{{\mathit{a}}_{\mathrm{2}}}{\mathrm{\int }}}\mathrm{d}\mathit{a}\mathit{s}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{\sigma }\mathrm{\left(}\mathit{\lambda ,a}\mathrm{\right)}$(1)where t is the spatial distribution of the H₂ number density, and λ the wavelength. Integrating t over space gives the total number of H₂ and, when multiplied with the H₂ mass, the assumed total core mass. Since ⟨τ⟩_s contains the line-of-sight (LoS) integral of the gas density, it can be related to the H₂ number column density N between the two points. The cross section σ refers to absorption, scattering, or extinction, which defines the optical depth ⟨τ_abs⟩_s, ⟨τ_sca⟩_s, and ⟨τ_ext⟩_s, respectively.

An important feature of the optical depth is that it depends not only on the core properties but also on the dust properties, so that for the same core, for example, two dust models can give very different central optical depths. This coupling of the core density structure and the dust properties with the optical depth is not easy to disentangle: basic properties like total core mass, its size, or visual extinction are often derived by interpreting observational data with a specific dust model and will change when the modeling requires a change in the dust properties.

In this paper, we use three approximations of the general 3D RT equation for this application based on the assumed maximum optical depth and give them simple names with the precise definitions provided in this section.

3.3. “Full RT”

For cores reaching optical depths for extinction beyond 1, shadowing effects and multiple scattering become important in the core. This is the most complex of the three modeling variants, and the stationary 3D RT equation in this case has the form

$\begin{array}{ccc}& & \frac{\mathrm{d}\mathit{I}\mathrm{\left(}{x}\mathit{,}{n}\mathit{,\lambda }\mathrm{\right)}}{\mathrm{d}\mathit{\ell }}\mathrm{=}\mathrm{-}{⟨{\mathit{\sigma }}_{\mathrm{ext}}⟩}_{\mathit{s}}\mathrm{\left(}\mathit{\lambda ,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathit{t}\mathrm{\left(}{x}\mathrm{\right)}\mathit{I}\mathrm{\left(}{x}\mathit{,}{n}\mathit{,\lambda }\mathrm{\right)}\mathrm{+}\mathit{S}\mathrm{\left(}{x}\mathit{,}{n}\mathit{,\lambda }\mathrm{\right)}\\ & & \mathrm{+}\mathit{t}\mathrm{\left(}{x}\mathrm{\right)}\underset{}{{\mathrm{\int }}^{\mathrm{4}\mathit{\pi }}}\mathrm{d}{\mathrm{\Omega }}^{\mathrm{\prime }}{⟨{\mathit{\sigma }}_{\mathrm{sca}}\mathit{p}⟩}_{\mathit{s}}\mathrm{\left(}\mathit{\lambda ,}{n}\mathit{,}{\mathrm{\Omega }}^{\mathrm{\prime }}\mathit{,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathit{I}\mathrm{\left(}{x}\mathit{,}{{n}}^{\mathrm{\prime }}\mathit{,\lambda }\mathrm{\right)}\mathit{.}\end{array}$

The equation describes the change in the intensity I given at the location x in the direction n along the path dℓ. The source term S contains the incident interstellar radiation field or nearby radiation sources. We have abbreviated the integrals over the grain sizes by ${⟨{\mathit{\sigma }}_{\mathrm{ext}}⟩}_{\mathit{s}}\mathrm{=}{\mathrm{\int }}_{{\mathit{a}}_{\mathrm{1}}}^{{\mathit{a}}_{\mathrm{2}}}\mathrm{d}\mathit{a}{\mathit{\sigma }}_{\mathrm{ext}}\mathrm{\left(}\mathit{\lambda ,a}\mathrm{\right)}\mathit{s}\mathrm{\left(}\mathit{a}\mathrm{\right)}$(4)and ${⟨{\mathit{\sigma }}_{\mathrm{sca}}\mathit{p}⟩}_{\mathit{s}}\mathrm{=}{\mathrm{\int }}_{{\mathit{a}}_{\mathrm{1}}}^{{\mathit{a}}_{\mathrm{2}}}\mathrm{d}\mathit{a}{\mathit{\sigma }}_{\mathrm{sca}}\mathrm{\left(}\mathit{\lambda ,a}\mathrm{\right)}\mathit{p}\mathrm{\left(}\mathit{\lambda ,a,}{\mathrm{\Omega }}^{\mathrm{\prime }}\mathit{,}{n}\mathrm{\right)}\mathit{s}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{.}$(5)The probability of scattering radiation from the solid angle dΩ′ into the considered direction n is given by the phase function p. This non-local integro-differential equation requires applying advanced numerical schemes like Monte Carlo or ray-tracing in order to solve for the intensity emerging in the direction of the observer (for a review see Steinacker et al. 2013b).

An essential difficulty in analyzing radiation received from a core with optical depths beyond 1 is that the structural information is mixed in at all LoS both by extinction effects and by multiple scattering. Since most cores show complex spatial structures, the error from assuming a wrong density model will be propagated with every ray or photon propagation into the derived SFB. It is important to note that this is not a problem of the available solvers. It is a problem of the modeling, which tries to avoid using too many free parameters that enter with the asymmetric core structure. Shadowing effects also amplify the impact of anisotropic illumination. The RT modeling of coreshine from the core L183 performed in Steinacker et al. (2010), for example, required the use of 100 3D Gaussian density clumps. An illustration of the RT variants is given in Fig. 2 with the full RT case being described in the righthand panel.

3.4. “Single-scattering RT”

When the optical depth for scattering is below 1 across the core, the impact of photons that scatter more than once is reduced. In this case we can use the single-scattering approximation. Instead of following rays or photons in arbitrary directions, we can separate contributions along the LoS (x-direction) in each PoSky position (y,z). Constructing the SFB from the background intensity I_bg that undergoes the extinction and radiation that is scattered into the LoS from all directions, we find (see, e.g., Steinacker et al. 2014a) $\begin{array}{ccc}& & \mathit{I}\mathrm{\left(}\mathit{\lambda ,y,z,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{ext}}\mathrm{\left(}\mathit{\lambda ,y,z,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}}{\mathit{I}}_{\mathrm{bg}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{+}{\mathit{I}}_{\mathrm{fg}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\\ & & \end{array}$with the foreground SFB I_fg, the direction-averaged product of incident radiation field, and phase function P now giving the probability of scattering radiation from the solid angle dΩ into the LoS

$\begin{array}{ccc}& & {⟨\mathit{IP}⟩}_{\mathrm{\Omega }}\mathrm{\left(}\mathit{a,x,y,z}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{ext;cell}\mathrm{-}\mathrm{obs}}\mathrm{\left(}\mathit{\lambda ,x,y,z}\mathrm{\right)}}\\ & & \mathrm{\times }\underset{}{{\mathrm{\int }}^{\mathrm{4}\mathit{\pi }}}\mathrm{d\Omega }\mathit{P}\mathrm{\left(}\mathit{\lambda ,}\mathrm{\Omega }\mathit{,a,x,y,z}\mathrm{\right)}{\mathit{I}}_{\mathrm{ISRF}}\mathrm{\left(}\mathit{\lambda ,}\mathrm{\Omega }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{ext;sky}\mathrm{-}\mathrm{cell}}\mathrm{\left(}\mathit{\lambda ,}\mathrm{\Omega }\mathit{,x,y,z}\mathrm{\right)}}\mathit{.}\end{array}$(8)

Interstellar radiation on its way to the core position (x,y,z) from a certain solid angle Ω undergoes extinction described by τ_{ext;sky−cell}. The radiation that is scattered by the dust at (x,y,z) and leaves towards the observer again undergoes extinction with the optical depth τ_{ext;cell−obs}.

The extinction of radiation still carries 3D density information from the entire core into the intensity of radiation being scattered into the LoS, and with it all errors of the density approximation. But the directions of the radiation transport are known and can be analyzed in advance.

The advantage of this approach is that it allows precalculation of all number column densities that are encountered by photons arriving from all directions at the considered LoS and then propagate to the observer. This speeds up the calculations when several wavelengths are considered and allows for the use of higher resolution ISRF maps, for example. We have created a code version that makes use of this advantage (used, e.g., in Andersen et al. 2013) and apply it also in this paper.

Owing to optical depth effects, the two bands at 3.6 and 4.5 μm can probe different regions in the outer parts of cores with masses of a few M_⊙ when the opacities change from 3.6 to 4.5 μm. For a standard core profile (flat followed by power law with an index of −1.5), a mass of 5 M_⊙, an outer radius of 10 kau, and an MRN-type distribution of silicate grains, among others, the locations of the τ = 1 layers differ by 14% of the core radius, and the probed column number densities differ by a factor of 1.8. For single-scattering RT like for full RT, the SFB in both bands is a result of the combined action of opacity and spatial density effects.

3.5. RT in the optically thin limit (“Thin RT”)

For some cores with low column densities and/or low-opacity grains, it is possible to treat RT in the single scattering and optically thin approximation even for grains with sizes of several μm in the wavelength range 3–5 μm. The limiting τ values for the thin approximation can be related, for example, to the considered order of the Taylor expansion of the exponential function (Steinacker et al. 2014a).

The SFB for an expansion to first order in τ_ext is

$\begin{array}{ccc}& & \mathit{I}\mathrm{\left(}\mathit{\lambda ,y,z,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{bg}}\mathrm{\left(}\mathit{\lambda ,y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{fg}}\mathrm{\left(}\mathit{\lambda ,y,z}\mathrm{\right)}\mathrm{\approx }\\ & & \mathrm{-}{⟨{\mathit{\sigma }}_{\mathrm{ext}}⟩}_{\mathit{s}}\mathrm{\left(}\mathit{\lambda ,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathit{N}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}{\mathit{I}}_{\mathrm{bg}}\mathrm{\left(}\mathit{\lambda ,y,z}\mathrm{\right)}\\ & & \end{array}$(9)

The lefthand side of Eq. (9) and the radiation behind the core I_bg can be determined observationally. Also for the incident field I_in, an approximation can be derived from observed allsky maps. Then we can calculate the righthand side with a model for the density and dust opacities. Unlike the cases of full or single scattering RT, the spatial and dust property problem decouple for thin RT since both terms on the righthand side of Eq. (9) contain the spatial variation in the form of N for each pixel or LoS. Since for optically thin scattering, both bands see the same column of grains, the spatial density variation is separated from the opacity variation. This makes the low-density cores best-suited to studying the optical properties of the grains as performed for translucent clouds and the diffuse ISM.

3.6. Modeling surface brightness ratios instead of surface brightnesses

For the modeling approach of the two images obtained at 3.6 and 4.5 μm for each source, we considered three methods partially based on the ratio, R, of the off-subtracted images.

i): Image modeling, where we assume a dust model and a densitystructure, calculate model images in the two bands, and comparemodel and real images at each wavelength. The modelingprocess is to vary the density structure and opacity until thedifferences are minimized.

ii): “Full RT” R and maximum surface brightness modeling, for which we build the ratio map from the observed images, select a region where the optical depth is expected to be small, and determine the observed R distribution. Then assuming a dust and density model, we calculate the model images in the two bands using “full RT”, build the ratio map, and find the theoretical R distribution in the selected region. From the images we find the maximum SFBs. The modeling process is to vary the opacity and the central column density until the differences are minimized.

iii): “Thin RT” R and maximum surface brightness modeling, by building a ratio map R from the observed images in the two bands, selecting a region where the optical depth is expected to be small, and determining the observed R distribution. Then assuming a dust model and a maximum column density, we can calculate the maximum SFB. The modeling process is to vary the opacity and the central column density value until we can reproduce the mean R and the maximum SFB.

To avoid either very complex spatial modeling with many parameters or the introduction of large errors from taking a density structure model that is too simple, we applied methods ii) and iii). To build the ratio map R = I_cs(3.6 μm) /I_cs(4.5 μm) from the SFBs above the background caused by scattered light in the core I_cs, we divided Eq. (9) for the two bands. For each pixel in the two images, the ratio on the lefthand side can be determined from the observed data ${\mathit{R}}_{\mathrm{obs}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{I}}^{\mathrm{3.6}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{bg}}^{\mathrm{3.6}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{fg}}^{\mathrm{3.6}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}}{{\mathit{I}}^{\mathrm{4.5}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{bg}}^{\mathrm{4.5}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{fg}}^{\mathrm{4.5}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}}\mathrm{·}$(10)Using the albedo ${\mathit{w}}^{\mathit{\lambda }}\mathrm{\left(}\mathit{a,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{sca}}^{\mathit{\lambda }}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{\mathrm{⟨}{{\mathit{\sigma }}_{\mathrm{ext}}^{\mathit{\lambda }}}^{\mathrm{⟩}}}_{\mathit{s}}\mathrm{\left(}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}}$(11)and the mean product of albedo and phase function ${\mathrm{⟨}{\mathit{w}}^{\mathit{\lambda }}{\mathit{p}}^{\mathit{\lambda }}\mathrm{⟩}}_{\mathit{s}}\mathrm{\left(}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathit{,}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}\underset{{\mathit{a}}_{\mathrm{1}}}{\overset{{\mathit{a}}_{\mathrm{2}}}{\mathrm{\int }}}\mathrm{d}\mathit{a}\mathit{s}\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{w}}^{\mathit{\lambda }}\mathrm{\left(}\mathit{a,}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}{\mathit{p}}^{\mathit{\lambda }}\mathrm{\left(}\mathrm{\Omega }\mathit{,a}\mathrm{\right)}\mathit{,}$(12)we can express the righthand side ratio for optically thin and single scattering as ${\mathit{R}}_{\mathrm{theo}}\mathrm{\left(}{\mathit{a}}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{⟨}{{\mathit{\sigma }}_{\mathrm{ext}}^{\mathrm{3.6}}}^{\mathrm{⟩}}}_{\mathit{s}}}{{\mathrm{⟨}{{\mathit{\sigma }}_{\mathrm{ext}}^{\mathrm{4.5}}}^{\mathrm{⟩}}}_{\mathit{s}}}\frac{\underset{}{{\mathrm{\int }}^{\mathrm{4}\mathit{\pi }}}\mathrm{d\Omega }{\mathrm{⟨}{\mathit{w}}^{\mathrm{3.6}}{\mathit{p}}^{\mathrm{3.6}}\mathrm{⟩}}_{\mathit{s}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}{\mathit{I}}_{\mathrm{in}}^{\mathrm{3.6}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{bg}}^{\mathrm{3.6}}}{\underset{}{{\mathrm{\int }}^{\mathrm{4}\mathit{\pi }}}\mathrm{d\Omega }{⟨{\mathit{w}}^{\mathrm{4.5}}{\mathit{p}}^{\mathrm{4.5}}⟩}_{\mathit{s}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}{\mathit{I}}_{\mathrm{in}}^{\mathrm{4.5}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{bg}}^{\mathrm{4.5}}}\mathrm{·}$(13)This equation contains only the dust properties (cross sections, phase functions, and size distributions) and the incident field (all-sky map I_in and background I_bg of the core).

Practically speaking, although some cores in the sample can be modeled with optically thin RT, we used single-scattering RT since the numerical effort is only marginally greater. The thin RT case is considered here because of its properties for separating dust opacity and structure, an advantage that is kept to some extent for higher optical depths.

To investigate the errors from using the three approximate RT schemes, we performed RT calculations for model cores with the masses typical of the range of low-mass cores 0.75, 1.7, and 2.4 M_⊙ and chose the core scale so that the optical depth for extinction through the center is 0.3, 0.7, and 1, respectively. Figure 3 compares thin and single-scattering RT in the upper panel. The images show the SFB error of the core at 3.6 μm | I_thin−I_single | /I_max in percent with the maximum SFB I_max of I_full. As expected, the maximum error for τ< 0.3 is small (about 2%). At τ< 0.7, the error is in the 5% range, for τ< 1 it reaches 12%. Comparing thin and full RT, the error up to τ< 0.7 is the same since the extinction effects and multiple scattering can be neglected. For τ< 1, both effects start to modify the resulting SFB, and the error rises to 15%.

3.7. The “limited optical depth” approach

To decide which core allows for which RT scheme, we need to specify the optical depth for the cores in the sample. However, neither the mass and density distribution in the core nor the optical properties of the dust grains are well known (the targets of investigation in this paper). The ranges in Table 1 give an indication of uncertainties estimated by authors and are often factors of a few.

Coreshine and scattered light from cores in general have a feature that can aid us in constraining the optical depth. When the central optical depth increases beyond one, the SFB profiles start to flatten and then show a central depression. The exact value when the depression is clearly visible depends on the properties of the core, the incoming field, and the background field, but in our experience, τ = 2 is a reasonable value to expect the occurrence of depressions in the images. Figure 1 contains two prominent cases (L1544 and Rho Oph 9) and a less pronounced one (L1262), but most of the cores have a pattern with a central maximum. In the modeling, we can use this observational evidence to consider only models with a central optical depth below a few when no central depression is visible at the considered wavelength. Andersen et al. (2013) used this argument to exclude grains with sizes beyond 1.5 μm since with the availability of NIR data, it is possible to find the wavelengths where depression sets in.

We therefore chose to apply a method which we call “limited optical depth” approach. It adjusts density and opacities in a way that they do not exceed a given τ limit. Fixing the density (as done, e.g., in Lefèvre et al. 2014) and then considering various dust models, τ may rise above a few and produce a central depression while the data show a maximum in SFB near the projected core center. Especially, models with steep density profiles and small kink radii (below which the profiles are assumed to be flat) accumulate most mass in the center leading to high central column densities. Furthermore, central column densities are usually derived from thermal emission measurements in the FIR/mm wavelength range based on a specific dust model, using the column densities with another dust model thus might be incorrect. The limited optical depth allows us to avoid this specification and the dust properties are consistently used only in the modeling of the coreshine.

In Table 1 we summarize the final maximum optical depth at 3.6 μm (2.2 μm in the case of L260) that we have used for the core modeling. Cores without central depression have values below 2. For the cores with τ_max> 2, we chose values that scale approximately with the average column density value given in the literature. The modeling is performed within the ranges of the mass and optical depth given in Table 1.

3.8. Properties of the theoretical surface brightness ratio R

In this section, we discuss how the cores can appear in the two bands, which basic physical properties are responsible for the appearance, and what the resulting range of theoretical SFB ratios is. Since the largest grains in the size distribution scatter the radiation more efficiently, our goal is to explore the maximum grain sizes when comparing the observed R distribution across the core with the theoretical R of a particular dust model.

Besides scattering the incident radiation toward the observer, the grains are also responsible for the extinction of background radiation along the LoS. This background extinction introduces several complications. The scattered light has to overcome the extinction of the background radiation by the core to arrive at a SFB above the background and to appear in excess of emission. The four possibilities are to see the core in extinction in both bands, in just one of them, or in emission in both bands. Since for our core sample we have analyzed only parts of the cores with excess of emission in both bands, we concentrate on the last case.

To keep the example simple, we consider the condition for excess of emission in the optically thin case. From (13), we find ${\mathit{I}}_{\mathrm{bg}}^{\mathit{\lambda }}\mathit{<}\underset{}{{\mathrm{\int }}^{\mathrm{4}\mathit{\pi }}}\mathrm{d\Omega }{\mathrm{⟨}{\mathit{w}}^{\mathit{\lambda }}{\mathit{p}}^{\mathit{\lambda }}\mathrm{⟩}}_{\mathit{s}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}{\mathit{I}}_{\mathrm{in}}^{\mathit{\lambda }}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}$(14)for both bands. The most prominent effect on R is the pole causing a singularity when the 4.5 μm denominator vanishes. Physically, this happens when the scattered light SFB is equal to the extincted background SFB. Increasing the maximum grain size, the SFB in the two bands increases due to the higher scattering efficiency of larger grains. In the picture of Rayleigh scattering, the condition 2πa ≪ λ is less fulfilled for larger grains. The 4.5 μm band SFB benefits stronger from the increased maximum size due to this condition and R decreases.

Figure 4 shows in panel A the numerator and denominator of R for the example of L260. We have used ice-coated silicate grains and an MRN-type size distribution ranging from 0.01 μm to a₂. We follow the approach described in Andersen et al. (2013) to calculate the incoming fields I_in from zodiacal-subtracted DIRBE² maps. The calculation is done in single-scattering approximation.

The left yellow regions indicate where the core is seen in extinction in both bands or just the 4.5 μm band. The observed coreshine SFB defines the limit of this region when becoming zero. The thick dashed lines give the maximum observed coreshine SFB, and we also show the uncertainty in the value due to the choice of background (see Appendix A) by thin dashed lines forming the light blue and green error ranges. To indicate the resulting range of maximum grain sizes, we label the crossing points with the observed values and their error range derived from the off-measurement uncertainty by circles. Bars near the x-axis aid in reading the a₂-range corresponding to the agreement of observed and theoretical coreshine. As is visible in the figure, the bars partially overlap, indicating that by extending the distribution up to sizes around 0.84 μm, the model reproduces the observed data.

Andersen et al. (2013) have modeled the 2.2, 3.6, and 8 μm SFB along a cut through L260, and found maximum grain sizes around 1 μm, which is close to what we find for our ice-coated silicate and carbonaceous grains. Their calculations required considering an enhanced incident DIRBE field boosted by a factor of 1.7, which might have been due to an error in the older RT code and is no longer needed in the new code version.

In panel B, we show the ratio of the blue and green curves from panel A, which is the theoretical band ratio R with the same notations, and the bullets indicate the crossing points of the model with the observed maximum coreshine SFBs. We show just the bullets and not the dashed lines from panel A to be able to plot more than one R curve in the following figures without confusing the plots. As expected, the curves show the increase due to the pole in R when the 4.5 μm SFB cannot overcome the background extinction, and the core turns to an extinction appearance at smaller maximum sizes.

In the optically thin case, the curve in panel B should be independent of assumptions about the spatial density distribution within the core, while the SFB profiles in panel A assume a central column density or optical depth and therefore rely on spatial parameters. For optical depths higher than unity both panels will be affected by the assumptions made within the density distribution model.

We note that R also depends on the position of the core since the ISRF is anisotropic with a maximum at the GC, and convolving it with the anisotropic phase function will lead to different values for different core positions. As a technical note we add that when referring to the theoretical quantities like the optical depth, we skip the <>_s notation, indicating the integral over grain size to improve readability from this point on.

4. Coreshine modeling of the core sample

The general approach taken in this section is to compare the measured R distribution in a region of the core with the theoretical ratio R_theo derived from a model core. The core properties and the maximum grain size are varied without exceeding the limiting optical depth (Sect. 3.7). Since two sets of 3.6 and 4.5 μm images can have the same ratio map but might differ strongly at each wavelength, we also verify the agreement of the SFBs. For cores with τ< 2, it is sufficient to consider a single SFB value. We have chosen the maximum observed SFB and not some average for two reasons. First averaging considers more regions with worse S/N, increasing the overall error, and second the averaging would bring in more assumptions about the way the averaging is done. Even if the core has regions with τ> 2, we try to select a region in the observed and model image that has optical depth below that limit. In this way we take advantage of the optically thin approximation being valid, which removes spatial dependences from the problem. For cores with τ> 2, we perform “full RT” modeling by calculating images, building ratio maps and an R distribution, and fitting the observed R distribution and the central column density within the optical depth limit.

Throughout the section, we make use of the approximate core properties summarized in Table 1.

4.1. Modeling details explained: the core L260

To describe the modeling approach in detail, we have chosen the core with the brightest coreshine in the sample: L260. According to Steinacker et al. (2014a), the core is located in a region favorable to coreshine detection because it is near the GC in the PoSky, but the latitude is high enough to prevent a strong background. Still it has the second strongest 3.6 μm background of the cores modeled in this paper.

Coreshine is a weak extended feature, and especially the 4.5 μm SFB is only factors 2–8 larger than the sum of instrumental and point source noises. Correspondingly, we have selected special regions in the images to minimize the impact of noise and stellar point spread function (PSF) on the analysis. The two background SFBs ${\mathit{I}}_{\mathrm{bg}}^{\mathrm{3.6}}$ and ${\mathit{I}}_{\mathrm{bg}}^{\mathrm{4.5}}$ entering R_theo are taken from the DIRBE maps by subtracting the stellar flux measured by WISE³. We did not subtract a diffuse foreground component of the radiation field used as background of the core since the current models give only indications (for a discussion see Lefèvre et al. 2014). We briefly address its effect in Sect. 5. For 3.6 μm, the bands of IRAC and DIRBE are comparable, for 4.5 μm we interpolate from the other bands (see also Andersen et al. 2013). The sum of the background and foreground contributions can be determined from the SFB near but off the core (“off” characterizes quantities derived in this way) ${\mathit{I}}_{\mathrm{off}}^{\mathit{\lambda }}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{bg}}^{\mathit{\lambda }}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{+}{\mathit{I}}_{\mathrm{fg}}^{\mathit{\lambda }}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathit{.}$(15)The details of how this subtraction is carried out are explained in Appendix A. We note that the subtraction removes the possible zodiacal light contribution.

Figure 5 shows the steps in the analysis for the core L260. The off-subtracted SFB maps at 3.6 and 4.5 μm in a 5.5′ × 5.5′ region around the core are shown in the panels A and B, respectively. The cold Spitzer 8 μm map is given in panel C to indicate the location and the shape of the core (“Coreplanets_cores”: Program ID 139, PI N. Evans). L260 is seen in emission in 3.6 and 4.5 μm. The white square indicates the region where SFB ratios are measured. Since the core shows no sign of central depression in the two bands, the square is placed in the central region of extinction in the 8 μm map. The size and precise position of the square have been adopted to avoid the impact from nearby stars and their PSF, to minimize the impact from column pull-down in the image, and to recover the grain information in the part of the core where the coreshine S/N is maximum in both bands. Alternative shapes of the measuring area would have been possible, but an investigation of the ratio maps (panel G, see farther down) showed that the results are not sensitive to the shape of the area. All remaining pixels that show impact from a stellar PSF or column pull-down are masked and not modeled.

To check the distribution of optical depths for extinction in the region where we measure R, we show the direction–averaged optical depth map of a core with the overall properties of L260 in panels D and E. The direction average takes into account that radiation is extincted while reaching the LoS from outside and on its way to the observer. As simple core model, we use an ellipsoidal core with two of the three axes identical and a radial density profile that flattens in the inner part while following a power law with index −1.5 outside. Core mass and outer radius are taken from Table 1, and the kink radius is assumed to be one-third of the outer radius. We note that both the density power-law index and the kink radius will depend on the evolutionary stage of the core, but more evolved cores might have a smaller kink radius, for example, and for L1544 we have used a different radial profile. For L260, the mean optical depth values reach about 0.4 for 3.6 μm. We there fore use single-scattering RT to model the core. We note that the models shown in panels D and E for the core sample do not aim to model the spatial density distributions but only serve to estimate the optical depth variation, with the exception of the cores with maximum optical depth above 2 where we create images in full RT based on the shown spatial model.

To read approximate values of the SFB band ratio, the noise level and the contamination by stars, we show in panel F horizontal (thick) and vertical (thin) cuts through the center of the white line frame for 3.6 (4.5) μm with black (green) lines. While there are local gradients along the cuts that vary the surface brightness ratio, the figure already allows reading an approximate ratio of about 2 to 4. The cuts also indicate the level of stellar radiation impact and the overall noise of of about 0.005 MJy/sr. Figure 4 shows the two maximum surface brightness values and their error bars as dashed vertical lines. The error bar contains not only the error read from panels A, B, and F in Fig. 5, but also the error in the off-measured surface brightness that is subtracted from the measured surface brightness. The entire map of off-subtracted surface brightness ratios R is shown in panel G. Most variations in the white square are due to stellar sources. We masked all pixels affected by stellar PSFs and by column pull-down. In panel H, we show the R pixel number distribution across all unmasked pixels in the white square. It has a maximum at R = 2.37, and as a measure of its size, we use the approximate width of the distribution at the half maximum value indicated by thin vertical lines (the full width at half maximum, FWHM, is about 1).

In panel I, the range of observational ratios R_obs is compared to the theoretical ratio R_theo as a function of a₂ as derived in Sect. 3, with the background values given in Table 2 and using the core sky position for the phase function integration of the scattering integral. The limited optical depth approach is applied, which prevents the maximum optical depth for extinction to exceed 1 at λ = 2.2μm (see Table 1). The green color-coding indicates the observed R-distribution from panel H.

The crossing points of the theoretical R curve with the observed coreshine range at 3.6 μm (blue) and 4.5 μm (green) are shown as dots as in Fig. 4 with small dots indicating the error range and large dots giving the crossing point with the exact value. There is an overlap of the a₂ ranges of the two bands indicated as blue and green horizontal bars at the bottom of the figure and the black bar marking the a₂-range where R_obs and R_theo agree. Thus the model is able to meet all three observational constraints for maximum grain sizes between 0.78 and 1.02 μm.

We have performed extensive tests on varying the spectral slope of the size distribution as suggested by Weingartner & Draine (2001). The lower grain size limit has no effect on the coreshine since the scattering is dominated by the largest grains as shown in Steinacker et al. (2010). The effect of the size distribution slope on the derived grain size limits is as expected: flattening the distribution will decrease the limits, as the number of bigger grains increases with a flatter slope. Since they are responsible for most of the observed coreshine, less big grains are needed to meet the observed flux. The impact on the limits is moderate since all opacities are integrated over the size distribution. For a slope of −2 instead of −3.5, the limits move by about 0.15 μm. More importantly, the overlap between the a₂-ranges derived from the two observational constraints remains the same. As a result, we argue that the derived grain size limits may vary by several tens of percent when changing the slope of the size distribution power law in a reasonable range of [−2 to −5], but this will not affect the conclusions. For the models presented in this section, we have used a spectral slope of −3.5.

4.2. ecc806 (PLCKECC G303.09-16.04)

The second brightest core in the sample, PLCKECC G303.09-16.04 (we will use the shorter name ecc806 here), is surrounded by more filaments with coreshine that are located outside the area shown in Fig. 6. As for L260, the position favors the detection of coreshine because it is close to the GC in longitude and far below (b = −16.04°) the Galactic plane. Moreover, it has the lowest background in the sample in both bands. There is no sign of depression in the images and in the cut in panel F. We therefore use single-scattering RT. The chosen white frame encloses the core center with optical depths below the τ limit. In panel C, the 250 μm map observed by Herschel (observation id 1342213209) as part of the Gould Belt Survey (PI: André) is shown for comparison. As with the extinction pattern at 8 μm shown for the other sources, the thermal emission pattern also matches nicely the coreshine pattern. The SFB profiles are remarkably flat in the central part. The R distribution is only slightly broader than that of L260 (FWHM about 1.2) but with a maximum at larger R near 3. The fit indicates that the maximum grain size in ecc806 is around from 0.4 to 0.52 μm.

4.3. L1262

Modeling of the starless core in the binary core L1262 (CB 244) is complex owing to the nearby YSO with its outflows creating H2 emission in the shocked gas in both bands. L1262 also has a remarkably extended cloudshine that is visible at shorter wavelength and is the most distant core in the sample (about 200 pc). The YSO and its warm-dust surroundings may give rise to an additional radiation field. Since the density distribution is complex, it is hard to model this component in order to derive the intensity and directional dependency at the location of the starless core. With 60° away from the GC, the phase function of the largest grains will not increase coreshine by beaming the strong GC radiationforward .

The core position in the PoSky is elevated enough (b = 12.4°) to have a low background, especially at 4.5 μm, which allows the detection of the weak coreshine signal in this band. The mass of the core is rather uncertain with a mean value of 5 M_⊙ and an error bar of 2 M_⊙. However, the column density range given in Table 1 also points to τ> 2, and the 3.6 μm image in panel A reveals a depression right where the extinction against the background has a maximum in the cold 8 μm image in panel C (“darkclouds_IRAC”: Program ID 94, PI C. Lawrence). We therefore have chosen maximum optical depth values around 4.5 (scaled with the mean column density values of the other cores) and use full RT. We have chosen a white frame away from the YSO and its extended tail and excluded the part of the square where the extinction seen in panel C of Fig. 7 indicates the center of the core. As mentioned before, this improves the modeling because at lower optical depths, some of the spatial dependencies cancel out. The observed R distribution is about a factor of 2 narrower than for L260 with a lower mean value of 1.75. Aside from the color-coded observational R-distribution, panel I shows the distribution of the theoretical R-values with a thick line for the mean value and thin lines for the FWHM. The a₂ ranges overlap between 0.52 and 0.71 μm.

4.4. L1517A

The L1517 cloud in Taurus contains four filaments with five embedded cores (Hacar & Tafalla 2011). For the modeling we have chosen core A since it was visible in both bands in our new deep observation images. The coreshine of all Taurus cores (in our sample L1517A and L1512, L1544, L1506C, and L1498 with l around 180° discussed later) benefit from a side maximum in the dust- model phase function of large grains that predicts enhanced backward scattering of the strong GC radiation (see Steinacker et al. 2014a).

The core is located at the border of the 3.6 μm image (Fig. 8, see also panel C, “darkclouds_IRAC”: Program ID 94, PI C. Lawrence), but were able to select a region that includes the core center and most of the core. We considered using cold 3.6 μm data instead, but decided to stay with the warm data for the comparison to minimize the impact of systematic effects between the observing runs because the data show no strong increase in noise compared to the warm 4.5 μm data. There is no indication of any depression, and the mass and column densities point toward a maximum optical depth value for extinction below 1. We therefore use single-scattering RT. The cuts show a clear constant factor between the two bands in both cut directions, as one would expect from optically thin scattered light where the SFB is proportional to the column density. This should give rise to a narrow R distribution. On the contrary, panel H shows a broader R distribution than L260 (FWHM is about 1.5, peak at 2.75). We attribute this to the source having the strongest background in the sample. The distribution also is asymmetric with a strong tail toward higher R. Panel G shows an increase in R values in the upper part of the white frame Probably because of the loss of the 4.5 μm coreshine signal, and we attribute the asymmetry in panel H to this signal loss. The model is able to meet all observational constraints for a₂ values above 1.1 μm.

4.5. L1512

This Taurus core is located close to the Galactic plane and has an enhanced background (second strongest 3.6 μm background in the sample). As for L1517A, we use the new warm data showing the core at the image border instead of using center-image 3.6 μm cold Spitzer data in Fig. 9. The cold Spitzer map at 8 μm shown in panel C is from the “darkclouds_IRAC” program (ID 94, PI C. Lawrence). Based on the core property estimates from Table 1, the optical depth analysis indicates that the central part of an ellipsoidal representative of L1512 would see optical depth above 1, and therefore extinction and second scattering make full RT necessary. The image in panel A and cuts like the one given in panel F show little central depression though, but strong stellar PSFs make this analysis difficult. We therefore performed full RT modeling. The 4.5 μm signal is just about a factor 2 above the noise with the image showing just hints of excess of emission. The R distributions are correspondingly broad (FWHM of about 2.6 around R = 3.6) and conclusions from the modeling must be made with caution. The full RT distribution for R overlaps with the observed R values for maximum grains between 0.48 and 0.61 μm.

4.6. L1544

The Taurus core L1544 is one of the few known cores across the sky that show a prominent inner depression of SFB in both bands. While the mass seems to be comparable to that of L1512 showing no clear depression, the central SFB of L1544 very likely drops because of optical depth effects of the centrally condensed radial profile. The 8 μm map shows a strong extinction pattern at the core center (“cores2deeper”: Program ID 20386, PI P. Myers). This would mean that we can expect moderate optical depth in the outer parts of the elliptically shaped core and a strong rise in τ near the center. Using a profile with a power-law index of −2.5 in the outer parts and a flattening starting at the kink radius of 1 kau, the mean optical depth maps shown in panels D and E in Fig. 10 provide such a pattern. Scaling τ with the estimated central column density of L1544, we use a τ limit of 6 at λ = 3.6μm. The white square region was placed outside the depression, and we use full RT. The derived ratios are low with a mean value of 1.25, and the FWHM is comparable to that of L260 (about 0.95). But as is visible in panel F, the data are extremely noisy with a dense stellar background field interfering and with a loss of the 4.5 μm signal in the lower white frame part that causes a wing in the R distribution at high R. The model suggests large grains beyond a₂> 1.5 μm to reproduce the SFBs and ratios.

4.7. L1506C

This extended Taurus core with low density, low turbulence, but strong depletion (Pagani et al. 2010a) was discussed in Steinacker et al. (2014b) as a core that either is still contracting to become a prestellar core (as suggested by Pagani et al. 2010a) or has passed through a core phase and is dissolving again. Since the core is invisible at 8 μm, we show in Fig. 11 instead a MAMBO II map (reproduced from Pagani et al. 2010a) at 1.2 mm revealing the same pattern (in thermal emission) as the coreshine. The large Galactic latitude of −17.57° helps to make the coreshine visible. The 4.5 μm signal is just a factor 2 over the noise but the entire core pattern is nicely visible and exhibits no strong impact of potential noise confusion. From the noisy cuts displayed in panel F, it is difficult to judge on potential central depressions, since the images in panels A and B do not reveal a darkening, and the basic parameters suggest optical depths below 1 in the center. As a result, we placed the white frame near the core center, thereby avoiding the two large stellar PSFs, and use single-scattering RT. The R map in panel G shows the interesting feature that the entire clump visible in the image seems to have R values comparable to the mean value measured in the white frame as 1.7 (FWHM about 1).

The model is able to account for the SFBes based on grains with maximum sizes around 0.4 μm, while the ratios produced by such grains of about 7 would be too high to match the observed low values between R = 1.2 and 2.2. This steep ratio points toward grains with sizes larger than 0.7 μm. The analysis of L1506C presented in Steinacker et al. (2014b) using a grain growth model and simple spatial mixing of the size distribution pointed toward large grains beyond 1 μm in size. However, including the spatial mixing in the growth process would likely reduce the efficiency of the growth process, leading to smaller maximal sizes.

4.8. L1439

The core L1439 (CB26) is located only 6° above the Galactic plane in the PoSky and is therefore seen against a crowded star field. The 4.5 μm signal is only a factor 2 above noise (Fig. 12). At its rim, a YSO creates a jet, and it is not clear whether the stellar component or the warm dust adds local radiation to the ISRF. The cold Spitzer map at 8 μm is shown in panel C (“darkclouds_IRAC”: Program ID 94, PI C. Lawrence). The core parameters listed in Table 1 make central optical depths below 1 likely. We therefore placed the white frame right on the core and used single-scattering RT. The modeling with a variable τ limit led to values around 0.5 for 3.6 μm. The core is located at the border of the 3.6 μm frame so that the dark righthand part of the R map displayed in panel G is meaningless. The R pixel number distribution shown in panel H of Fig. 12 contains an asymmetric wing from the map parts with a weak 4.5 μm signal. It has a maximum around 1.7 and a FWHM of 1.4. Panel I reveals that the model produces no common size range to meet the observed SFBs and the observed low R_obs ratios.

4.9. L1498

L1498 in Taurus is discussed to be on the verge of becoming a core (Lada et al. 2004), and only a faint extinction pattern is visible at 8 μm in panel C of Fig. 13 (“darkclouds_IRAC”: Program ID 94, PI C. Lawrence). Positioning of the white frame was guided by the 3.6 μm image border so that it was placed in one half of the ellipsoidal shape of the core. The R maps in panel G show that this choice does not largely affect the results since the other parts of the core seem to have similar values to those of the core L1506C. The core parameters suggest central optical depths in the range of L1506C, so we use the same optical depth limit and single-scattering RT. Given the low 3.6 μm signal, L1498 shows a surprisingly strong 4.5 μm SFB yielding to low ratios around 2 (FWHM is 1.55). The model grains produce too much coreshine when using the grains with size above 0.77 μm as suggested by the model to fit the R_obs distribution so that the a₂ ranges do not overlap.

4.10. Rho Oph 9

The core Rho Oph 9 is located in the most complex environment of all cores in the sample (Gal. coord. 354.37 +16.17). The cold Spitzer data shown in Pagani et al. (2010b) show nearby PAH regions in all four bands with reduced SFB at 4.5 μm owing to the lack of a main PAH feature in this band.

The excess of emission that is seen near the strong extinction pattern is also present at 4.5 μm similar to the pattern seen in L1544. Unlike the PAH emission, it disappears for larger wavelength in the cold Spitzer data as expected from scattered light. Figure 14 shows the 3.6 and 4.5 μm warm Spitzer images. The overall pattern is similar to L1544 with a strong central SFB depression and an ellipsoidal emission region around. But disentangling the coreshine emission from variations in the background and PAH emission turned out to be difficult. We have not modeled the core since we could not determine an off-region that was clearly free of PAH emission and near the core, and hope to return to the source when more data are available. We show the maps and the ratio map for completeness.

The peak and background SFBes of the modeled sources are summarized in Table 2 along with the derived observed SFB ratios R_obs.

5. Discussion

In Fig. 15 we summarize the a₂ ranges derived from the modeling with the same notation as in Fig. 5: blue, green, and black bars for agreement with the 3.6 μm maximum SFB, with the 4.5 μm maximum SFB, and with the observed SFB ratio R_obs, respectively.

For L260, ecc806, L1262, L1517A, L1512, and L1544, the bars overlap so that the model is able to reproduce the observations with a single grain size distribution. The derived maximum grain sizes for these sources are around 0.9, 0.5, 0.65, 1.5, 0.6, and >1.5 μm, respectively. The value derived for L260 of about 0.9 μm is consistent with the results derived in Andersen et al. (2013).

The grain size limit of about 0.5 μm in ecc806 is surprisingly low given that it is bright in coreshine (about 0.075 MJy/sr) and that large grains dominate the scattering. At a longitude of about 303°, radiation from the GC already needs large scattering angles to reach us. Because grains larger than 0.5 μm scatter mainly in the forward direction, their contribution to the observed signal is reduced. However, the low background due to its high elevation (b about −16°) aids for seeing a strong coreshine.

The structure around the starless core in the binary core system L1262 is complex and not represented well by an ellipsoid. Correspondingly, the derived images carry the errors from the spatial density because of both shading and multiple scattering. Nevertheless, maximum grain sizes of 0.65 μm agree with the measured ratios and maximum SFBs in the chosen region outside the core center. They may actually be smaller since the nearby YSO could contribute to the local radiation field both by direct irradiation and by causing hot dust and PAH emission in its vicinity. Further exploration with a better spatial model may answer that question.

Strong coreshine (about 0.055 MJy/sr) also comes from L1517A but on top of a much stronger background. Correspondingly the grain size limit a₂> 1.5μm derived by the model is larger than, for example, that of ecc806 with similar maximum SFB values to enable a stronger scattering. Like all cores in Taurus, L1517A benefits from side maxima of the phase function supporting back scattering of the strong GC signal.

For the cores with the lowest coreshine SFB L1506C, L1439, and L1498, the three different a₂ ranges do not overlap in one coherent way. The maximum grain sizes to explain the ratios are in all cases larger with just a lower limit. To investigate whether this systematic trend arises from a simplification of the model, we consider a local radiation field in addition to the ISRF. Since the cores are located in star-forming regions where warm dust is present near the YSOs, such a local radiation field is not unlikely. Moreover, the DIRBE map reveals radiation on top of the stellar component from the regions where cores are located, with a spectrum increasing with wavelength in the considered bands. While the radiation sources near the core will also contribute to the DIRBE flux from that region, if it is not shielded along the LoS, sources very close to the core will make an increased contribution to the local field due to the squared distance dependency of the flux.

The directional dependence of the local field is unknown and is likely to be different from that of the ISRF. To explore the impact of a local field, we assume a simple field with no directional variation and an intensity of 0.05 MJy/sr (about 1/6 of the mean DIRBE field) at 3.6 μm. We modify the spectral shape by a factor (λ/ 2.2 μm)^δ so that a variation in δ takes the rise in the local field intensity with wavelength into account. Figure 16 shows three cases for the core L1439. The thin line gives the theoretical R without a local field. When adding a constant field without spectral variation (δ = 0, medium thick line), more radiation is scattered toward the observer. The grains can thus be smaller to explain the observed SFB, but the overlap between the ranges for each band, and the ratio is not improved. Assuming a local field with an increase described by δ = 1.5, this trend continues to lead to no improvement in the fit. For cores with an YSO in the vicinity or even embedded, a local anisotropic radiation field is added to the incident field, and its impact on the coreshine ratios was discussed in Lefèvre et al. (2014). They find higher coreshine SFBes and lower observed coreshine SFB ratios I_3.6/I_4.5 for embedded cores than for starless cores. This would make it more difficult in our model to converge to a common maximum grain size. Alternatively, the size distributions explored in that paper and in Steinacker et al. (2014b) mark the beginning of a more systematic search of the large size distribution parameter space for a dust model that can explain the presence or absence of coreshine in the various observational bands.

For L260, Andersen et al. (2013) also used the observed Ks band to constrain the size distribution. At shorter wavelengths, the central depression is more likely. It needs to be verified that peaks in the grain size distribution at sizes around 1 μm will be compatible with the coreshine and extinction observed at shorter wavelengths. This emphasizes the importance of deep Ks band and NIR observations in general of cores with coreshine. Moreover, Lefèvre et al. (2014) argue that 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 about 2 μm in meaningful quantity.

A word of caution is needed for the uncertainties in the background estimates. A better estimate of the foreground diffuse radiation field would improve the current uncertainty that enters coreshine modeling. Weaker background would increase the coreshine SFB in both bands so that smaller grains would be sufficient. The cores where the model fails would require the opposite. Moreover, the background values derived by Lefèvre et al. (2014) using a different method can vary from one region to the other by a factor of up to two. This can modify the results for the cores with weak coreshine or strong background. For future 3D modeling of a single source, it could be important to put more effort into modeling the foreground and background contributions.

One particular finding about L1506C and L1498 is interesting. Both cores show little variation in R_obs across the entire core as far as 4.5 μm coreshine is detectable. For both cores this ratio is small, as expected from the presence of larger grains. Both are discussed as being on the verge of becoming low-mass cores with typical properties, and also show central gas depletion that requires a density that might not be reached in these cores.

For both cores the existence of larger grains would be surprising if the grains are to be formed by coagulation since their densities are not high enough to coagulate them, as indicated by studies of Steinacker et al. (2014b) for L1506C. One of the hypotheses presented in that work was that L1506C went through a phase of higher density and stronger turbulence to explain the existence of larger grains and strong gas depletion. With the similarities between the two cores also emerging from the work presented here, this could also be the case for L1498.

6. Summary

This work is based on deep 3.6 and 4.5 μm IRAC Spitzer warm mission data taken during Cycle 9. The sample is composed of ten cores, nine of which show coreshine at both 3.6 and 4.5 μm, and one where coreshine is unconfirmed owing to strong PAH emission. The observations and the work in this paper were motivated by the expectation that coreshine observed at MIR wavelengths might be dominated by dust grains in the micron size range and that opacity models predict a stronger impact of the largest grains at longer wavelengths.

Given the vast parameter space of the general grain scattering problem in cores caused mainly by the unknown density structure, the presented model approach makes use of the fact that in the optically thin limit, both the observed coreshine intensity and the intensity losses by extinction through the core are proportional to the column density of grains. By considering ratio maps of 3.6 and 4.5 μm intensities, the impact of the spatial density distribution is expected to be reduced for the cores with moderate optical depth <1. Moreover, we use the constraints from the observed SFBes based on the expected column densities. Since the core properties are known only within factors of a few, and the opacity changes caused by increasing the grain size can be substantial, it is essential for the modeling not to exceed the observed optical depth limits, which are provided in the form of depression of the central SFB for cores with optical depth for extinction above about 2. When the column density is varied within its uncertainties, we make sure that this optical depth limit is never exceeded.

We applied a model based on the opacities (cross sections and phase functions) of ice-coated silicate and carbonaceous spheres and a size-distribution following a power law (index −3.5) ranging from 0.01 to a variable maximum grain size for both species. The incident field is approximated by DIRBE maps, and the radiation right behind the core is constructed from the absolutely calibrated DIRBE map with point source contribution subtracted using the WISE catalogue.

We outlined three ways to perform RT and describe the gain and potential errors by comparing them. We illustrated that the spatial and opacity information is blended into the derived SFB by the action of multiple scattering and extinction for increasing optical depth. For an example core, we found that the error of the optically thin approximation is about 14% compared to full RT as long as the optical depth for extinction stays below 1. The theoretical coreshine ratios are described and applied for a model core with the properties of L260. They have a singularity for the limiting grain size where the net effect of core extinction and scattering is zero and then decline toward R_theo-values around 1 for grains beyond 2 μm in size.

We described the modeling procedure in detail for L260, also using longer wavelength maps like cold Spitzer data at 8 μm to construct a simple spatial model. The model was used to estimate the optical depth variation. Selecting a region in the ratio map of the two bands that primarily shows optical depth smaller 1, the pixel distribution as a function of the ratio serves to characterize the observed ratios. Then we applied either single-scattering or full RT depending on the assumed maximum optical depth to derive theoretical R values or distributions to compare them to R_theo along with the constraints from the observed maximum coreshine SFBes.

For the cores L260, ecc806, L1262, L1517A, L1512, and L1544, the model is able to account for the observed central SFBes and the ratio maps with maximum grains size around 0.9, 0.5, 0.65, 1.5, 0.6, and >1.5 μm, respectively. The low grain size limit for ecc806 is remarkable given the bright coreshine from this core, and might be related to the low background.

The maximum grain sizes found agree with the findings of earlier studies listed in the introduction, suggesting that coreshine indicates the presence of grains larger than found in the diffuse ISM. No obvious correlation with basic core properties is evident for our sample of nine modeled cores.

The model fails to simultaneously reproduce the SFBes and ratio maps of L1506C, L1439, and L1498. The grain size limits derived from the ratio maps are larger than the limits from the coreshine intensity. Since the coupling to the spatial structure is moderate (τ_ext< 1), this possibly points to assumptions that are too simplified or to the action of an additional component.

Assuming a constant local field with a rising spectral shape on top of the ISRF described by the DIRBE all-sky map does not improve the modeling for L1439. The flat R distribution of the cores L1506C and L1498 gives rise to speculations that both cores host larger grains across the core, either as a primordial component or created in a former more dense and more turbulent phase as discussed in Steinacker et al. (2014b) for L1506C.

Rho Oph 9 is in a complex environment with strong PAH emission. We did not model it in this paper and showed the maps and the ratio map for completeness.

The presented results encourage further exploration of the size distribution parameter space. Observational efforts to gain information at shorter wavelengths will be important for supplying further constraints. Finally, full 3D modeling based on more realistic density structures is required with the option of also fitting thermal emission maps.

Online material

Appendix A: Determination of the off-core surface brightness

In this Appendix, we describe how we have estimated the surface brightness near to but off the core and which variation in the determined ratios can be expected when this off-core measurement is performed differently.

The off-core surface brightness needs to subtracted in Eq. (10) ${\mathit{R}}_{\mathrm{obs}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{I}}^{\mathrm{3.6}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{off}}^{\mathrm{3.6}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}}{{\mathit{I}}^{\mathrm{4.5}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{off}}^{\mathrm{4.5}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}}\mathit{,}$(A.1)since the IRAC measurements are not absolute and contain instrumental, as well as back- and foreground, contributions. Since I_off(x,y) is needed at any PoSky location of the core, but where it cannot be measured, an approximate value needs to be determined for each point (x,y). Ideally, the region where I_off is measured is chosen (i) to be near the core (to represent the surface brightnesses at core location as closest as possible); (ii) to avoid outer core parts (which emission should not be subtracted); (iii) to contain no stellar contributions; and (iv) to enable interpolation of surface brightness variations across the core.

In former work, constant I_off were determined by circular averages around the core (e.g., Nielbock et al. 2012) or choosing the local region near the core with the lowest surface brightness (e.g., Stutz et al. 2009). For cuts through the image, variations have also been used, e.g., by linearly interpolating I_off from locations left and right from the core (e.g., Steinacker et al. 2010; Andersen et al. 2013).

As is visible from the extinction features of cores at 8 μm (see, e.g., Pagani et al. 2010b), almost all cores have a shape deviating from simple spherical symmetry. Judging where outer parts of the core gas are located and where we are looking at variations in the foreground or background gas is therefore difficult, especially for cores near the Galactic plane when the LoS likely crosses other regions.

In the following, we use the core L260 with the strongest coreshine surface brightness in our sample to illustrate how we have performed the off-core measurement and how the results depend on the choice of the background.

In Fig. A.1, the two IRAC images are shown for L260 in the bands 3.6 μm (top) and 4.5 μm (bottom), respectively. The white frame shows the region around the core that is used for the top panels in Fig. 5. As is visible from both images, there is a large scale horizontal gradient across the background of the core.

We have chosen three spots near the core to measure I_off and indicate their location in panel A of Fig. A.2 as numbered white frames. We also give for comparison the frame 0 in which we measured R. In panel B we show the number distributions of pixels as a function of their surface brightness for the entire IRAC image. The distributions contain the stellar sources that are visible in the high SFB wing of the distribution as they are above the mean SFB of each frame: the frame around Core (A), the off-frames 1-3, and the R-measurement frame 0. The gradients visible in the IRAC image results in a mean surface brightness shift from 1 to 3. We therefore interpolate I_off from 1 and 3 which leads to a value close to the mean surface brightness seen in frame 2.

The situation changes for 4.5 μm. Using the same frames as indicated in panel C, the surface brightness distributions in frames 2 and 3 are almost identical. Nevertheless, the gradient between 1 and 3 remains, and we also interpolate I_off at 4.5 μm from both frames.

To estimate the uncertainty in the derived range of observable R, we compare the mean R from interpolating between 1 and 3 and between 2 and 3. Using the mean surface brightnesses of the four frames, we get R = 2.3 for frame 1 to 3 off measurement, and R = 2.41 for frames 2 to 3. We performed this procedure of testing the variation in various frames for all sources discussed here.

Appendix B: Mie calculations

The absorption and scattering coefficients of large grains can be exactly computed for spherical particules only using Mie theory (Bohren & Huffman 1983). Other grain shapes rely on approximate numerical models, or their validity is restricted to small grain sizes.

In this study the water optical constants are extracted from the database of Hudgins et al. (1993). The silicate data are taken from Draine & Lee (1984), and the amorphous carbon data are provided by Zubko et al. (1996).

To simulate porous grains, we use an effective medium formulation where the inclusions are made of vacuum. Amorphous carbon is also considered as inclusions in the silicate matrix.

Consider a particulate composite consisting of a matrix and including various sizes and shapes made of material other than the matrix. Under certain conditions, the composite can be homogenized; i.e., the composite can be replaced by a homogeneous dielectric medium with the same macroscopic electromagnetic response and a certain effective permittivity. Landau & Lifshitz (1960) and independently Looyenga (1965) proposed an effective medium formulation that take inclusion connections for all volume fractionss into account (hereafter LLL model). The effective dielectic function is a volume-fraction weighted-average of the spherical constituents for the composite and is correct to the second order in the differences in permittivities, although dipole-dipole interaction is still not taken into account. For N constituents, the effective dielectric function is ${\mathit{ϵ}}_{\mathrm{eff}}^{\mathrm{LLL}}\mathrm{=}{\mathrm{⟨}{{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{1}\mathit{/}\mathrm{3}}}^{\mathrm{⟩}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{=}{\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{f}}_{\mathit{i}}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{1}\mathit{/}\mathrm{3}}\right)}^{\mathrm{3}}\mathit{,}$(B.1)

where ϵ_i is the dielectric function of the material i, and f_i are the volume fractions of arbitrary shape. Here, N is the total number of inclusions of a given composition. The sum of all volume fractions has to be lower than 1. For two constituents, the formulation is extremely simple: ${\mathit{ϵ}}_{\mathrm{eff}}^{\mathrm{LLL}}\mathrm{=}{\mathrm{[}\mathrm{(}{\mathit{ϵ}}_{\mathrm{inc}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{-}{{\mathit{ϵ}}_{\mathrm{mat}}^{\mathrm{1}\mathit{/}\mathrm{3}}}^{\mathrm{)}}{\mathit{f}}_{\mathrm{inc}}\mathrm{+}{{\mathit{ϵ}}_{\mathrm{mat}}^{\mathrm{1}\mathit{/}\mathrm{3}}}^{\mathrm{]}}}^{\mathrm{3}}\mathit{,}$(B.2)where ϵ_mat is the dielectric function of the matrix. It is clear that the LLL model is symmetric with respect to the constituents. Other formulations of the effective medium dielectric function exist (Maxwell Garnett 1904; Bruggeman 1935), each of them with their own strengths and weaknesses. We chose to use the LLL model for its validity for all volume fractions.

Water ice is assumed to form a mantle on the top of the porous silicate + amorphous carbon core. Absorption and scattering cross-sections, as well as the phase-functions for ice-coated porous spherical grains, are computed using the latest version of the Mie routine for coated spheres provided by Toon & Ackerman (1981).

---

Acknowledgments

J.S., M.A., and W.F.T. acknowledge support from the ANR (SEED ANR-11-CHEX-0007-01). M.J. and V.-M.P. acknowledge support from the Academy of Finland grant 250741. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

References

All Tables

All Figures