© ESO 2018

1 Introduction

The exact composition and size distribution of interstellar dust are still open questions. The properties of dust are known to evolve from diffuse interstellar medium (ISM) to dense regions (Bernard et al. 1999; Ridderstad et al. 2006). It has been assumed that in diffuse regions the properties of dust are uniform. However, the Planck-HFI observations show significant variations in dust spectral energy distribution (SED), in the long-wavelength opacity spectral index all over the sky (Planck Collaboration XXIV 2011; Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014; Planck Collaboration Int. XXIX 2016). Ysard et al. (2015) have shown that the variations of the SED cannot be explained by small changes in the radiation field, but rather they can be explained by variation or evolution of the dust properties.

Evidence of dust grains growing in size in dense regions of the ISM was reported by Stepnik et al. (2001), using dust emission observations of a filament in the Taurus cloud complex covering a wavelength region between 200 and 600 μm. Stepnik et al. (2001) concluded that the observations could only be explained by taking into account larger grains (Rawlings et al. 2005; Ysard et al. 2012, 2013; Martin et al. 2012).

Scattering at near-infrared (NIR) wavelengths can be used to study the composition of interstellar dust. The first observations of surface brightness at NIR wavelengths were reported by Lehtinen & Mattila (1996). Foster & Goodman (2006) showed that it is possible to obtain large maps of the surface brightness, and named the phenomenon “cloudshine”. The history of light scattering observations in dark clouds is covered in more detail by Juvela et al. (2006). The sensitivity to the size of the scattering particles makes observations of scattered light potentially useful in the study of the evolution of grain sizes.

Steinacker et al. (2010) studied molecular cloud L183 and reported an excess surface brightness in the mid-infrared (MIR) towards the dense core. The MIR excess was also discovered in a larger sample of cloud cores by Pagani et al. (2010) and Juvela et al. (2012). The surprisingly high surface brightness in the Spitzer 3.6 μm band was named “coreshine”. It was interpreted as tracing light scattering caused by larger grains than indicated by the classical size distribution presented by (Mathis et al. 1977; hereafter MRN). Andersen et al. (2013) studied the scattered light of cloud core L260, and came to the conclusion that grain sizes up to 1 μm were required in order to explain the shapes of the surface brightness profiles. Other recent studies on coreshine have considered the location on the sky and the strength of the coreshine effect (Steinacker et al. 2014a) and the effect of turbulence as a cause for enhanced grain growth and coreshine (Steinacker et al. 2014b). The detection of coreshine is considered direct evidence of dust grains growing in size in the ISM, because confirming the growth from thermal emission observations requires detailed modelling (but see Ysard et al. 2016).

Coreshine has been detected from numerous sources above and below the galactic plane (Steinacker et al. 2014a), and in some cases near the galactic anti-centre. Detection of coreshine towards the bright galactic plane is difficult due to the strong background combined with the extinction caused by the core (Steinacker et al. 2014a). Additionally, the unknown properties of the local radiation field cause problems in the reliable confirmation of coreshine. In a recent paper, Lefèvre et al. (2016) have studied the possibility of using longer wavelengths, up to 8 μm, to place constraints on dust properties. The authors show that the surface brightness of the scattered light at 8 μm is brighter than was previously thought, and scattering at these longer wavelengths should also be taken into account.

If the grain sizes are originally below 0.2 μm, the grain size should affect the NIR colours. In this study we explore the limitations of using NIR and NIR combined with MIR observations to probe the properties of interstellar dust. Apart from the Spitzer Space Telescope and until the launch of the James Webb Space Telescope, there are no space-borne observatories capable of observing at MIR wavelengths. Since ground-based MIR observation are difficult, the more readily available NIR observations could be a crucial tool in the study of the dense ISM.

Radiative transfer computations can be used to determine how the properties of dust grains and of the local radiation field translate into observed surface brightness. From the modelled surface brightness we can estimate confidence limits for dust parameters, for example, the maximum grain size and the power-law exponent of the size distribution. In this study we use radiative transfer calculations to derive the surface brightness for model clouds in the J, H, K, and 3.6 μm bands. The calculations use a dust model with varying size distributions. We use a one-dimensional spherical model and a three-dimensional ellipsoid cloud model. The results of our radiative transfer computations are analysed with the Markov chain Monte Carlo (MCMC) method to estimate confidence regions for the parameters of the grain size distribution. As an example of real observations and to compare the simulations to the observations, we utilise NIR and MIR observations of the filament TMC-1N (Malinen et al. 2012, 2013) in the Taurus molecular cloud complex, covering the J, H, K, 3.6 μm and 4.5 μm bands. Malinen et al. (2012) studied the possibility of using NIR observations to derive filament profiles and Malinen et al. (2013) studied the possibility of using scattered NIR light on a large scale to study the properties of filamentary structures compared to dust emission studies, for example.

The content of this paper is as follows. In Sect. 2, we give an overview of the radiative transfer and MCMC methods and present the cloud models we use in our study. We present our main results in Sect. 3 and discuss the results in Sect. 4. Finally, in Sect. 5 we summarise our findings. The TMC-1N observations are summarised in Appendix A.

2 Methods

We investigate how the observed surface brightness depends on the local radiation field and the dust properties. In the following we describe the cloud models and the radiative transfer methods used in this study.

2.1 Cloud models

We use a spherically symmetric one-dimensional cloud model and a three-dimensional ellipsoid model. In the one-dimensional spherical cloud model the density distribution is set according to the Bonnor–Ebert model (ξ = 4.5, M_⊙ = 25) (Bonnor 1956; Sipilä et al. 2015). The model optical depth τ is varied by a direct multiplication of the cloud densities. The model is divided into 50 shells such that the innermost 20 shells contain approximately 40% of the mass.

The three-dimensional elliptical cloud model is a prolate ellipsoid with an axis ratio of 3:1, which is discretised onto a Cartesian grid of 68³ cells. The density profile is defined as $${\rho }_i=\exp (-k_1{X}^2-k_2(Y^2+Z^2)),$$(1)

where X, Y, and Z are the Cartesian coordinates of the point i and $k=4\times \ln (2)/L_{\text{ax}}^2$, where L_ax is the length of the axis.

2.2 Radiative transfer

We assume that the surface brightness consists of scattered radiation and of background radiation seen through the cloud. In particular, we assume that there is no notable thermal emission. The assumption follows the arguments provided by Steinacker et al. (2010); integrating the infrared emission model of Draine & Li (2007), we would expect ten times more emission in the 5.8 μm band than in the 4.5 μm band. In observations, including our data on the Taurus filament (see Appendix A), the intensity of the 4.5 μm band is observed to be lower than the intensity of the 5.8 μm band and thus the contribution of emission is not significant. The observed surface brightness excess relative to the background sky is $$I(\lambda )=I_{\text{sca}}(\lambda )+I_{\text{bg}}(\lambda ){\text{e}}^{-\tau }-I_{\text{bg}}(\lambda ),$$(2)

where I_sca(λ) is the intensity of the scattered light, I_bg(λ)e^−τ is the intensity of the background radiation coming through an optical depth of τ, and I_bg (λ) is the intensity of the background seen around the cloud. Thus, we compare the surface brightness of the cloud against the background intensity. In this paper we do not consider the complication that part of the I_bg (λ) originates inregions between the cloud and the observer.

In our simulations, dust properties are defined by the absorption and scattering efficiencies Q_sca, Q_abs, and the asymmetry parameter of the Henyey–Greenstein scattering function g from Weingartner & Draine (2001). Furthermore, we assume a simple power-law size distribution, n(a) ∝ a^−γ, where a is the grain size (with a constant minimum grain size 1.0 nm and a maximum grain size A_max) and a power-law exponent γ that can be varied, compared to the reference case with A_max = 0.25μm and γ = 3.5.

In order to ensure comparability between the models, we normalise the extinction of the J band of both silicate and carbon grains. As a default, we assume that the fraction of the J-band optical depth produced by silicate grains is r_Si = 0.5. To examine the sensitivity to differences in the dust properties, we also test cases with r_Si = 0.3 and r_Si = 0.7. Furthermore, because we assume that emission is negligible, we remove the polycyclic aromatic hydrocarbons (PAHs) from all dust models. To study the effect of increased or decreased scattering efficiency, we scale the albedo α, of the grains by ±10%, and in order to study the effect of asymmetry the value of the asymmetry parameter g is changed by ±10%.

We use the radiative transfer code CRT (Juvela & Padoan 2003; Juvela 2005) to compute surface brightness maps for the J, H, K, and 3.6 μm bands for each combination of A_max and γ, for each band.

2.3 Parameter estimation

The MCMC method is used to draw samples from a probability distribution without having to know the explicit shape of the distribution. This is achieved by a random walk. Based on observations (real or synthetic) and model-predicted surface brightness values, we use the MCMC method to sample the probability distributions of the dust parameters.

In the MCMC procedure, we take a random step in the parameter space and compare the probabilities in the new and old positions. If the difference between the probability of the new and old parameter combinations is larger than ln (R), where R is a random number between 0 and 1, the new parameter combination is accepted, otherwise we keep the old parameters. The current parameter values are saved before repeating the procedure. With enough steps, the saved parameter values trace the probability distributions of the dust parameters.

The radiative transfer computations are performed only for a fixed grid of parameter values, thus we uselinear interpolation to estimate the intensities of scattered light at any value of A_max and γ. The MCMC routine goesthrough the probability distribution step-by-step by taking the interpolated intensity values corresponding to a single value of A_max and γ at a time. The probability p is $$\ln (p)=-\frac{1}{2}{\displaystyle \sum {\left(\frac{I_{\text{obs},i}-[I_{\text{sca},i}{K}_{\text{ISRF}}-I_{\text{BG},i}(1-e^{-{\tau }_i})]}{{\sigma }_i}\right)}^2},$$(3)

where I_{obs, i} are the observed background-subtracted intensities, I_{BG, i} are the unattenuated background intensities, τ_i are the optical depths, σ_i are the assumed observation uncertainties, I_{sca, i} are the scattered intensities, and i runs over the bands. Representing the scaling of the ISRF, we use a fourth parameter K_ISRF to multiply the I_{sca, i} values.

The χ² values are $${\chi }^2={\displaystyle \sum _{i=1}^N{\left(\frac{I_{\text{obs}}-I_{\text{mod}}}{{\sigma }_i}\right)}^2},$$(4)

where I_obs are the “observed” (input) intensities, I_mod are the corresponding intensities predicted by our model, σ_i are the uncertainties, and N is the number of bands.

The models are illuminated either by an isotropic radiation field corresponding to a reference interstellar radiation field (ISRF; Mathis et al. 1983) or by an anisotropic radiation field derived from the DIRBE all-sky maps (see Malinen et al. 2013 for details).

We use flat prior distributions: 0.28–2.5 μm for A_max, 2.0–5.0 for γ, and 0.3–3.0 for K_ISRF. For τ_J, we use a value of 6.0 or 2.0 (corresponding to our reference observations, see Sect. 2.4). The optical depths of the other bands are determined by the extinction curve of the dust model. The radiative transfer calculations were carried out with a fixed external radiation field, but as the scattered light depends on this linearly, the result can be scaled to any value of ISRF. We use the same K_ISRF value for all bands.

In the fitting, we assume a 10% uncertainty for the NIR bands. For the 3.6 μm band the uncertainty is 25%, which takes into account the typically higher uncertainties on the background estimation (see Appendix A). We also examine a case with 20% uncertainty for all four bands. We refer to these two cases as σ_NIR = 10% and σ_NIR = 20%, respectively.

Shown in Fig. 1 is a summary of the work flow from dust parameters to probability distributions.

Fig. 1 Flow chart describing the steps from dust parameters to probability distributions.

2.4 Taurus observations

Our study is based mainly on simulated data. However, as an example of real observations, we analyse surface brightness measurements of a filament in Taurus, TMC-1N. The filament was examined by Malinen et al. (2012, 2013) using the Wide Field CAMera (WFCAM) instrument of the United Kingdom InfraRed Telescope (UKIRT). The observations consist of data in the J, H, and K bands. In addition, weuse data from the Spitzer InfraRed Array Camera (IRAC; Fazio et al. 2004). A short overview of the observation is provided in Appendix A. For a more detailed description of the observations and data reduction see Malinen et al. (2012, 2013).

The twochosen positions A and B of the Taurus filament are shown in Fig. A.1. The surface brightness and optical depth values are listed in Table A.1.

3 Results

Before looking at the synthetic observations in Sect. 3.2, we examine how the dust cross sections and asymmetry parameters depend on the parameters γ and A_max of the grain size distribution.

3.1 Basic effects of grain size distribution

Shown in Fig. 2 are the cumulative sums of the absorption and scattering efficiencies of graphite and silicate grains as a function of grain size. For each grain size the absorption and scattering efficiency is defined as $$\kappa =\pi a^{2}n(a)Q(\lambda ),$$(5)

where a is the grain size, n(a) is the number density of grains of size a, and Q(λ) is the absorption or scattering efficiency of the grains of size a at wavelength λ. Light scattering from grains smaller than ~0.1 μm is insignificant; however, the absorption caused by small grains (a < 0.01μm) can reach ~25% with higher values of γ.

The extinction efficiency, albedo, and asymmetry parameter for pure silicate and graphite grains and for our three models with r_Si = 0.5, r_Si = 0.3, and r_Si = 0.7 are shown in Fig. 3. The Q_ext is almost constant at NIR wavelengths; however, the albedos and asymmetry parameters vary in the ranges ~0.5–0.8 and ~0.3–0.6, respectively.

In Fig. 4, we show the values of g, Q_abs, and Q_sca, and the resulting albedo α and extinction efficiency Q_ext, for different combinations of A_max and γ, and assume r_Si = 0.5. As expected, the extinction decreases and the scattering efficiency relative to absorption increases when moving from J to K band. Furthermore, the maxima of the scattering efficiency Q_sca shift from A_max ~ 1.5μm to A_max ~ 2.25μm. For both J and K bands, the scattering efficiency decreases with increasing γ values.

The differences of the asymmetry parameters indicate that the relative brightness between the J and K bands depends strongly on the location of the cloud with respect to the source of illumination. For example, if an optically thin cloud is seen towards the galactic centre, the J band should be bright relative to the K band; instead, if the cloud is seen towards the anti-centre, the K band should be brighter. The effect is caused by the lower value of the asymmetry parameter of the K band, resulting in stronger back scattering for all dust species with sizes smaller than ~1.2 μm. However, in dense clouds the optical depth of the cloud will have a stronger effect on the relative brightness between the bands. Higher optical depth reduces both the amount of scattered light (relative to dust mass) and the difference between the attenuated and unattenuated background radiation.

Fig. 2 Cumulative sums of absorption and scattering efficiencies for graphite (solid lines) and silicate (dotted lines) grains at J, K, and 3.6 μm bands. The red lines correspond to γ = 2.5, the blue lines to γ = 3.5, and the black lines to γ = 4.5.

Fig. 3 Extinction efficiency Qext, albedo α, and asymmetry parameter g, for different dust mixtures. The optical properties are normalised so that for J band Qext = 1.

Fig. 4 Grain properties as a function of the parameters Amax and γ of the grain size distribution for the J and K bands and for the ratio J∕K. The panels show the asymmetry parameter g; albedo α; and the scattering, absorption, and extinction efficiencies Qsca, Qabs, and Qext, respectively. The plot corresponds to rSi = 0.5, and the colour bars show the values of the dust parameters.

3.2 Synthetic observations

In this section, we present the results of synthetic observations produced by radiative transfer calculations with one-dimensional models, assuming r_Si = 0.5 and using the reference radiation field described in Sect. 2.3.

Figure 5 shows the surface brightness and optical depth profiles computed trough the centre of the model cloud for a dust model with A_max = 1μm and γ = 3.5. The density of the cloud is set so that the average optical depth of the J band computed over the central r∕10 region is τ_J = 2.0 or 6.0, where r is the radius of the model. The Monte Carlo noise in the computed surface brightness and optical depth values is less than ~1%, significantly less than the typical uncertainties on the observations.

For the τ_J = 6 case, the J-band surface brightness at the centre of thecloud is ~30% lower than on the outer edge because of the high optical depth and the resulting saturation of the surface brightness. The K-band intensity is ~20% lower in the centre. The 3.6 μm band is starting to saturate and traces more closely the density profile of the cloud. For the τ_J = 2 case, only the J band shows saturation. In Fig. 6, we show the intensity of the J band (left), the intensity ratios of the K to the J band (centre), and of the 3.6 μm to the J band (right) as a function of γ and A_max. As in Fig. 4, because of the small number of large grains, values of γ > 4.5 produce significantly less intensity regardless ofthe value of A_max. Given that the J-band optical depth is fixed, the J-band intensity reaches its maximum value at A_max ~ 1.1μm and γ ~ 2.22. Compared to the MRN model, the distribution contains more large grains because the slope of the size distribution is not as steep.

The spectral energy distributions with average τ_J = 2 and τ_J = 6 are shown in Fig. 7, calculated with A_max = 1.0μm and γ = 3.5. The optical depth of the J band is varied between 1.4–2.4 and between 4.85–7.85, respectively. The J band begins to saturate when the optical depth approaches ~2, as can be seen from the panel on the left. With higher optical depth (right panel) the intensity of the J, H, and K bands decreases when the optical depth increases, but the intensity of the 3.6 μm band increases.

The surface brightness values of the J and K bands and the ratio of the J to the K band as a function of the optical depth of the J band are shown in Fig. 8. The curves correspond to different A_max and γ values. We show curves for a fixed value of γ = 3.5 with A_max varying from 0.1 to 1.8 μm, and for a fixed value of A_max = 1.0μm with γ varying from 2.0 to 4.5.

When A_max is constant, the peak intensity is reached at τ_J ~ 3.8 for the J band and at τ_J ~ 5.0 for the K band, regardless of the value of γ. Comparing the ratios of the K to the J band intensities, it is possible to determine the value of γ only if τ_J < 6 or γ > 4.

With a constant γ, the peak intensity of the J band shifts towards higher optical depth when A_max increases, whereas the peak intensity of K band shifts towards lower optical depth when A_max increases. It is possible to determine the value of A_max from the ratio of K and J band intensities in almost all cases, except if A_max < 0.2 and the optical depth is low, τ_J < 2.

Fig. 5 Results for one-dimensional radiative transfer models with rSi = 0.5, Amax = 1.0μm, and γ = 3.5 with averageτJ = 2 (left) and τJ = 6 (right). The panels show the optical depth of the J band and the surface brightness of the J, K, and 3.6 μm bands. The surface brightness values and optical depth have been normalised.

Fig. 6 Left: intensity of the J band as a function of γ and Amax. Centre and right: the intensity ratios of the K to the J band, and of the 3.6 μm to the J band, respectively. A one-dimensional model is used with rSi = 0.5 and τJ = 6. The colour bars show the intensity and intensity ratios.

Fig. 7 Spectral energy distributions resulting from a computation with rSi = 0.5. The maximum grain size is set to Amax = 1.0μm and γ = 3.5.

Fig. 8 Surface brightness of the J and K bands, and the ratio of the K to the J band as a function of optical depth. The colours of the curves correspond to different values of Amax and γ. For the first column γ = 3.5, and for the second column Amax = 1.0μm.

3.3 MCMC on synthetic observations

To determine the accuracy to which dust parameters can be deduced from the observed intensities, we apply the MCMC procedure to synthetic observations of one-dimensional and three-dimensional clouds where the true dust and ISRF parameters are precisely known. In rough correspondence to the Taurus observations discussed in Sect. 2.4, the one-dimensional synthetic observations are created using randomly generated values values τ_J = 6.04, K_ISRF = 1.53, A_max = 1.11μm, and γ = 3.49.

In Fig. 9, we show the cross sections of the χ² distribution using the synthetic observations. The white star marks the location of the minimum at A_max = 1.09, γ = 3.53, K_ISRF = 1.58, and τ = 6.36. The small difference in the input parameters is caused by the use of discretised parameter grids and the shallow χ² valley running almost parallel to the τ_J axis. Each parameter grid has 50 points, for A_max from 0.25 to 2.5 μm, for γ from 2.5 to 5.0, for K_ISRF from 0.3 to 3.0, and for τ from 4.8 to 7.4.

Figure 10 shows the marginalised probability distributions of a one-dimensional model cloud with τ_J = 6. An isotropic radiation field is used in the scattering computations. The distributions are normalised and we only compare the relative probabilitiesof the different parameter combinations. Figure 10 indicates a strong dependence between γ and K_ISRF. The dependence between parameters is also seen in the (γ, τ) projection as a strong cut-off at high γ. Thus, to be able to determine either the γ or K_ISRF, the values of the other parameters must already be known with good precision. The four-dimensional parameter space is not easily represented with two-dimensional correlation plots. The maximum of the marginalised probability and the projected location of the minimum in the four-dimensional χ² are both shown.

The assumed error estimates of the intensity values have little effect on the results, apart from the broadening of the distributions. Small variations in the intensities of individual bands do affect the location of the χ² minimum, as can be seen from the locations of the green and red symbols in Fig. 10. Even a moderate change of ±10% in the intensity of one band can change the recovered χ² minimum values of γ and K_ISRF by up to ~60%. In the fit a change of optical depth can be compensated by a change of the radiation field intensity. For further computations we use a constant value of τ_J = 6 or τ_J = 2, since it is possibleto derive a reasonable estimate for the optical depth using other methods, for example with observations of background stars (Lombardi & Alves 2001).

Using an anisotropic radiation field in the scattering computations (Fig. B.1) produces well-constrained distributions for the parameters γ and K_ISRF, but the A_max distribution is considerably broader and shifted to higher parameter values.

The analysis of the surface brightness data depends not only on the line-of-sight optical depth, but also on the cloud morphology. We derive synthetic observations from a three-dimensional prolate ellipsoidal model cloud (see Sect. 2.1). The dust parameters were set to A_max = 1.18, γ = 3.67, and K_ISRF = 1.2 and an anisotropic radiation field corresponding to the Taurus case (see Malinen et al. 2013) is used in the scattering computations. The optical depth along the line of sight is set to τ_J = 6.0.

For the analysis of these synthetic observations, the relation between dust parameters and the surface brightness is needed. This can be derived with radiative transfer modelling, but it also depends on the assumptions of the cloud shape. We compare two cases where the relationship is derived with radiative transfer modelling; a three-dimensional prolate ellipsoidal cloud is viewed either along its major axis or along its minor axis (see Fig. 12). In the following we refer to these as the major axis and minor axis radiative transfer models, respectively. In these radiative transfer calculations the model cloud has an axis ratio of 3 and the line-of-sight optical depth is set to τ_J = 6 in both cases. When the radiative transfer model is used to calculate surface brightness along its minor axis, it suffers from more extinction in the direction perpendicular to the line of sight (where the full optical depth through the cloud is τ_J = 18). This affects the radiation field inside the radiative transfer model and thus also leads to a different mapping between the dust parameters and the predicted surface brightness. This is expected to bias the estimates derived for the dust parameters. In the scattering computations we use the anisotropic radiation field. The synthetic observations correspond to the case where the radiative transfer model is viewed along its minor axis. The results with σ_NIR = 10% are shown in Fig. 11, which shows that the assumptions of the radiative transfer models have a clear effect on the estimated dust parameters.

The one-dimensional probability distributions (the minor axis case) of both γ and K_ISRF are not concentrated around the parameter values that were used to derive the synthetic observations. The maximum of the two-dimensional probability distribution of γ and K_ISRF are centred around γ = 4.1, K_ISRF = 1.0. However, the maxima of the one-dimensional A_max distribution and the χ² minimum match the input parameter values closely, but the mode of the marginalised probability distributions significantly differs from the location of the χ² minimum. In the major axis case, the probability of larger grains increases significantly as the A_max distribution is shifted to A_max > 1.5μm and the peak of the γ distribution is shifted from γ ~ 4.3 to γ ~ 4.1. On the other hand, the strength of the radiation field decreases by ~30%.

Increasing the uncertainty on the surface brightness observations, σ_NIR = 20%, only broadens the resulting probability distributions. Similarly, changes in the asymmetry parameter g do not produce noticeable changes compared to Fig. 11. If K_ISRF is fixed to a constant value (Fig. 13), the marginalised probability distributions and the χ² minimum match the used parameter values with much higher precision.

Shown in Fig. 14 are the results when the albedo is changed by ±10% in the scattering simulations; however, for the parameter estimation we use the unmodified albedo. Changes in the albedo have a clear effect on the parameter distributions. Increasing the albedo of the grains shifts the peak of the A_max distribution from A_max ~ 1.0μm to A_max ~ 1.5μm, and decreases the probability of γ < 3.5. Thus, the maximum grain size increases but the relative number of large grains decreases. On the other hand, with higher albedo, the strength of the required radiation field decreases from K_ISRF ~ 1.3 to K_ISRF ~ 0.7.

Fig. 9 Cross sections of the χ2 distribution through the location of the χ2 minimum for the synthetic observations. The white star indicates the χ2 minimum and the dashed lines show the dust parameter values used to simulate the synthetic observations. The colour bars correspond to the χ2 distribution.

Fig. 10 Marginalised one-dimensional probability distributions (panels A, F, K, and P) and two-dimensional projections of the probability of dust parameters for τJ = 6. The white star indicates the projected position of the χ2 minimum. The other symbols are the χ2 minima assuming a 10% higher or lower intensity for one of the bands, marked with green and red symbols, respectively. The circles correspond to J band, the diamonds to H band, the squares to K band, and the triangles to 3.6 μm band. The black contour shows the 1σ of the projection. The colour scale shows the normalised probability.

Fig. 12 Schematic of the major axis and minor axis radiative transfer models (right) used to analyse the synthetic observations(left).

Fig. 11 Marginalised probability distributions of dust parameters for an ellipsoidal cloud model. In the radiative transfer modelling (see text) the cloud is viewed either along its minor axis (red lines) or its major axis (purple lines). For both cases, the optical depth along the line of sight is τJ = 6 and the dashed lines show the input values used in the computations. The stars indicate the projected positions of the χ2 minima and the black contours show the 1σ and 2σ of the projections. The greyscale map corresponds to the red lines.

Fig. 13 As in Fig. 11, but the KISRF is set to a constant value of 1.2.

Fig. 14 Marginalised probability distributions of dust parameters for τJ = 6. As in Fig. 8, but with the grain albedo in the radiative transfer calculations changed by −10% (red lines) or +10% (purple lines). The stars indicate the projected positions of the χ2 minima, and the greyscale map corresponds to the red lines. The dashed lines show the input values used in the computations.

3.4 Taurus TMC-1N observations

In this section, we analyse the NIR observations of the Taurus filament described in Appendix A. We use the minor axis cloud model described in Sect. 3.3. For the optical depth of the J band we use values of τ_J = 6.0 and τ_J = 2.0. The surface brightness and optical depth values derived for the two cases are shown in Table A.1. Unless otherwise specified, we use the r_Si = 0.5 ratio.

3.4.1 MCMC results, τ_J = 6 position

The results of the MCMC simulations assuming σ_NIR = 10% or σ_NIR = 20%, where only the A_max and γ were varied during the scattering computations, are shown in Fig. 15. A clear peak can be seen at A_max ~ 1μm and a broad plateau around γ ~ 3.5. An increase in the σ_NIR produces smoother parameter distributions. The σ_NIR has no effecton the K_ISRF distribution; however, with higher σ_NIR the A_max distribution is pushed to higher values, A_max≳1.5μm.

The location of the χ² minimum matches the maximum of the marginalised probability distributions quite closely. However, the minimum seems to follow mainly the γ distribution as the minimum does not closely match the projections with K_ISRF or A_max. When only J, H, and K band observations are used, A_max = 0.8μm and γ = 3.2 (Fig. 16). The probability of A_max > 2.0μm has also decreased considerably.

According to Fig. 17, decreasing the albedo by 10% produces a sharper distribution around A_max ~ 1μm and increases the probability of small grain sizes, A_max < 1.3μm and γ < 3.6. A 10% higher albedo will decrease the estimated ISRF strength by ~20%. Furthermore, a higher albedo will also increase the distance between the χ² minimum and the peak of the marginalised probability. A larger σ_NIR will result in a broader and flatter parameter distributions (Fig. C.1).

Varying the asymmetry parameter g by ±10% produces similar results. Although, the projected position of the χ² minimum is shifted relative to the maxima of the marginalised probability distributions (Fig. 15). With larger σ_NIR, the probability distributions are again wider and flatter and their maxima match the χ² minima better.

Fig. 15 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. Forthe surface brightness data we assume σNIR = 10% (red lines) and σNIR = 20% (purple lines). The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines, and the contours show 1 σ and 2 σ for the two models.

Fig. 16 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. Forthe computations, only the J, H, and K bands were used. The white star indicates the projected position of the χ2 minimum and the black contours show the 1σ of the projection. The colour scale shows the normalised probability.

Fig. 17 As in Fig. 15, but the grain albedo in the radiative transfer calculations is changed by −10% (red lines) or +10% (purple lines). The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.

3.4.2 MCMC results, τ_J = 2 position

The results of the MCMC computations of the TMC-1N τ_J = 2 position areshown in Fig. 18. With lower optical depth the A_max distribution is limited from above only by the imposed upper limit A_max = 2.5. Likely for the same reason, the γ values are all concentrated around 3.5. Compared to the TMC-1N τ_J = 6 position, the ISRF estimate is lower with K_ISRF ~ 1.0.

Assuming larger uncertainties for the surface brightness data has a strong effect (Fig. 18B) shifting the distribution from γ ~ 3.0 to γ ~ 2.4, and decreasing K_ISRF from ~1.0 to 0.8. This suggests some tension between the observations and the model; for a more diffuse line of sight the dust grains should be less evolved with lower A_max and higher γ. In contrast, there is almost no effect on the A_max, with the exception of a slight broadening of the distribution. The lower optical depth has led to smaller offsets between the χ² minima and the maxima of projected probability.

The computations carried out by varying the albedo of the dust grains (Fig. C.3) produce similar effects to those in Fig. 17. Decreasing the albedo increases the required radiation field strength and increases the probability of large grains. Changing the asymmetry parameter does not effect the derived parameter distributions.

Fig. 18 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 2 position assuming σNIR = 10% (red lines) and σNIR = 20% (purple lines). The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.

3.4.3 Grain composition

In the previous sections, we used models that have equal J-band optical depth for silicate and carbon grains. In order to study how the chemical composition of the grains affects the scattering properties, we examine two alternative dust models with r_Si = 0.7 and r_Si = 0.3 ratios.

The model with r_Si = 0.7 prefers a lower value of γ than the model with r_Si = 0.3 (Fig. 19). The shapes of the two distributions could explain the broad γ distribution seen in the model computed with synthetic observations (Fig. 10). With r_Si = 0.7 the γ distribution is skewed towards γ < 3.5, but with r_Si = 0.3 the distribution is skewed towards γ > 3.5, thus a model with r_Si = 0.5 will have a broader γ distribution.

The A_max distributions of the two models have rather similar shapes, with a plateau at ~1.5 μm, although the distribution is broader for the r_Si = 0.3 case. The model with r_Si = 0.7 produces a sharp and narrow distribution around K_ISRF ~ 0.7, whereas the computation with r_Si = 0.3 is centred around K_ISRF ~ 1.0, with a tail towards higher values. The differences can be explained by carbon having a higher overall absorption efficiency, whereas silicates have higher scattering efficiency. With r_Si = 0.3, the χ² minimum does not match peak of the marginalised probability.

Fig. 19 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position with a rSi = 0.3 (red lines) and rSi = 0.7 (purple lines), respectively. For both computations σNIR = 10% is assumed. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.

4 Discussion

We have used radiative transfer computations to study how changes in the properties of interstellar dust grains, maximum grain size, and size distribution affect the observed scattered light at NIR and MIR wavelengths. We show that the scattered light can be used to constrain dust parameters. In our radiative computations we assume that the observed surface brightness consists of scattered radiation and of background radiation seen through the cloud. In particular, we assume that the thermal emission is negligible. We use a spherical one-dimensional cloud model and a three-dimensional ellipsoid model which are illuminated by an isotropic interstellar radiation field, or by an anisotropic radiation field. The properties of the dust grains are based on the model by Weingartner & Draine (2001).

In our simulations we have only used one of the Spitzer bands; however, as discussed by Lefèvre et al. (2016), light scattering even at 8 μm may not be negligible. Careful modelling of all Spitzer bands or even longer wavelengths, for example using far-infrared emission, can be used to place additional constraints on the dust parameters.

4.1 Constraints on dust parameters

Based on our simulations, strong correlations between the dust parameters are evident as changes in one parameter affect the other parameters. The K_ISRF and τ_J cannot be constrained simultaneously, as changes in the optical depth can be compensated by changes in the strength of the radiation field. Similarly, increasing the K_ISRF will shift the γ distribution to lower values as in many of our simulations both parameters tend to have long tails. On the other hand, the colour of the NIR-MIR spectra is affected by changes in the A_max and γ. Furthermore, assumptions of the morphology of the cloud can change the estimated dust parameters significantly.

To categorise the parameter distributions, we use relative uncertainty $$r(x)=\frac{Q_3-Q_1}{1.349\times Q_2},$$(6)

where x is the dust parameter, Q₁ and Q₃ are values of the first and third quartile, Q₂ is the median, and the factor 1.349 is the scaling between the interquartile range Q₃ –Q₁ and the standard deviation in the case of a normal distribution. For the probability distributions that are derived using the synthetic observations, we define the difference between the median of the distribution and the correct value of the dust parameter $$\Delta (x)=\frac{Q_2-x}{x}.$$(7)

In the computations with the synthetic observations, the probability distributions of A_max and K_ISRF have a clear peaknear the correct values used to derive the observations. For the model shown in red in Fig. 10, Δ (A_max) = 13% and r(A_max) = 25%, and Δ (K_ISRF) = −1.7% and r(K_ISRF) = 61%. However, the γ distribution of the model does not have a clear peak around the true value, however, the resulting Δ (γ) = −0.5% and r(γ) = 24%. The difference between assuming a major axis or minor axis radiative transfer model can change the estimated number of large grains by changing the A_max by up to ~30% and γ by up to ~10%, and can affect the estimated K_ISRF by up to ~30%.

Small variations in the optical properties of the dust grains are reflected to the derived dust parameters. If the albedo of the grains is increased by10% in the scattering simulations the probability of larger grain sizes increases as seen in Fig. 14. The peak of the A_max distribution shifts from A_max ~ 1μm to A_max ~ 1.5μm, and the probability of grain sizes in excess of 2 μm increases. However, the increase in albedo in the simulations also increases the probability of γ > 3.5 and decreases the strengthof the required radiation field. Thus, although the maximum grain size increases, the relative number of large grains decreases. Similarly, Δ(A_max) = 27% and r(A_max) = 26%, whereas for the γ distribution Δ(γ) = 5.2% and r(γ) = 18% and for the K_ISRF distribution Δ(K_ISRF) = 5.8% and r(K_ISRF) = 71%. Decreasing the albedo by 10% results in deviations between the correct values and the median values of Δ (A_max) = 1.4%, Δ (γ) = −0.8%, and Δ (K_ISRF) = 14% and the relative uncertainties are r(A_max) = 23%, r(γ) = 28%, and r(K_ISRF) = 38%.

The changes in the shape of the A_max and γ distributions are related to the way radiation is scattered by the grains. A higher albedo decreases the absorption efficiency of the dust grains resulting in a spectra that is bluer than the reference spectra, thus larger grain sizes are needed. However, because of the increased scattering efficiency, there are fewer large grains. Shown in Fig. 20 are the spectra computed from three different dust models by varying either the A_max or γ parameter. Changes in the A_max have a stronger effect on the resulting intensities than changes in γ.

An error of 10% in the intensity of a single band (see Fig. 10) affects the derived value of A_max by <10%, but the derived value of K_ISRF can change by up to 60%. Uncertainty on the observed intensity of the H or K band affects the derived parameter values more than changes in the J or 3.6 μm bands.

Many of our MCMC computations show a high probability for grain sizes in clear excess of 2 μm, but the estimated size distribution can also be changed by using different types of grains. We have used a simple model with bare silicate and carbon grains, thus in order to increase the surface brightness of scattered light at longer wavelengths the maximum grain size has to be increased. However, the scattering properties of the dust grains can be changed by including ice mantles or coagulation, or by changing the surface composition of the grains (e.g. hydrogen-rich carbon mantles), thus increasing the scattering efficiency without increasing the grain sizes (Ysard et al. 2016; Jones et al. 2016).

Coagulation and formation of ice mantles will change the size distribution, decreasing the number of small grains. As seen in Fig. 2, the small grains can cause a significant amount of absorption, and thus can reduce the number of photons available forscattering. We ran separate computations where we removed all grains smaller than 20 nm from the dust model. The results of the simulations, using the synthetic observations of the minor axis case (see Fig. 11) and assuming τ_J = 6, are shown in Fig. 21. Compared to the minor axis case (red lines), removing small grains from the model that is used to fit the synthetic observations (purple lines) has an effect on the resulting marginalised probability distributions, but it is not large. The peak of the γ distribution is shifted to a higher value of γ = 4.5 and the K_ISRF distribution has an upper limit of K_ISRF = 2.

The changes in A_max and γ affect the fraction of mass in large grains (Fig. 22) $$f_{\text{BG}}=\frac{m(a>0.25\mu \text{m})}{m},$$(8)

where m(a > 0.25μm) is the mass of grains larger than 0.25 μm and m is the total mass. In all cases, the mass fraction in large grains increases sharply when the value of γ decreases. In the case of lower optical depth, the f_BG is higher compared to the τ_J = 6 case. However, based on the MCMC results, the probability that a diffuse line of sight would have more large grains than an optically thicker line of sight is p(f_BG(τ_J = 6) > f_BG(τ_J = 2)) = 25%. The result is not statistically meaningful, thus the evidence for a high concentration of large grains towards the τ_J = 2 line of sight is not very strong.

In many of our computations, the estimated K_ISRF has a clear peak around ~1.1, but K_ISRF could also be significantly higher. However, if the radiation field and optical depth can be constrained, for example setting a constant value for K_ISRF and τ_J (Fig. 13) results in well-constrained parameter distributions for both A_max and γ with Δ (A_max) = 3.4%, r(A_max) = 20% and Δ (γ) = 0.2%, r(γ) = 3.2%. Assuming a prior 20% uncertainty on K_ISRF and τ_J, results in A_max and γ distributions that are broader (Δ(A_max) = 4.2%, r(A_max) = 18% and Δ (γ) = 7.9%, r(γ) = 17%, respectively).

The analysis of the computations where the synthetic observations are used shows that on average the r(A_max, γ) ~ 25%, whereas the r(K_ISRF) is ~45% on average. However, the r(K_ISRF) can reach values as high as ~70%.

As seen in the synthetic observations, the optical depth of the cloud is not well constrained by studying only the scattered light. This is related to the saturation of the surface brightness and to the degeneracy between radiation field intensity and dust opacity because changes in opacity can be negated with changes in the strength of the radiation field. When possible, the analysis should include external constraints on the opacity or the radiation field that can be obtained from dust emission or directly from extinction measurements.

If the radiation field strength and the optical depth of the cloud can be constrained, the surface brightness ratio of the J to the K bands can be used to place useful constraints on the dust properties. However, for diffuse regions, the location of the cloud with respect to the source of illumination should be taken into account since differences in the asymmetry parameter can change the relative brightness between the J and K bands; however, the effect of the asymmetry parameter also depends on the optical depth of the cloud.

Fig. 20 Near-infrared spectra derived from different dust models; a constant γ = 3.5 is used (left panel) and Amax = 1.0μm (right panel). The red and green lines are models where the albedo of the dust grains has been changed during the scattering computations by −10% (red) or by +10% (green), and the black line is a model where only the Amax and γ were varied during the scattering computations. Each symbol (circle, triangle, square) corresponds to a different Amax and γ value in the left and right panels, respectively.

Fig. 21 Marginalised probability distributions of the dust parameters for the minor axis case (red lines) (see Fig. 11), and the case with grains smaller than 20 nm removed from the dust model used to fit the synthetic observations(purple lines).

Fig. 22 Massfraction of grains larger than 0.25 μm for τJ = 6 (red lines) and τJ = 2 (purple lines). For panels A–C, the Amax and γ values of the model with rSi = 0.5 and σNIR = 10% are used. Forpanel D, we use the previous model but assume σNIR = 20%. For panels Eand F, we use the models with the albedo of the grains decreased by 10% and the asymmetry parameter decreased by 10%, respectively. The greyscale colour maps correspond to the red lines and the contours show 1σ and 2σ.

4.2 Limitations of the models

We assume that the background in Eq. (2) is truly a background and that there is no emission from between the cloud and the observer. Emission originating between the cloud and the observer causes error in the assumed transmitted radiation I_bg e^−τ and thus also in the deduced colour of the scattered component. A more detailed radiative transfer study is needed to explore the effects arising from the surrounding medium. In our simulations we have used smooth spherical and ellipsoidal cloud models, but have not considered effects of small-scale density inhomogeneities that could affect the way radiation penetrates the cloud. Furthermore, we have used a constant scaling for the strength of the radiation field K_ISRF, thus we have not considered possible changes in the shape of the background spectra.

The K_ISRF distribution tends to have a long tail towards higher values, while the power-law exponent of the grain size distribution γ has a tail towards lower values. The tails are probably caused by the chosen prior distributions of the parameters. Allowing the K_ISRF to have values higher than 3.0, would likely allow the γ distribution to reach even lower values.

The position of the χ² minimum does not usually coincide exactly with the maximum of the marginalised probability. The three-dimensional function χ² (A_max, γ, K_ISRF) is complex, and the minimum is shallow. Furthermore, the values of the χ² function are computed from simulations that have noise, thus the actual location of the computed χ² minimum can be affected by noise.

The marginalised probability distributions are projections of a three-dimensional distribution. Thus, it is possible that the χ² minimum of the three-dimensional distribution does not match the maximum of the projected two-dimensional probability. Furthermore, using higher error estimates for the surface brightness data can have a significant effect on the location of the χ² minima. The instability of the χ² minima seems analogous with the colour temperature and spectral index plane, T_C -β_spec. As discussed by Juvela & Ysard (2012), even a slight change in the weighting of the frequency points or in the noise can significantly shift the location of the global minimum. High uncertainties make the χ² minima shallow and, together with noise, contribute to ambiguity in the actual location of the χ² minima.

4.3 Taurus observations

As an example, we have analysed observations of a filament in the Taurus molecular cloud, TMC-1N. Our cloud model is an ellipsoid with a 3:1 ratio of the major to the minor axes, illuminated by an anisotropic radiation field derived from the DIRBE observationsand viewed along one of the minor axes. The optical depth of the model cloud along the line of sight is set to τ_J = 6 or τ_J = 2. Based on our analysis for the TMC-1N τ_J = 6 position, a maximum grain size of A_max ~ 1.2μm and γ > 4 are very likely, with an average relative uncertainty of r(A_max) = 36%, r(γ) = 22%, and r(K_ISRF) = 38%. However, there is a considerable possibility of grains larger than A_max > 2.0μm. The mass weighted average grain size is ⟨a_m⟩ = 0.113, which is larger than the MRN value of ⟨a_m⟩ = 0.089. The anisotropic radiation field derived from the DIRBE observations is sufficient to reproduce the observed surface brightness values.

In Fig. 23 a cut of the marginalised probability distribution is shown around the χ² minimum for the TMC-1N τ_J = 6 position with r_Si = 0.5 and assuming σ_NIR = 10% (see also Fig. 15). In all of the panels, the χ² minimum matches the maximum of the marginalised probability.

If we assume larger observational uncertainties (σ_NIR = 20%) the resulting parameter distributions are smoother and broader, with an average relative uncertainty of r(A_max) = 26%, r(γ) = 24%, and r(K_ISRF) = 30%. However, the increased σ_NIR of both TMC-1N positions (τ_J = 6 and τ_J = 2) can also change the interpretation of the dust properties. For example, with a lower σ_NIR the A_max distributions in Fig. 15 show a clear peak around ~1 μm, a tail towards higher parameter values, and a r = 37%. Adopting a higher σ_NIR shifts the peak of A_max distribution to higher value (~1.5 μm), and causes a significant increase in A_max > 1.5μm, and decreases r(A_max) to 25%. With lower uncertainties the γ distribution has a plateau or a bump around 3.5 and r(γ) = 22%, whereas with higher σ_NIR the distribution is smoother and r(γ) = 24%.

The 3.6 μm band will place more constraints on the probability distributions, significantly increasing the probability of A_max > 1.5μm with r(A_max) = 30%, r(γ) = 25%, and r(K_ISRF) = 56%. However, caution should be taken with the 3.6 μm band. To derive reliable estimates for the dust parameters the fraction of emission in the observed surface brightness should be known to a reasonable accuracy.

With lower optical depth (i.e. τ_J = 2.0), the parameter distributions are significantly narrower with average r(A_max, γ, K_ISRF) < 10%. The γ values are lower than ~3.8 and the values of A_max are concentrated to values higher than 2 μm. However, the opposite would be expected, since for a more diffuse line of sight, the dust grains are thought to be less evolved, and thus lower values of A_max and higher values of γ would be expected. Increasing the σ_NIR pushes the γ distribution to a lower value (γ ~ 2.4). The change in the γ distribution is likely related to the shape of the χ² space, and changes in the assumed error estimates can also strongly affect the way the priori affect the posterior probability distributions.

The K_ISRF distribution is slightly lower in the case where τ_J = 2.0, although the radiation field there should be less attenuated. A region with lower optical depth could be used to constrain the strength of the radiation field. For case τ_J = 2.0, the intensity of the 3.6 μm band is low, which introduces additional uncertainty on the derived dust parameters.

Fig. 23 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. The white star indicates the projected position of the χ2 minimum. Panels A, D, and G and C, F, and I show the marginalised probability integrated along the third parameter axis within ±5% of the χ2 minimum, respectively. Panels B, E, and H show the marginalised probability on the χ2 minimum.

5 Conclusions

We have used radiative transfer computations of near-infrared scattered light with the Markov chain Monte Carlo method to study the optical properties and size distributions of interstellar dust. We have examined how NIR and MIR observations of scattered light can be used to constrain properties of the dust grains and the radiation field illuminating dense clouds.

The main results of our study are the following:

• Our analysis of the synthetic observations indicates that for surface brightness measurements with ~10% errors, the relative uncertainty on the main dust parameters A_max and γ are ~25%. The relative uncertainty of K_ISRF is higher ~45%. The difference between the median of the posterior probability distribution and the correct parameter value is typically <15%.

• Tests with ellipsoidal models with aspect ratios of 1:3 showed that the assumed morphology of the cloud can change the estimated K_ISRF and A_max by up to ~30% and the estimated value of γ by ~10%.

• An error of 10% in the surface brightness of one band does not significantly affect the derived value o A_max. However, the uncertainty affects the derived values of γ and K_ISRF by up to 30% and 60%, respectively, depending on which of the wavelength bands has the uncertainty. Uncertainty on the H or K band has a stronger effect on the derived dust parameters than on the J or 3.6 μm band, as they affect the shape of the spectra more.

• The strength of the radiation field K_ISRF and the optical depth τ_J cannot be well constrained simultaneously. However, other independent methods can be used to derive constraints for τ_J, for example by extinction studies. If K_ISRF and τ_J are known a priori to an accuracy of ~20%, the uncertainty on A_max and γ is decreased to ~18%.

• The prior used in our MCMC computations cause long tails in the γ and K_ISRF distributions.

• In our tests, when synthetic observations were analysed with radiative transfer models with 10% higher albedo, the A_max and γ distributions shift to higher values (A_max > 1.3μm and γ > 3.5). However, the required radiation field strength decreases to K_ISRF ~ 0.8.

• The χ² minima and the maximum of marginalised probability do not correspond exactly to the used parameter values in most of our computations. The three-dimensional function χ²(A_max, γ, K_ISRF) is complex, and the minimum of the χ² distribution can be very shallow.

• Considering the TMC-1N observations for a line of sight with τ_J = 6, a maximum grain size A_max > 1.5μm and a size distribution γ > 4.0 have high probability, with an average relative uncertainty of ~25%. However, there seems to be a strong peak around A_max = 1μm. The background radiation field derived from the DIRBE observations is sufficient to reproduce the observed surface brightness with K_ISRF ~ 1 and the average relative uncertainty is ~30%. The mass weighted average grain size is ⟨a_m⟩ = 0.113, compared to the value of ⟨a_m⟩ = 0.089 given by the MRN distribution.

• We have used a simple dust model by Weingartner & Draine (2001) that only includes bare silicate and carbon grains. Using a more refined model with additional properties like ice mantles will change the scattering properties of the grains and can result in smaller grain sizes.

Acknowledgements

MS and MJ acknowledge the support of the Academy of Finland Grant No. 285769. JMa acknowledges the support of ERC-2015-STG No. 679852 RADFEEDBACK.

Appendix A Near-infrared observations

The near-infrared observations in the J, H, and K band, corresponding to 1.25, 1.65, and 2.22 μm, respectively are from the Wide Field CAMera (WFCAM) instrument of the United Kingdom InfraRed Telescope (UKIRT). The images are centred at RA(J2000) 4 h39 m36 s, DEC(J2000)+26°39′32″ covering an area of 1° × 1° corresponding to approximately (2.4 pc)². The data reduction is described in Malinen et al. (2013), but in addition we carried out the subtraction of point sources.

The 3.6 μm and 4.5 μm band observations are from the Spitzer InfraRed Array Camera (IRAC; Fazio et al. 2004). The images were obtained from the NASA/IPAC Infrared Science Archive (IRSA). The images are centred on the TMC-1N filament in the Taurus molecular cloud complex, the coordinates of the centre position are RA(J2000) 4h37m45s, DEC(J2000)+26°57’24”. The data (observation numbers 11230976 and 11234816) are from the Taurus Spitzer legacy project (PI D. Padgett). We used the level-2 data.

For background subtraction we use images from the DIRBE satellite. The DIRBE images were obtained from the Sky View virtual observatory. For our data analysis, we only needed the first three DIRBE bands, which correspond to wavelengths of 1.25, 2.2, and 3.5 μm. We used the Zodi-Subtracted Mission Average (ZSMA) images, from which the modelled interplanetary dust signal has been subtracted. We have rescaled the DIRBE magnitudes to 2MASS magnitudes; the first two bands were scaled using the filter conversion described by Levenson et al. (2007) and the third band was scaled using the conversion factors given by Flagey et al. (2006).

The DIRBE bands 1, 2, and 3 can be used directly, after the conversion to 2MASS magnitudes, to specify the intensity at the corresponding wavelengths 1.2, 2.2, and 3.5 μm. The H-band intensity was computed by multiplying the average of the J- and K-band intensities with the intensity ratio I_H∕ < I_J, I_K > derived from the Mathis et al. (1983) ISRF model.

The background intensities I_BG for the J, H, and K bands were adopted from Malinen et al. (2013), who used DIRBE and WISE observations to derive the intensities. However, the value of the background intensity of the 3.6 μm band indicated by Spitzer (0.1315 MJy sr⁻¹) is higher than the value derived by Malinen et al. (2013; 0.0765 MJy sr⁻¹). Thus, we have used the average of the Spitzer and DIRBE values. The background indicated by Spitzer was obtained by taking an average value over a section in the 3.6 μm map close to the filament (RA(J2000) 4h41m30s, DEC(J2000)+26°25′0″ covering an area of ~30″ × 30″) and the value indicated by DIRBE was adopted from Malinen et al. (2013). The obtained background intensities are listed in Table A.1. In addition to deriving the background intensities, we have used the DIRBE ZSMA maps to produce an anisotropic radiation field for our radiative transfer computations.

In order to investigate the intensity of the diffuse signal, point sources must be subtracted. We have used both SExtractor (Bertin & Arnouts 1996) and PSFex (Bertin 2011) to subtract the point sources from the Spitzer and WFCAM images. The subtraction was done in three steps. In the first step, we used the Spitzer and WFCAM images as input for SExtractor. For the initial run, we wanted to detect only the brightest point sources, and so we used a high detection threshold

corresponding to 15 times the estimated noise in the images. The size of the background mesh was set at 64 pixels and the photometric parameters were set to the correct values for the corresponding Spitzer and WFCAM bands. With these parameters, SExtractor produces a catalogue of the brightest point sources and their positions on the image.

For the second step, we used the catalogue as input for PSFex, which computes a point-spread function (PSF) and constructs a small model image of the PSF, for each source in the catalogue. To avoid loss of information, PSFex uses a sampling that is different from the original pixels depending on the observations, a rougher sampling for the oversampled observations and a finer sampling forthe undersampled observations.

In the final step, we used the PSF file and the original images as input of a second SExtractor run. For this run, we used a significantly lower detection threshold, 1.5 times the estimated noise, to detect and later to remove even faint sources. The other parameters discussed above were not changed. The second run produces an image of objects which was subtracted from the original images resulting in an image in which stars appear as smooth holes.

A map of optical depth at 250 μm, derived from the Herschel submillimetre observations, and a map of optical depth of the J band, derived from the colour excess of background stars using the NICER method (Lombardi & Alves 2001), are shown in Fig. A.1. The maps are at 40′′ resolution. On the right is a masked map of the J-band surface brightness showing the area of the filament with τ(J) > 0.2. The images cover an area of 0.5° × 0.5°.

The observed intensities for the J, H, K, and 3.6 μm bands were acquired by taking average values over two small areas of the filament, one from the densest part of the filament and one from a part with less optical depth. The locations are marked in the middle panel of Fig. A.1 with white and blue plus signs, respectively. The obtained intensities are listed in Table A.1.

Fig. A.1 0.6° × 0.6° hydrogen column density map derived from the Herschel data (left), and the optical-depth map derived using the surface brightness of the scattered light (centre); right; a Ks band intensity map (MJy sr−1; Malinen et al. 2013). The positions used to represent the filament in the CRT and MCMC simulations are shown with a white plus sign (τJ ~ 6) and with a blue plus sign (τJ ~ 2).

The above method was also used to acquire values for the observed optical depths. We utilised two optical depth maps derived by Malinen et al. (2013). The first was τ₂₅₀ derived from the Herschel colour temperature and intensity maps with spectral index γ = 1.8. The second map, τ_J, was derived from the near-infrared colour excess following the NICER method described by Lombardi & Alves (2001). The obtained optical depths were τ_J = 6.10 and 2.05.

Appendix B Shape ofthe ISRF

In order to distinguish between the effects caused by the ISRF model and the cloud model, we compute a simulation with a

one-dimensional cloud model combined with a three-dimensional anisotropic ISRF model and assuming σ_NIR = 10%. The resulting marginalised probability distributions are shown in Fig. B.1.

Fig. B.1 Marginalised probability distributions of dust parameters for the case τJ = 6. For the radiative transfer computations, a one-dimensional cloud model and a three-dimensional radiation field is used. The white star indicates the projected position of the χ2 minimum and the black contour shows the 1σ of the projection. The colour scale shows the normalised probability.

Appendix C Additional MCMC simulations

The marginalised probability distributions for the TMC-1N τ_J = 6 position areshown in Figs. C.1 and C.2, where the albedo or the asymmetry parameter of the dust grains is increased (red lines) or decreased (purple lines) by 10% in the radiative transfer computations. For Fig. C.1 σ_NIR = 20% is assumed and for Fig. C.2

Fig. C.1 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. The grain albedo in the radiative transfer calculations is changed by −10% (red lines) or +10% (purple lines), compare to Fig. 15. For the surface brightness data σNIR = 20% is assumed. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.

Fig. C.2 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. The asymmetry parameter of the grains in the radiative transfer calculations is changed by −10% (red lines) or +10% (purple lines), compare to Fig. 15. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.

Fig. C.3 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 2 position. The grain albedo in the radiative transfer calculations is changed by −10% (red lines) or the asymmetry parameter has been changed by +10% (purple lines), compare to Fig. 15. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show the 1σ and 2σ.

σ_NIR = 10% is assumed. For both computations r_Si = 0.5 is used.

The marginalised probability distributions for the TMC-1N τ_J = 2 position is shown in Fig. C.3, where the asymmetry para- meter of the dust grains is changed by +10% (red lines) or −10% (purple lines) in the radiative transfer calculations. For both computations we use r_Si = 0.5, and σ_NIR = 10% is assumed.

References

All Tables

All Figures

	Fig. 1 Flow chart describing the steps from dust parameters to probability distributions.
In the text

	Fig. 2 Cumulative sums of absorption and scattering efficiencies for graphite (solid lines) and silicate (dotted lines) grains at J, K, and 3.6 μm bands. The red lines correspond to γ = 2.5, the blue lines to γ = 3.5, and the black lines to γ = 4.5.
In the text

	Fig. 3 Extinction efficiency Qext, albedo α, and asymmetry parameter g, for different dust mixtures. The optical properties are normalised so that for J band Qext = 1.
In the text

	Fig. 4 Grain properties as a function of the parameters Amax and γ of the grain size distribution for the J and K bands and for the ratio J∕K. The panels show the asymmetry parameter g; albedo α; and the scattering, absorption, and extinction efficiencies Qsca, Qabs, and Qext, respectively. The plot corresponds to rSi = 0.5, and the colour bars show the values of the dust parameters.
In the text

	Fig. 5 Results for one-dimensional radiative transfer models with rSi = 0.5, Amax = 1.0μm, and γ = 3.5 with averageτJ = 2 (left) and τJ = 6 (right). The panels show the optical depth of the J band and the surface brightness of the J, K, and 3.6 μm bands. The surface brightness values and optical depth have been normalised.
In the text

	Fig. 6 Left: intensity of the J band as a function of γ and Amax. Centre and right: the intensity ratios of the K to the J band, and of the 3.6 μm to the J band, respectively. A one-dimensional model is used with rSi = 0.5 and τJ = 6. The colour bars show the intensity and intensity ratios.
In the text

	Fig. 7 Spectral energy distributions resulting from a computation with rSi = 0.5. The maximum grain size is set to Amax = 1.0μm and γ = 3.5.
In the text

	Fig. 8 Surface brightness of the J and K bands, and the ratio of the K to the J band as a function of optical depth. The colours of the curves correspond to different values of Amax and γ. For the first column γ = 3.5, and for the second column Amax = 1.0μm.
In the text

	Fig. 9 Cross sections of the χ2 distribution through the location of the χ2 minimum for the synthetic observations. The white star indicates the χ2 minimum and the dashed lines show the dust parameter values used to simulate the synthetic observations. The colour bars correspond to the χ2 distribution.
In the text

Fig. 10 Marginalised one-dimensional probability distributions (panels A, F, K, and P) and two-dimensional projections of the probability of dust parameters for τJ = 6. The white star indicates the projected position of the χ2 minimum. The other symbols are the χ2 minima assuming a 10% higher or lower intensity for one of the bands, marked with green and red symbols, respectively. The circles correspond to J band, the diamonds to H band, the squares to K band, and the triangles to 3.6 μm band. The black contour shows the 1σ of the projection. The colour scale shows the normalised probability.

In the text

	Fig. 12 Schematic of the major axis and minor axis radiative transfer models (right) used to analyse the synthetic observations(left).
In the text

Fig. 11 Marginalised probability distributions of dust parameters for an ellipsoidal cloud model. In the radiative transfer modelling (see text) the cloud is viewed either along its minor axis (red lines) or its major axis (purple lines). For both cases, the optical depth along the line of sight is τJ = 6 and the dashed lines show the input values used in the computations. The stars indicate the projected positions of the χ2 minima and the black contours show the 1σ and 2σ of the projections. The greyscale map corresponds to the red lines.

In the text

	Fig. 13 As in Fig. 11, but the KISRF is set to a constant value of 1.2.
In the text

	Fig. 14 Marginalised probability distributions of dust parameters for τJ = 6. As in Fig. 8, but with the grain albedo in the radiative transfer calculations changed by −10% (red lines) or +10% (purple lines). The stars indicate the projected positions of the χ2 minima, and the greyscale map corresponds to the red lines. The dashed lines show the input values used in the computations.
In the text

	Fig. 15 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. Forthe surface brightness data we assume σNIR = 10% (red lines) and σNIR = 20% (purple lines). The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines, and the contours show 1 σ and 2 σ for the two models.
In the text

	Fig. 16 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. Forthe computations, only the J, H, and K bands were used. The white star indicates the projected position of the χ2 minimum and the black contours show the 1σ of the projection. The colour scale shows the normalised probability.
In the text

	Fig. 17 As in Fig. 15, but the grain albedo in the radiative transfer calculations is changed by −10% (red lines) or +10% (purple lines). The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.
In the text

	Fig. 18 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 2 position assuming σNIR = 10% (red lines) and σNIR = 20% (purple lines). The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.
In the text

	Fig. 19 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position with a rSi = 0.3 (red lines) and rSi = 0.7 (purple lines), respectively. For both computations σNIR = 10% is assumed. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.
In the text

Fig. 20 Near-infrared spectra derived from different dust models; a constant γ = 3.5 is used (left panel) and Amax = 1.0μm (right panel). The red and green lines are models where the albedo of the dust grains has been changed during the scattering computations by −10% (red) or by +10% (green), and the black line is a model where only the Amax and γ were varied during the scattering computations. Each symbol (circle, triangle, square) corresponds to a different Amax and γ value in the left and right panels, respectively.

In the text

	Fig. 21 Marginalised probability distributions of the dust parameters for the minor axis case (red lines) (see Fig. 11), and the case with grains smaller than 20 nm removed from the dust model used to fit the synthetic observations(purple lines).
In the text

Fig. 22 Massfraction of grains larger than 0.25 μm for τJ = 6 (red lines) and τJ = 2 (purple lines). For panels A–C, the Amax and γ values of the model with rSi = 0.5 and σNIR = 10% are used. Forpanel D, we use the previous model but assume σNIR = 20%. For panels Eand F, we use the models with the albedo of the grains decreased by 10% and the asymmetry parameter decreased by 10%, respectively. The greyscale colour maps correspond to the red lines and the contours show 1σ and 2σ.

In the text

	Fig. 23 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. The white star indicates the projected position of the χ2 minimum. Panels A, D, and G and C, F, and I show the marginalised probability integrated along the third parameter axis within ±5% of the χ2 minimum, respectively. Panels B, E, and H show the marginalised probability on the χ2 minimum.
In the text

Fig. A.1 0.6° × 0.6° hydrogen column density map derived from the Herschel data (left), and the optical-depth map derived using the surface brightness of the scattered light (centre); right; a Ks band intensity map (MJy sr−1; Malinen et al. 2013). The positions used to represent the filament in the CRT and MCMC simulations are shown with a white plus sign (τJ ~ 6) and with a blue plus sign (τJ ~ 2).

In the text

	Fig. B.1 Marginalised probability distributions of dust parameters for the case τJ = 6. For the radiative transfer computations, a one-dimensional cloud model and a three-dimensional radiation field is used. The white star indicates the projected position of the χ2 minimum and the black contour shows the 1σ of the projection. The colour scale shows the normalised probability.
In the text

Fig. C.1 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. The grain albedo in the radiative transfer calculations is changed by −10% (red lines) or +10% (purple lines), compare to Fig. 15. For the surface brightness data σNIR = 20% is assumed. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.

In the text

	Fig. C.2 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 6 position. The asymmetry parameter of the grains in the radiative transfer calculations is changed by −10% (red lines) or +10% (purple lines), compare to Fig. 15. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show 1σ and 2σ.
In the text

Fig. C.3 Marginalised probability distributions of dust parameters for the TMC-1N τJ = 2 position. The grain albedo in the radiative transfer calculations is changed by −10% (red lines) or the asymmetry parameter has been changed by +10% (purple lines), compare to Fig. 15. The stars indicate the projected positions of the χ2 minima. The greyscale map corresponds to the red lines and the contours show the 1σ and 2σ.

In the text