© ESO, 2013

1. Introduction

The structure of molecular clouds can be studied via a number of methods. These include molecular line mapping, observations of dust emission at far-infrared/sub-millimetre wavelengths, star counts in the optical and near-infrared (NIR) wavelengths, and measurements of colour excesses of background stars.

All techniques have their own drawbacks. For example, line and continuum emission maps are subject to abundance variations (gas and dust, respectively) and variations in the physical conditions, most notably the excitation and kinetic temperatures. Mass estimates based on dust emission can also be biased because of line-of-sight temperature variations, especially in high density clouds where star formation is potential (see, e.g., Malinen et al. 2011). The colour excess method provides column density estimates for extremely narrow lines of sight toward background stars. However, the intrinsic colours of the stars are usually unknown and this introduces significant noise, especially at low column densities. A full extinction map is obtained only after spatial averaging. This means that for all the listed methods the spatial resolution is usually some tens of arc seconds or worse. See, e.g., Juvela et al. (2006) for a more thorough review of the methods, Goodman et al. (2009) for a comparison of several methods, and Malinen et al. (2012) for a comparison of filament properties derived using NIR extinction and Herschel observations of dust emission.

Surface brightness caused by scattered NIR light provides another means of studying cloud structure. The first observation of scattered NIR light from molecular clouds illuminated by the normal interstellar radiation field (ISRF) was by Lehtinen & Mattila (1996). Later the observations of Nakajima et al. (2003) and Foster & Goodman (2006) (who also named the phenomenon “cloudshine”) have shown that it is now possible to obtain large maps of the surface brightness of normal interstellar clouds illuminated by the normal ISRF. NIR scattering can therefore be a new, complementary tool for studying the structure of dark clouds. See, e.g., Juvela et al. (2006) for a more complete review of the history of scattered light observations in dark clouds.

Padoan et al. (2006) presented a method to determine the cloud column density from the intensity of the near-infrared scattered light. Juvela et al. (2006) analysed the method in more detail using simulations, and developed a method to combine the surface brightness with extinction, to reduce errors caused by wrong assumptions of radiation field or dust properties. Starting with the known properties of the interstellar dust and ISRF, the papers made predictions for the visibility of the cloudshine and for the accuracy of the resulting column density estimates. They also demonstrated, independently of the direct evidence given by the Foster & Goodman (2006) data, that such observations are well within the capabilities of modern wide-field infrared cameras. The main advantage of the new method is the potentially extremely good spatial resolution. When the J, H, and K bands are used, the method remains accurate in regions with A_V even beyond 10 mag.

The same NIR observations provide data for the colour excess method, which means that the results of the methods can be compared at lower spatial resolution. Although both methods depend on NIR dust properties, the main sources of error are different. Comparison of the results can be used to study the values and variation of near-infrared dust properties (e.g., albedo and the shape of the extinction curve) and spatial variations in the strength of the local radiation field. Furthermore, correlations between wavelength bands give a direct way to estimate the point at which the contribution of dust emission from stochastically heated grains becomes significant. The level of the emission of these so-called Sellgren grains depends on the radiation field (Sellgren et al. 1996), but is still uncertain in normal interstellar clouds. On the other hand, the amount of the scattered light depends heavily on the grain size (Steinacker et al. 2010). These dependencies have implications for models of interstellar dust.

Juvela et al. (2008) used scattered NIR light to derive a column density map of a part of a filament in Corona Australis, and continued the analysis in a larger area in Juvela et al. (2009). They also compared the NIR data with Herschel sub-millimetre data in Juvela et al. (2012). Nakajima et al. (2008) applied a similar method to convert NIR scattered light to column density. They used the colour excess of individual background stars to calibrate an empirical relationship between surface brightness and column density, instead of an analytical formula, as used in Juvela et al. (2006, 2008). Scattered surface brightness from dense cores in the mid-infrared (MIR) was reported by Pagani et al. (2010) and Steinacker et al. (2010). Steinacker et al. (2010) named this phenomenon “coreshine” as a counterpart to “cloudshine”, which is also observed in the outer parts of the clouds.

In this paper, we study a filament in the Taurus molecular cloud using observations of NIR light of a 1° × 1° field observed with WFCAM instrument (Casali et al. 2007). Distance to the Taurus molecular cloud is ~140 pc, making it one of the closest relatively high latitude clouds, and consequently one of the most studied star-forming regions (see, e.g., Cambrésy 1999; Nutter et al. 2008; Kirk et al. 2013; Palmeirim et al. 2013). In Malinen et al. (2012), we compared the properties of this filament derived using NIR extinction and dust emission observed with Herschel. Here, we construct maps of the diffuse surface brightness, determine the intensity of the NIR scattered light, and derive the optical depth based on scattered NIR light using the method presented in Padoan et al. (2006). We compare the scattered light images with the other tracers, NIR extinction and sub-millimetre dust emission and, using these results, draw some conclusions regarding the intensity and spectrum of the local ISRF. We also perform radiative transfer modelling to compare observations with the level of NIR and MIR scattered light that is expected using standard ISRF levels and standard dust models.

The contents of the article are the following: We present observations and data processing in Sect. 2. We describe the method for deriving optical depth from scattered NIR surface brightness in Sect. 3. We derive NIR surface brightness maps and optical depth maps based on observations of dust emission, NIR extinction, and NIR scattered surface brightness and compare the results in Sect. 4. We describe radiative transfer modelling of a filament seen in scattered light in Sect. 5. We discuss the results in Sect. 6 and present our conclusions in Sect. 7.

2. Observations and data processing

2.1. WFCAM

We have used the Wide Field CAMera (WFCAM) (Casali et al. 2007) of the United Kingdom InfraRed Telescope (UKIRT) to observe a 1° × 1° field in the NIR J, H, and K bands (1.25, 1.65, and 2.22 μm, respectively). The field, which we call TMC-1N (Malinen et al. 2012), is in the Taurus molecular cloud complex north of TMC-1. The central coordinates of this field are RA (J2000) 4h39m36s and Dec (J2000) +263932. At a distance of 140 pc, this corresponds to an area of ~(2.44 pc). The Galactic latitude of this area is approximately −13.3°. The target field was chosen based on the Taurus extinction maps of Cambrésy (1999) and Padoan et al. (2002). According to Rebull et al. (2011), there are no strong young stellar objects (YSOs) in TMC-1N that could cause additional scattered light and therefore complicate the analysis. There is mainly just one continuous filament in the field. The A_V range of the area is suitable for our method: there are some regions with A_V ~ 20^m, but the mean value remains well below 10 mag, even in most parts of the filament.

We applied sky correction using offset fields to be able to measure faint surface brightness features. The observations were made during 13 nights between 2006–2008 using 2  × 2 pointings towards the selected field. Because WFCAM consists of four separate CCD arrays, the result is a 1° × 1° image consisting of 4 × 4 subimages. We used four separate OFF fields. The standard data reduction was conducted in accordance with the normal pipeline routine¹, including, e.g., dark-correction, flatfielding (including internal gain correction), decurtaining, sky correction, and cross-talk. Standard decurtaining methods² were needed to remove the stripes caused by the instrument. The details of the observations are shown in Malinen et al. (2012) where we compared the filament properties derived using NIR extinction and sub-millimetre dust emission observed with Herschel. There, the NIR data were calibrated to magnitudes with the help of 2MASS catalogue stars to derive an extinction map with the NICER method (Lombardi & Alves 2001). Here, in order to study the surface brightness, we calibrated the data from magnitudes to MJy/sr units using aperture photometry of several stars in each frame.

The reduced images contained residual gradients in each sub-image, shown as brightening of the signal towards the frame edges. These are probably of instrumental origin. The size of the gradients as a percentage of overall sky level varies, but is typically between 0.002–0.008. Dye et al. (2006) and Warren et al. (2007) report several types of artefacts in WFCAM data and note that removing sky subtraction residuals is a complex problem, especially when observing near the Moon. However, during our observations the Moon was always further than 40 degrees away. In addition, the moonlight artefacts that cause most problems are local scattered light from dust on the optics, not gradients. The master twilight flats may have low level gradients present in them, given the size of the WFCAM field-of-view. The dark correction can also leave low level reset anomalies, particularly near the detector edges. It is possible that our data are showing comparatively low-level residual issues that the general surveys have not encountered.

We first removed the stars using the Iraf Daophot-package, masked the remaining bad pixels (noisy borders and residuals of bright or saturated stars) and performed the surface brightness calibration. We modelled the gradients with a method described in the following subsection (Sect. 2.1.1). We convolved the calibrated frames with a 2 Gaussian beam and resampled the data onto 0.8 pixels. The obtained gradients were then subtracted from the frames. The maximum intensity of the filament is a few times higher than the typical magnitude of the instrumental gradients.

We further median-filtered the images to ~16.8 resolution to diminish the effect of residual stars and stellar artefacts, and finally convolved them to ~40 resolution, for later comparison with Herschel observations at that resolution. We combined the frames to a full map using Montage³ and resampled the data onto 8 pixels for the analysis. For comparisons with NICER data, we further convolved the intensity maps to 60 resolution. We subtracted the background from the maps, using an area of low column density outside the filament as a reference, see Fig. 9 (middle frame).

A map of visual extinction, A_V, of the TMC-1N area was already derived and presented in Malinen et al. (2012), using a Cardelli et al. (1989) extinction curve with R_V = 4.0. In this paper, we mainly use the J band optical depth τ_J, which is in practice independent of the assumption of the R_V value. For comparison, a Cardelli et al. (1989) extinction curve with R_V = 4.0 gives the relation τ_J = 0.2844 × A_V.

Roy et al. (2013) discuss the effect of the finite width of the filters on the observed extinction. In their Fig. 10, they show the relation of colour excess E(J − K_S) obtained with 2MASS filters to the corresponding monochromatic colour excess as a function of N_H. In TMC-1N, having an A_V of mostly less than 20, the effect is small (less than 5% even in the densest parts), and we have not made this correction.

2.1.1. Gradient modelling

During the data reduction, separate analysis was carried out to model and subtract the residual artefacts of the surface brightness data. To be able to model the large scale gradients, we first filtered the masked and calibrated data with a median filter to  ~40 resolution. The gradients were strongest near the borders, which meant that the overlap between the frames could not be reliably used to aid the gradient removal. To model the gradient, therefore, we assumed a similar relation between the column density and the surface brightness as the one described in Sect. 3. We used the Herschel column density map described in Sect. 2.2 as an independent tracer of the filament structure in the fitting of the gradients. In principle, the extinction map obtained from NIR reddening of background stars could also have been used, but we chose Herschel data because of its lower noise.

We modelled the gradients with a third order surface describing the artificial gradient in the image, with an additional column density dependent term. We used least squares fitting to minimize, separately for each band i in (J, H, K) and each frame, the residual (1)where P is the complete third order polynomial for a 2-dimensional surface (10 terms), I_i is the observed intensity map in band i, and N the column density derived from Herschel data.

The parameter b describes the saturation of the relation between I and N, and it mainly depends on the grain properties (Juvela et al. 2006). If we presume that the dust grains in the cloud are similar as in the model of Juvela et al. (2006), the same b values should apply as long as the optical depth is not high. Therefore, we used the parameter values b_J = 0.34, b_H = 0.23, and b_K = 0.15, (in 1/mag units, as the model used A_V values instead of N) taken from the model data. The parameter b does not depend on the strength of the incoming radiation, but in some amount on the direction distribution. The model was based on a large cloud with an inhomogeneous density distribution and an isotropic radiation field. TMC-1N, on the other hand, has a single, filamentary structure, and a larger part of the radiation comes in from the observer’s side. We do not expect the model to describe the filament perfectly, but the model parameter values served as a good starting point.

The coefficients of the polynomial are free parameters and are used to describe the artificial gradients that will be removed in the subsequent analysis. The parameters a_i are related to the actual signal from the cloud and are the same for all frames in the same band. They are also kept as free parameters to minimise the effect the cloud structure has on the gradient fits. The exponential term M = a_i(1−e^−b_iN) is only an approximation of the column density dependence of the scattered light. However, because gradients are described only as third order surfaces over each frame and scattered light is visible only in some parts of each frame, the obtained gradient model is not expected to be strongly dependent on the reference map, here the column density derived from Herschel. In particular, the scaling between column density and the surface brightness is a free parameter and thus no assumption is made of the expected level of NIR surface brightness. Furthermore, the non-linear approximation represented by the term M becomes important only at A_V ~ 10^m. To obtain a better fit for the central part of the frames used in the following analysis, the gradients were fitted excluding a 10% wide border in each frame. For comparison, we also performed a fit to the whole area of the frames, in order to better fit the borders of the frames. The fitting provides the model P (for each frame and band) of the gradients that we subtract from the original surface brightness images. Note that the exponential term M is not subtracted as that is not part of the artefact.

We tested the reliability of the gradient removal method in several ways. We calculated the residual between the used model and the corrected surface brightness map, S = I − P, with the function Res = a_i(1−e^−b_iN) − S_i, where b_i are the constants used in the gradient modelling and a_i are the values obtained from the fit. N is the map used in the gradient modelling. We used surface brightness maps in ~2.2″ resolution. The mean values for the residual are –0.0052, 0.00035, and 0.000074 MJy/sr for the J, H, and K bands, respectively. The standard deviation for the residuals of the same bands are 0.017, 0.031, and 0.021 MJy/sr. Noise, especially near the borders, increases the obtained standard deviation values. The absolute values for the mean residuals for all bands are less than 0.006 MJy/sr, indicating that the minimisation of the residual in Eq. (1) has been effective. The filament is not apparent in the residuals nor in the removed gradients.

We also compared the surface brightness maps derived using different areas for the fitting of the gradient, full area or excluding 10% wide borders of each frame. After the gradient removal, we compared the surface brightness maps (in 40″ resolution) by calculating the difference between the maps in the final masked area shown in Fig. 2. The means of the difference of the final masked area are 0.0026, 0.0026, and –0.0014 MJy/sr for the J, H, and K bands, respectively. The standard deviations for the difference of the same bands are 0.0044, 0.0072, 0.0052 MJy/sr. This indicates that small differences in the fitting area do not cause a significant difference in the central parts of the corrected frames. Near the borders of the frames, the difference can be larger.

2.2. Herschel

The Taurus molecular cloud has been mapped with Herschel (Pilbratt et al. 2010) as part of the Gould Belt Survey⁴ (André et al. 2010), see Kirk et al. (2013) and  Palmeirim et al. (2013). We used SPIRE (Griffin et al. 2010) 250 μm, 350 μm, and 500 μm maps of the TMC-1N field. We obtained the data from the Herschel Gould Belt Survey consortium. The observation identifiers of the data are 1342202252 and 1342202253. The details of the Herschel data and the derivation of column density map are presented in Malinen et al. (2012).

For easier comparison with other maps and to avoid the assumptions needed to make column density maps, we used optical depth τ₂₅₀ instead. We use the notation τ₂₅₀ to mean optical depth τ_ν at wavelength 250 μm (ν = 1200 GHz). We derived τ₂₅₀ maps from Herschel intensity and colour temperature maps, using a value of 1.8 for the spectral index β. Recent studies with Planck (Planck Collaboration 2011) show that this is a good estimate for Taurus. For comparison, we also made the analysis with the commonly used value of β = 2.0.

2.3. Spitzer

The Taurus filament is covered by observations of the Spitzer IRAC instrument (Fazio et al. 2004). The IRAC data (observation numbers 11230976 and 11234816) are from the Taurus Spitzer legacy project (PI D. Padgett). The archival IRAC maps show some checkered pattern. We therefore started with the pipeline-reduced, artefact-mitigated images (cbsd images). The residual offsets were corrected by subtracting the median value from each frame and by estimating the residual offsets with respect to the median of all frames⁵. The final mosaic image was then made with the MOPEX tool (Makovoz et al. 2006).

3. Deriving optical depth from scattered NIR surface brightness

Padoan et al. (2006) presented a method for converting scattered surface brightness into column density. Juvela et al. (2006, 2008) discussed the method and its reliability in more detail. We review the main points of the method given in these articles and formulate it to suit this study.

Based on the one-dimensional radiative transfer equation in a homogeneous medium, the relationship between surface brightness (I) and column density (N) can be approximated with the formula (2)where the parameters a and b are positive constants defined separately for each band. The parameter b scales the column density to optical depth, τ. The function depicts the radiation coming from the cloud. We presume that the observed radiation is mainly caused by scattering. Also other tracers of cloud structure, such as extinction A_V or optical depth τ can be used instead of column density N. In that case, the value and unit of the parameter b must be changed accordingly. The parameters b mainly depend on the NIR dust properties and the parameters a on the radiation field, although the interpretation is not this straightforward (Juvela et al. 2006). Equation (2) is exactly valid only for a homogeneous cloud. In real clouds it works only as an approximation. It is expected that both parameters change if the dust properties, radiation field or cloud geometry change. Both can, however, be treated just as empirical parameters.

In areas of low optical depth, the intensity of NIR scattered light is directly proportional to the column density as shown in the relation (3)The product ab gives the scattered intensity per column density. With higher column densities (A_V ~ 10^m), the NIR intensity values begin to saturate, starting from the shorter wavelengths. Thus, the relation becomes more non-linear with increasing column density.

Equation (2) can be considered as an empirical description of the relation between surface brightness and column density, and it is not necessarily the optimal functional form to describe this relation. However, as long as the data points follow this relation, it is possible to derive approximations for column density. Based on numerical simulations (Padoan et al. 2006; Juvela et al. 2006), this function is rather reliable in areas of low to medium visual extinction (A_V = 1−15^m), where the saturation is not strong. If no independent column density estimates are available, only the ratios between different bands can be derived from surface brightness observations. However, NIR observations also provide the colour excesses of background stars, and an extinction map can therefore be calculated and used to derive the necessary parameter values (Juvela et al. 2006).

To derive an optical depth map from the surface brightness maps, we obtain a and b parameter values for each band from correlations by fitting Eq. (2) to the data, using optical depth instead of column density N in the exponential part. Thus, Eq. (2) is used to establish a non-linear scaling between the extinction that increases linearly with column density and the surface brightness for which the increase is non-linear. Once the parameters of this mapping have been determined, we can convert surface brightness observations directly to estimates of column density, possibly at a resolution higher than what is possible using background stars alone. Using the obtained parameter values, we calculate a value for each pixel of the map by minimising the squared sum of residuals (4)where I_i are the pixel values of the intensity maps, and (5)for each i in (J, H, K). Minimisation is needed, as the observed intensities do not follow Eq. (2) exactly due to noise and possible model errors.

4. Results

Surface brightness map in J band with 40 resolution is shown in Fig. 1. In this figure, the gradient removal is made using a fit to the whole area of each frame. In the analysis, bad pixels, border regions suffering from imperfect gradient removal, and the main artefacts remaining after star removal were all masked. We also limited the analysis to areas where Herschel column density values are above 1.0 × 10²¹cm. The masked J map used in the analysis is shown in Fig. 2. We used the same mask for all maps (J, H, K, NICER, and Herschel) in the analysis. The reference area used for background subtraction is shown in Fig. 9 (middle frame). The median τ_J in this area is ~0.1.

A combined three-colour image of WFCAM J, H, and K band intensity maps with 40 resolution is shown in Fig. 3. Here, the removed gradients were fitted to the frames excluding 10% wide borders of each frame. The remaining gradients can most clearly be seen near the corners of the frames. Borders of the frames have been masked.

4.1. Correlations between NIR surface brightness and optical depth derived from extinction and dust emission

We show correlations between the NIR surface brightness in J, H, and K bands for the main filament area in Fig. 4. The main cloud of points is concentrated to values below 0.08 MJy/sr in I_J, 0.17 MJy/sr in I_H, and 0.12 MJy/sr in I_K. The correlation between I_J and I_K is approximately linear up to I_K ~ 0.06 MJy/sr, after which I_J starts to saturate. Similarly, the correlation between I_H and I_K is approximately linear up to I_K ~ 0.08 MJy/sr, after which I_H starts to saturate.

Compared to Fig. 3 in Juvela et al. (2008) for Corona Australis, our WFCAM data for TMC-1N show similar linear relations between I_J and I_K and between I_H and I_K in the low end of the intensity scale, although our data show slightly lower values. In TMC-1N, the main cloud of points is concentrated to I_K values below 0.12 MJy/sr, whereas in Corona Australis there are also points up to I_K ~ 0.7 MJy/sr. In Corona Australis, the correlations turn to negative between I_K ~ 0.4–0.7 MJy/sr.

Correlations between surface brightness, I_ν, and optical depth in J band derived from extinction, , are shown for bands J, H, and K in Fig. 5. We fit Eq. (2) to the data, using optical depth instead of column density N in the exponential part. As a result, we obtain for each band the a and b parameters needed for the derivation of the optical depth map. The fitted values are a_J = 0.08, b_J = 0.90, a_H = 0.20, b_H = 0.44, a_K = 0.19, and b_K = 0.21. For comparison, we plot in the same figures the data of Corona Australis from Juvela et al. (2008). We have scaled those data to I_ν units (MJy/sr), as a function of τ_J. Compared to Corona Australis, we obtain notably lower values for I_ν as a function of τ_J in all three bands. We also fit Eq. (2) to these data and show the obtained parameters in the figure. For comparison, we also show correlations between I_ν and optical depth derived from Herschel maps, τ₂₅₀, in Fig. 6.

We made error estimates for the a and b parameters obtained from the fitting of Eq. (2) using two different methods. First, we used a standard bootstrap method to estimate the errors caused by sampling. We made the fit using 100 different samples from the data, all samples having the size of the full dataset. The standard deviations for the a and b parameters and their product a × b, that is the slope of the linear part of the function, are shown in Table 1. The standard deviation for parameter a can be up to 0.004, and for parameter b up to ~0.006, but for the slope the uncertainty is rather small,  ~0.0002 for K, ~0.0003 for J, and ~0.0004 for H band.

Secondly, we tested the possible systematic effect of changing the upper limit of the fitted range between τ_J values 1.5–6.0. In Table 2 we show the relative change in the parameter values when the upper limit of the fitted range is changed from τ_J = 2.0 to τ_J = 6.0, as between these values the change is mainly systematic. Even though the change in the a parameter can be up to –33% and in the b parameter up to 61%, the change in the slope of the linear part is less than 10% for each band. Below a τ_J value of 2.0, the parameter values can change in a more unpredictable way.

We derive emissivity, the ratio of FIR dust emission to column density, by using τ_J obtained from NIR extinction map as an independent tracer of column density. Correlation between τ₂₅₀ and with β = 1.8 is shown in Fig. 7 for the main filament area shown in Fig. 2. The slope of a straight line fitted to the range –4 is  ~0.0013. We also tested the possible change of slope in the low and middle range. The derived slope is 0.0013 for range 0–2, and 0.0015 (more precisely 0.00149) for range 2–4. If we change β from 1.8 to 2.0, the slope between 0–4 increases ~32% to the value ~0.0018. We also made a similar fit using all the data in the maps. The derived slope is 0.0013 for range 0–2, and 0.0014 (more precisely 0.00144) for range 2–4, indicating that there is no significant change when compared to the masked area.

For comparison with other studies, we convert our slope of τ₂₅₀/τ_J = 0.0013 to opacity or dust emission cross-section per H nucleon (6)where τ_ν is optical depth, N_H is the total H column density (H in any form), μ is the mean molecular weight per H (1.4), m_H is the mass of H atom, and κ_ν is the mass absorption (or emission) coefficient (cm/g) relative to gas mass, also often called opacity. We again use wavelength instead of frequency in our notation: σ_e(250) = σ_e(250 μm) = σ_e(1200 GHz).

We derive the τ₂₅₀ map directly from the Herschel observations, and use our WFCAM NIR extinction map as an independent tracer of the column density. We use the conversion factor for diffuse clouds N(HI + H₂)/E(B − V) = 5.8 × 10²¹ cm/mag of Bohlin et al. (1978). We derive the relation E(B − V)/E(J − K) from Cardelli et al. (1989) extinction curves, and obtain values 1.999 (with R_V = 3.1) and 1.413 (with R_V = 4.0), leading to relations N(H) = 11.59 × 10²¹E(J − K) (with R_V = 3.1) and N(H) = 8.196 × 10²¹E(J − K) (with R_V = 4.0). Martin et al. (2012) have observed a similar relation, N(H) ~ 11.5 × 10²¹E(J − K_S), for regions of moderate extinction in Vela. Cardelli et al. (1989) extinction curves also give the relation A_J/E(J − K) ~ 1.675 (with both R_V values 3.1 and 4.0). Converting magnitudes to optical depths (A = 2.5lg(e)τ ~ 1.086τ) gives the relation E(J − K) ~ τ_J/1.54. This leads to the relations σ_e(250) = 1.33 × 10^-22τ₂₅₀/τ_J cm/H (with R_V = 3.1) and σ_e(250) = 1.88 × 10^-22τ₂₅₀/τ_J cm/H (with R_V = 4.0).

The value τ₂₅₀/τ_J = 0.0013 leads to values σ_e(250) = 1.7 × 10^-25 cm²/H (with R_V = 3.1) or σ_e(250) = 2.4 × 10^-25 cm²/H (with R_V = 4.0). These can be converted to κ_ν values 0.07 cm/g or 0.10 cm/g, respectively.

The map of the ratio is shown in Fig. 8. The map shows no clear evidence for systematic increase of the ratio τ₂₅₀/τ_J with increasing density inside the densest filament. Some high value areas can be attributed to imperfections in the two maps, τ₂₅₀ and τ_J. For instance, the high value area next to the densest part of the filament seems to be caused by the different shape of the filament in these two maps. There the value of the ratio is high, because the filament is slightly more narrow in τ_J, possibly caused by the lack of background stars seen behind the densest filament.

4.2. Optical depth derived from scattered light

As described in Sect. 3, in order to derive an optical depth map from the surface brightness maps, we obtain a and b parameter values for each band from correlations between I_ν and shown in Fig. 5. We calculate a value for each pixel of the map by minimising Eq. (4).

We use 40 and 60 maps in the following analysis, but also show a higher resolution (~2.2) map in Fig. 9 (right frame). In the same figure, we also show maps of Herschel τ₂₅₀ and NICER . In Fig. 10, we show close-ups of the same maps.

Correlations between (resolution 40) and (resolution 60) are shown in Fig. 11. The reference area used for background subtraction is shown in Fig. 9 (middle frame). The fitted values for the slopes are 808 for and 1.019 ~  1 for . This is as expected, since was derived based on the correlation between I_ν and . The correlations are linear up to . Above that, the values saturate strongly. The areas where are marked with contours in Fig. 9 (middle frame), indicating that these form only a small area in the densest clumps inside the cloud.

4.3. Filament cross-sections

Median profiles of the filament are shown in Fig. 12, for bands J, H, and K (upper frame) and for and (lower frame). The cross-sections are taken from the area between the two red lines shown in Fig. 9 (middle frame), where we have continuous data in all WFCAM maps. This is not the densest part of the filament, but a moderately dense area next to it. The cross-sections show that the derived gives rather similar results to the map for the filament profile, except for an extra peak in the profile. The profiles seen in J and K bands are rather similar, whereas the H band profile has approximately two times stronger peak than the other two bands. The filament width or FWHM is ~3′ ~0.1 pc, as shown in Malinen et al. (2012).

4.4. Spitzer data

We compare the IRAC 3.6 μm surface brightness with the Herschel optical depth τ₂₅₀, to study also the correlations between MIR scattered light and sub-millimetre dust emission. The surface brightness at  ~3.6 μm can still be dominated by light scattering, while at longer wavelengths, the scattering decreases and the signal is expected to be dominated by dust emission. Some contribution of dust emission cannot be excluded even around 3.6 μm and we conservatively consider the Spitzer data only as an upper limit on the intensity of the scattered light.

The point sources are a major problem in estimating the level of the extended emission. Because the examined area is small, we can produce an extended MIR surface brightness map by manually masking all the visible point sources. The masks altogether cover 24 arcmin. After the masking, we carry out median filtering. The filter calculates the 25% percentile of all unmasked pixels within a given radius that is fixed to 5. If the area contains less than ten unmasked pixels, the value is left undefined. The images are then convolved with a Gaussian beam to produce final surface brightness maps at the resolution of 40 that corresponds to the resolution of the Herschel data. The convolution ignores the undefined values.

To reduce the noise (and because of the uncertainty of the fidelity of the large scale flat-fielding), data are correlated only around the main column density peak. The column density threshold of 4 × 10²¹ cm defines an area of ~77 arcmin.

The Spitzer 3.6 μm map and the area used in the correlations are shown in Fig. 13. The resulting correlations between MIR surface brightness and optical depth τ₂₅₀ are shown in Fig. 14. In the figure, we also show the median surface brightness that is calculated for τ₂₅₀ bins with a width of 0.002. We fit a robust least squares line to the individual surface brightness values in pixels where τ₂₅₀ is between 0.0025–0.0065. The fit is performed iteratively, discarding points falling further than 2.5σ from the fitted line. The value of the fitted slope is 2.40 MJy/sr.

4.5. Spectral energy distributions

Spectral energy distributions are shown in Fig. 15. We calculate the values for I_ν/τ_J in TMC-1N with Eq. (3), using τ_J instead of column density N. In areas of low optical depth, the product ab gives the ratio I_ν/τ_J. The obtained values for J, H and K bands are 0.074, 0.086, and 0.038 MJy/sr, respectively.

We compare our results with the values obtained for Corona Australis in Juvela et al. (2008) (filled squares in their Fig. 13). We have scaled the Corona Australis values to I_ν/τ_J units (with R_V = 3.1). We have also marked Mathis ISRF model (Mathis et al. 1983) values (with R_V = 4.0) on the figure for comparison. The values are obtained by multiplying the Mathis intensities with the scattering cross-sections of the Draine (2003) dust model.

The values of J, H and K bands in TMC-1N are approximately one third of the values in Corona Australis. Mathis model values for H and K are ~1.5 times as high as in TMC-1N. For J band, the Mathis model value is already ~3 times as high as the TMC-1N values. The value we obtain for Spitzer 3.6 μm is also only approximately one third of the value given by the Mathis model. In both TMC-1N and Corona Australis, I_ν/τ_J is lower in the J band than in the H band. However, in the Mathis model, the value for the J band is notably higher than for the H band.

5. Radiative transfer modelling

We carry out radiative transfer calculations of the ISRF light that is scattered from the cloud in the wavelength range 1.2–3.5 μm. We construct a realistic three-dimensional model of the density distribution of the densest part of the filament. Based on the density structure, we calculate predictions for the surface brightness using the Diffuse Infrared Background Experiment (DIRBE) observations on the Cosmic Background Explorer (COBE) as a template of the sky brightness. The radiative transfer calculations were carried out with a Monte Carlo radiative transfer program (Juvela & Padoan 2003).

The filament is illuminated by an anisotropic radiation field that, in conjunction with the scattering phase function, affects the strength and spatial distribution of the scattered intensity. We define the intensity of the incoming radiation using the DIRBE all-sky maps⁶ (Zodi-Subtracted Mission Average (ZSMA) Maps). The DIRBE bands 1, 2, and 3 can be used directly to specify the intensity at the corresponding wavelengths 1.2 μm, 2.2 μm, and 3.5 μm. The H band intensity is obtained by multiplying the average of the J and K band intensities with the intensity ratio I_H/ ⟨ I_J,I_K ⟩  taken from the Mathis et al. (1983) ISRF model. Note that, averaged over the whole sky, the total intensities themselves are ~50% higher than in the Mathis et al. (1983) model.

The density distribution of the model cloud is based on the column density map that is derived from the Herschel sub-millimetre observations (see Malinen et al. 2012). This is based on the colour temperature of the emission and on the assumption that dust opacity follows the law 0.1 cm/g (ν/1000 GHz), with β = 2.0. The actual cloud model corresponds to 17′ × 17′ area that we cover with 68 × 68 pixels, 15″ in size. Our three-dimensional cloud model is correspondingly a cube that consists of 68³ cells and is viewed along one of its major axes.

The Herschel column density map constrains the mass distribution only in the plane of the sky. To construct a three-dimensional filament, we start with a cylindrical structure where the radial density profile follows the Plummer function (7)with parameters ρ_C = 4 × 10⁴ cm, R_flat = 0.03, and p = 3.0. These are close to the values previously determined from the fitting of the Herschel observations (Malinen et al. 2012). This initial cylindrical structure therefore has properties close to those of the average filament. We modify this initial model by requiring that, for each line-of-sight, the column density of the model exactly matches that of the Herschel column density map. We calculate for each pixel the ratio of the Herschel column density and the initial column density in the model. The densities along the same line-of-sight are then multiplied by this value. The procedure modifies the filament so that it is no longer cylinder symmetric. However, the deformations remain relatively small so that the ratio of the major and minor axis of the 2D filament cross sections always remains below two.

The absolute value of the background affects the contrast between the cloud and the background. The surface brightness relative to the brightness of the background sky is calculated as (8)where I is the total observed signal, I_sca the scattered signal coming from the cloud, I_bge ^− τ the background signal coming directly (without scattering) through the cloud with optical depth τ, and I_bg the background signal. The radiative transfer program takes into account the incoming ISRF and the absorption and scattering of those photons inside the cloud. I_sca is the amount of light coming out of the cloud after this process.

The term I_bg was estimated using WISE (Wright et al. 2010) and DIRBE data. The WISE map of a one degree diameter area around the column density peak was first adjusted to the absolute level indicated by the DIRBE data. We read the I_bg value from the corrected WISE map, using the area that corresponds to the area where the model column density is approximately zero (corresponding to the upper part of Fig. 16). To exclude the effect of point sources in this area, we use the 40% percentile of pixels in this reference area. The values obtained for the background I_bg are 0.122 MJy/sr (J), 0.0878 MJy/sr (K), and 0.0765 MJy/sr (3.4 μm). We estimate a value 0.10 MJy/sr for the H band. The cosmic infrared background (CIB) between the J and 3.5 μm bands is notably lower than these values,  ~0.01–0.02 MJy/sr (see, e.g., Cambrésy et al. 2001; Wright & Reese 2000).

The dispersion of DIRBE values over the area used for comparison with WFCAM and WISE is  ~10% of the estimated I_bg value. Because this includes not only noise but also real surface brightness variations, the statistical error of the mean DIRBE value is significantly lower. However, the dispersion is also a measure of the possible difference between the background value I_bg derived for the reference area and the actual background at the position of the filament. The final I_bg is calculated as 40% percentile of the reference area within WFCAM and WISE maps. The standard deviations of all pixels falling below the 40% value (assumed not to be significantly contaminated by point sources) are 0.011, 0.015, and 0.005 MJy sr for J, K, and 3.4 μm bands, respectively. The zero point uncertainty of the zodiacal light model subtracted in DIRBE ZSMA maps is ~0.006 MJy sr in the J band and less at longer wavelengths (Kelsall et al. 1998), but differences between zodiacal light models can be even larger (Wright 2001). Furthermore, Taurus is located at low ecliptic latitude and two thirds of the total signal observed by DIRBE consists of zodiacal light. A relative error of 10% in the zodiacal light model would therefore amount to ~0.04, 0.03, and 0.02 MJy sr for J, K, and 3.5 μm, respectively, and could be the dominant source of error.

Diffuse emission that originates between the filament and observer does not contribute to Eq. (8), which describes the surface brightness contrast between the source and the background. In that equation, I_bg corresponds to that part of the diffuse component that truly resides behind the filament. We can assume that most of the diffuse material is associated with the Taurus cloud but can reside either in front of or behind the filament. We therefore make the assumption that I_bg corresponds to half of the values derived above.

If the diffuse component is closely associated with the filament, consistency requires that its effect is also taken into account when calculating the component I_sca. This term corresponds to photons scattered from the filament, excluding photons that may be scattered in some envelope around it. In this case, the filament (i.e., the densest part that is being modelled) is not illuminated by the full ISRF but by an ISRF that is attenuated by e^−0.5τ_d, where τ_d is the optical depth of the envelope. Similarly, when a photon scatters within the model, it must travel through a similar layer a second time (half of the full τ_d) before reaching the observer. Thus, the observed signal is reduced by a factor e^−τ_d.

We estimate that for the diffuse cloud , and calculate values 0.34, 0.2, and 0.0998 for , , and , respectively, using Cardelli et al. (1989) extinction curves. We compare three cases: (1) no background I_bg; (2) with background I_bg; and (3) with background I_bg and attenuation of the scattered light in the envelope. We show the corrections in function (9)where C is either 0 (case 1) or 0.5 (cases 2 and 3), indicating the amount of background, and use value 0 for τ_d for cases 1 and 2, and the τ_d values derived above for case 3.

In the calculations, we use the NIR dust properties of normal Milky Way dust (Draine 2003) and the tabulated scattering phase functions that are publicly available⁷. We do not seek a perfect match to the observations but explicitly assume the column density distribution and dust properties as described above. The comparison with the observations is thus a test for the consistency of these assumptions.

The τ_J map of the model is shown in Fig. 16 and simulated scattered surface brightness maps in Fig. 17. Correlations of I_ν and τ_J from simulations are shown in Fig. 18 for intensity in bands J, H, and K, using the three cases for the background correction described above. For comparison, we plot in the same figures the data from our TMC-1N observations. We have subtracted the background from both the simulated and observed data using as reference area the area in which the model column density is lowest (corresponding to the upper part of Fig. 16). In the WFCAM data, the median τ_J in this area is  ~0.5. Note the difference from Fig. 5, in which we used as reference area the low optical depth area marked in Fig. 9. We fit Eq. (2) to both of the data, using optical depth instead of column density N in the exponential part. As a result, we obtain the a and b parameters for each band.

In the H band, the simulations give ~1.5 times as high intensity values as the observations, without background (case 1). The observed values are mostly settled between the cases 2 and 3 in the simulations. In the K band, the simulations give nearly two times as high values to the observations, when no background is included (case 1). In case 3, with the background and attenuation in the envelope, the simulations give similar values as the observations. However, in the J band, the simulated values are approximately three times as high as in our observations, when no background is included (case 1). Even with background (case 2), the simulations still give approximately two times as high values. Only in case 3, with the attenuation in the envelope, the simulations give similar values to the observations. Even when the maximum intensities for observations and simulations are approximately the same, the form of the fitted Eq. (2) can be different, leading to notably different values for the parameters.

Correlation between I_ν and τ_J for the 3.5 μm band is shown in Fig. 19, similarly using the three cases. At this wavelength, case 1 still gives two times as high values as case 2, but the difference between cases 2 and 3 is only ~10%.

6. Discussion

We have studied a filament in the Taurus molecular cloud using NIR images of scattered light. The observations carried out with WFCAM instrument cover an area of 1° × 1° corresponding to ~(2.44 pc) at the distance of 140 pc, making this, to our knowledge, the largest NIR map where the surface brightness is analysed in detail. We have analysed the faint surface brightness, determined its intensity (attributed to light scattering), and used it to derive an optical depth map based on the method of Padoan et al. (2006). We have compared the data derived from NIR scattering, NIR extinction and Herschel dust emission.

The signal-to-noise value, S/N, of our WFCAM data is lower than expected. There were significant artefacts, presumably of instrumental origin, namely the curtain effect and large scale gradients. It is possible that the excess noise is caused by these effects and the processes used to remove them. Due to this, the data could not be used at as high a resolution as expected. In principle, scattered light can be measured at arcsecond resolution. However, as shown here, it can be difficult to obtain reliable maps at the full resolution.

Correlations between I_ν and in Fig. 5 suggest that the radiation field in TMC-1N is notably lower in all three NIR bands when compared to similar observations of Corona Australis (Juvela et al. 2008). This reaffirms the results obtained with dust emission models of Juvela et al. (2012), suggesting that the radiation field in Corona Australis is at least three times that of the normal ISRF.

NIR intensity as a function of column density followed roughly the functional form assumed in Eq. (2). We have estimated the possible errors in obtaining the parameters a and b for each band from fitting Eq. (2) to the correlation between I_ν and . Errors caused by sampling are small, but changes to the fitting limit can cause up to ~60% deviations in a single parameter. However, the change in the product a × b, that is the slope of the linear part of the Eq. (2), is only up to ~10% for each band. As discussed in Sect. 3, Eq. (2) and the parameters a and b can be thought of as an empirical model that represents the relation between surface brightness and column density. The model may not be optimal and no direct physical interpretation can necessarily be attached to its parameters. Bias in the fitted parameter values could cause bias also when minimising the Eq. (4) to obtain an optical depth map.

We obtained the value 0.0013 for the slope with β = 1.8, using values in the range τ_J = 0–2 (and the same for range 0–4), in Fig. 7. When we fit data in the range 2–4, the slope increased ~15% to 0.0015, possibly suggesting a small increase in grain size in the denser areas, similarly to Martin et al. (2012), and Roy et al. (2013). With β = 2.0, and using range 0–4, the slope increased ~38% to 0.0018.

For comparison with other studies, we convert our slope of to dust absorption cross-section per H nucleon, also called opacity. We use the conversion factor N(HI + H₂)/E(B − V) = 5.8 × 10²¹ cm/mag of Bohlin et al. (1978). In principle, this relation is valid only for diffuse areas, and can be used only as an approximation for the densest parts of our filament. We derive the extinction relations, E(B − V)/E(J − K) and A_J/E(J − K), using the extinction curves of Cardelli et al. (1989), either with R_V = 3.1 or R_V = 4.0. This leads to σ_e(250 μm) values 1.7 × 10^-25 cm²/H (with R_V = 3.1) or 2.4 × 10^-25 cm²/H (with R_V = 4.0). In terms of mass absorption (or emission) coefficient per gas mass, κ_ν, the values are 0.07 cm/g or 0.10 cm/g, respectively. Changing the assumed R_V value in E(B − V)/E(J − K) from 3.1 to 4.0 causes ~40% increase in the obtained values. It is not obvious which assumptions can be made in areas having both diffuse and denser parts, but it is important to notice that the assumptions made will have notable differences in the results.

For high-latitude, diffuse ISM, the standard value for dust opacity is σ_e(250 μm) ~ 1.0 × 10^-25 cm²/H (see, e.g., Boulanger et al. 1996; Planck Collaboration 2011). In denser areas, 2−4 times larger values have been derived (see, e.g., Juvela et al. 2011; Planck Collaboration 2011; Martin et al. 2012; Roy et al. 2013). Our estimates for σ_e(250 μm) are ~1.7–2.4 times as large as the previous results for diffuse areas, at the lower limit of the values for denser areas.

As the width (or FWHM) of our filament is ~0.1 pc, which for Taurus corresponds to  ~150, we have enough resolution to look at the details of the filament. The map shown in Fig. 8 gives no indication of systematic growth of the towards the densest filament. However, as discussed above, all the assumptions made in the process may not be valid in the densest part of the filament. As an additional confusing factor, some areas of high ratio can be caused by small differences in the τ₂₅₀ and maps, due to, e.g., the lack of background stars seen behind the densest filament, instead of real changes in the dust properties.

In Sect. 5, we constructed a realistic three-dimensional model of the density distribution of the densest part of the filament. We calculated predictions for the surface brightness using radiative transfer calculations of the ISRF light that is scattered from the cloud in the wavelength range 1.2–3.5 μm. We used DIRBE observations as a model for the sky brightness. As the absolute value of the background affects the contrast between the cloud and the background, we calibrated our simulation data to the same absolute level with DIRBE. We also considered the role of a potential diffuse envelope that would affect the radiation impinging on the dense part of the filament, i.e., the part included in our model. We take into account that a scattered photon has to pass through the diffuse cloud twice. We compared three cases of using different background correction: (1) no background; (2) with background; and (3) with background and attenuation of the scattered light in the envelope surrounding the filament.

In the simulations, we have assumed that the filament is in the plane of the sky and is not prolate (long along the line-of-sight). If the filament is prolate, this will increase the short wavelength surface brightness relative to the longer wavelengths, and the saturation of the J band will be smaller. If the filament is not in the plane of the sky (i.e., our line-of-sight is not perpendicular to the filament axis), the effect is similar. The dust opacities used in the simulations are appropriate in high density environments (Hildebrand 1983; Beckwith et al. 1990) but could, in diffuse regions, overestimate κ (underestimate column density) by a factor of two (see Boulanger et al.1996).

We find that the results change notably, especially for the J band, depending on what fraction of the background intensity is assumed to be really behind the filament and what fraction between the filament and the observer. In the J band, case 1 gives ~1.5 times as high values to the relation between I_ν and τ_J as case 2, and ~3 times as high values as case 3. In longer wavelengths the differences get smaller, but at 3.5 μm, case 1 still gives ~2 times as high values as case 2. At 3.5 μm, the difference between cases 2 and 3 is only ~10%.

In the H and K bands, the simulations give rather similar results to our observations in TMC-1N. However, in the J band, the simulations give over two times as high values as the observations, unless case 3 with the background correction and the correction for the scattered light is used. In the H band, case 3 also gives values that are lower than the observed values. The three test cases can be seen as rough error limits for the modelled cloud. Compared to the simulations, our observations do not suggest that there is any notable emission in the K band in addition to the scattered light.

7. Conclusions

We have used WFCAM NIR surface brightness observations to study scattered light in the TMC-1N filament in Taurus Molecular Cloud. We have presented a large NIR surface brightness map (1° × 1° corresponding to ~(2.44 pc)) of this filament. We have converted the data into an optical depth map and compared the results with NIR extinction and Herschel observations of sub-mm dust emission. We have also modelled the filament by carrying out 3D radiative transfer calculations of light scattering.

• We see the filament in scattered light in all three NIR bands, J, H, and K.

• In all three NIR bands, our WFCAM observations in TMC-1N show lower intensity than previous results in Corona Australis, indicating a lower radiation field in this area. This reaffirms the previous findings, that the radiation field in Corona Australis is at least three times that of the normal ISRF.

• We derive a value 0.0013 for the ratio . This leads to values σ_e(250 μm) ~ 1.7−2.4 × 10^-25 cm²/H, depending on the assumptions of the extinction curve (R_V = 3.1 or 4.0) which can change the results by over 40%. These σ_e(250 μm) values are twice the values reported for diffuse medium, at the lower limit of the values for denser areas.

• Changing β from 1.8 to 2.0 increases the ratio by ~30%.

• We see no indication of systematic growth of the ratio towards the densest filament. However, all the assumptions made in the process may not be valid in the densest part of the filament. Also, some areas of high ratio can be caused by imperfections in the τ₂₅₀ and maps, due to, e.g., the lack of background stars seen behind the densest filament, instead of real changes in the dust properties.

• 3D radiative transfer simulations predict surface brightness that is in intensity close to the observed values, especially in the H and K bands. In the J band, the model predictions can be over two times larger than observations, if no background correction is made. However, using background correction can change the results notably.

• We see no clear evidence for emission in the K band, in addition to the scattered light, based on the observations and simulations.

• NIR surface brightness can be a valuable tool in making high resolution maps, also at large scales.

• NIR surface brightness observations can be complicated, however, as the data can show comparatively low-level artefacts, that are still comparable to the faint surface brightness. This suggests caution when planning and interpreting the observations.

• It is possible to remove most of the effects of instrumental gradients, provided that they only affect large scales.

---

Acknowledgments

We thank the referee for useful comments. We thank CASA for carrying out the standard data reduction of the WFCAM data, and Mike Irwin for useful comments. J.M. and M.J. acknowledge the support of the Academy of Finland Grants No. 250741 and 127015. M.G.R. gratefully acknowledges support from the National Radio Astronomy Observatory (NRAO), the Joint ALMA Observatory and the Joint Astronomy Centre, Hawaii (UKIRT). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The United Kingdom Infrared Telescope is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the U.K. This paper makes use of WFCAM observations processed by the Cambridge Astronomy Survey Unit (CASU) at the Institute of Astronomy, University of Cambridge. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research made use of Montage, funded by the National Aeronautics and Space Administration’s Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. Montage is maintained by the NASA/IPAC Infrared Science Archive. This work is based in part 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. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

References

All Tables

All Figures