{'fraction_non_alphanumeric': 0.03524947632292847, 'fraction_numerical': 0.023848768929638404, 'mean_word_length': 4.236737996763172, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mapping of interstellar clouds with infrared light scattered from dust: TMC-1N

Mapping of near-infrared (NIR) scattered light is a recent method for the study of interstellar clouds, complementing other, more commonly used methods, like dust emission and extinction. Our goal is to study the usability of this method on larger scale, and compare the properties of a filament using NIR scattering and other methods. We also study the radiation field and differences in grain emissivity between diffuse and dense areas. We have used scattered J, H, and K band surface brightness WFCAM-observations to map filament TMC-1N in Taurus, covering an area of 1dx1d corresponding to ~(2.44 pc)^2. We have converted the data into optical depth and compared the results with NIR extinction and Herschel observations of submm dust emission. We see the filament in scattered light in all three NIR bands. We note that our WFCAM observations in TMC-1N show notably lower intensity than previous results in Corona Australis using the same method. We show that 3D radiative transfer simulations predict similar scattered surface brightness levels as seen in the observations. However, changing the assumptions about the background can change the results of simulations notably. We derive emissivity by using optical depth in the J band as an independent tracer of column density. We obtain opacity sigma(250um) values 1.7-2.4x10^-25 cm^2/H, depending on assumptions of the extinction curve, which can change the results by over 40%. These values are twice as high as obtained for diffuse areas, at the lower limit of earlier results for denser areas. We show that NIR scattering can be a valuable tool in making high resolution maps. We conclude, however, that NIR scattering observations can be complicated, as the data can show relatively low-level artefacts. This suggests caution when planning and interpreting the observations.


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.
Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. 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 nearinfrared 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 magnitudes.
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 near-infrared dust properties, the main sources of error are different. Comparison of the results can be used to study the values and variation of nearinfrared 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 ). 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. (2006Juvela et al. ( , 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.

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) +26 • 39 32 . At a distance of 140 pc, this corresponds to an area of ∼(2.44 pc) 2 . 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 magnitudes, 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 1 , including, e.g., dark-correction, flatfielding (including internal gain correction), decurtaining, sky correction, and cross-talk. Standard decurtaining methods 2 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 Daophotpackage, 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 3 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 m , the effect is small (less than 5% even in the densest parts), and we have not made this correction. 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 where P is the complete third order polynomial for a 2dimensional 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 ). 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 −biN ) 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 −biN ) − 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.

Herschel
The Taurus molecular cloud has been mapped with Herschel (Pilbratt et al. 2010) as part of the Gould Belt Survey 4 (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 τ 250 instead. We use the notation τ 250 to mean optical depth τ ν at wavelength 250 µm (ν =1200GHz). We derived τ 250 maps from Herschel intensity and colour temperature maps, using a value of 1.8 for the spectral index β. Recent studies with Planck (Planck Collaboration et al. 2011) show that this is a good estimate for Taurus. For 4 http://gouldbelt-herschel.cea.fr comparison, we also made the analysis with the commonly used value of β = 2.0.

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, artefactmitigated 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 5 . The final mosaic image was then made with the MOPEX tool (Makovoz et al. 2006).

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. (2006Juvela et al. ( , 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 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 ). 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 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 nonlinear 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 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 .
To derive an optical depth τ SB J map from the surface brightness maps, we obtain a and b parameter values for each band from I ν /τ N icer J correlations by fitting Eq. 2 to the data, using optical depth τ N icer J 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 τ SB J map by minimising the squared sum of residuals where I i are the pixel values of the intensity maps, and 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.

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 21 cm −2 . 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.  The removed gradients were fitted to the frames excluding 10% wide borders of each frame.

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 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, τ N icer J , are shown for bands J, H, and K in Fig. 5. We fit Eq. 2 to the data, using optical depth τ N icer J 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, τ 250 , 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 Table 1. Standard deviations (σ) for the a and b parameters and their product a × b.  Table 2. Relative change in the a and b parameters and their product a × b when the upper limit of the fitted range is changed from τ J = 2.0 to τ J = 6.0. 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 τ 250 and τ N icer J 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 τ N icer J =0-4 is ∼0.0013. We also tested the possible change of slope in the low and middle τ N icer J 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 τ 250 /τ J = 0.0013 to opacity or dust emission crosssection per H nucleon 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 2 /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 (1200GHz).
We derive the τ 250 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 2 )/E(B − V ) = 5.8 × 10 21 cm −2 /mag of Bohlin et al. (1978). We derive the relation E(B − V )/E(J − K) from Cardelli et al. (1989)  The map of the ratio τ 250 /τ N icer J is shown in Fig. 8. The map shows no clear evidence for systematic increase of the ratio τ 250 /τ J with increasing density inside the densest filament. Some high value areas can be attributed to imperfections in the two maps, τ 250 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.

Optical depth derived from scattered light
As described in Sect. 3, in order to derive an optical depth τ SB (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 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.  Fig. 9. Maps of optical depth τ : Herschel τ 250 (resolution 40 ), NICER τ N icer J (resolution 60 ), and optical depth in J band, τ SB J , based on J, H, and K surface brightness maps (resolution ∼2.2 ). In the NICER map (middle frame) the densest areas where τ N icer J > 4 are marked with a blue contour. The white rectangle marks the reference area used for background subtraction in the analysis. The profile cross-sections shown in Fig. 12 are taken as a median from the area between the two red lines.

Spitzer data
We compare the IRAC 3.6 µm surface brightness with the Herschel optical depth τ 250 , 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 2 . 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 21 cm −2 defines an area of ∼77 arcmin 2 .
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 τ 250 are shown in Fig. 14. In the figure, we also show the median surface brightness that is calculated for τ 250 bins with a width of 0.002. We fit a robust least squares line to the individual surface brightness values in pixels where τ 250 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.

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 (Interstellar Radiation Field) model (Mathis et al. 1983)

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 threedimensional 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 6 (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) (Mathis et al. 1983) values are marked with blue crosses and a dashed line.
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 2 /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 3 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 with parameters ρ C = 4×10 4 cm −3 , 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 where I is the total observed signal, I sca the scattered signal coming from the cloud, I bg e −τ 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 −1 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 −1 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 −1 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 τ J d = 0.5, and calculate values 0.34, 0.2, and 0.0998 for τ H d , τ K d , and τ 3.5µm d , 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 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 7 . 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 τ N icer J 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  Fig. 19. Simulated MIR surface brightness I ν in 3.5 µm band as a function of τ J in the three test cases described in the text. The data are fitted with Eq. 2, and marked with the same notation as in Fig. 18. The black dashed line shows the upper limit for the data used in the fit. The obtained parameter values are marked in the figure. 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%.

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) 2 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 τ N icer J 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 τ N icer J . 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 τ 250 /τ N icer J 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;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 τ 250 /τ N icer J = 0.0013 to dust absorption cross-section per H nucleon, also called opacity. We use the conversion factor N (HI + H 2 )/E(B − V ) = 5.8 × 10 21 cm −2 /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 2 /H (with R V =3.1) or 2.4 × 10 −25 cm 2 /H (with R V = 4.0). In terms of mass absorption (or emission) coefficient per gas mass, κ ν , the values are 0.07 cm 2 /g or 0.10 cm 2 /g, respectively. Changing the assumed R V value in E(B−V )/E(J −K) from 3.1 to 4.0 causes ∼40% increase  Comparison of observed and simulated NIR surface brightness I ν in J, H, and K bands as the function of optical depth in J band, τ J . WFCAM data for TMC-1N are shown with a 2D histogram, the colour scale corresponding to the density of points, and a black line shows the fitted Eq. 2. For comparison, simulated data are shown in the same figures using the three test cases described in the text. We plot the simulated data points only for case 1 (black dots). The fitted function is marked with a green line (case 1: without background), blue line (case 2: with background), and red line (case 3: with background and attenuation of the scattered light in the envelope). In J band, black dashed vertical line shows the upper limit of fitting of both WFCAM and simulation data. In H and K bands, the black dashed vertical line shows the upper limit of fitting of WFCAM data, and green dashed vertical line the upper limit of fitting of simulated data (all three cases). To get a better fit, the limit for the observations is slightly higher. The obtained parameter values are marked in the figures with the same colours as the fitted lines.
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.
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 τ 250 /τ N icer J map shown in Fig. 8 gives no indication of systematic growth of the τ 250 /τ N icer ament. 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 τ 250 /τ N icer J ratio can be caused by small differences in the τ 250 and τ N icer J 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.

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) 2 ) 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 τ 250 /τ N icer J . This leads to values σ e (250µm) ∼ 1.7 − 2.4 × 10 −25 cm 2 /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 τ 250 /τ N icer J by ∼30%. -We see no indication of systematic growth of the τ 250 /τ N icer J 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 τ 250 /τ N icer J ratio can be caused by imperfections in the τ 250 and τ N icer J 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.
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.