© ESO 2020

1 Introduction

Dust is central to the physics of the interstellar medium (ISM) and important for the heating of the gas, the formation of H₂ molecules, and the shielding of more complex molecules from the interstellar radiation field. Dust emission is used as a tracer of the ISRF, star formation (SF) activity, and ISM mass, which are all affected by the dust properties and its abundance. It is thus essential to know the properties of interstellar dust that affect their light scattering and thermal dust emission.

Coreshine and cloudshine refer to excess signal detected in the infrared. Cloudshine is caused by the scattering of the interstellar radiation field (ISRF) from the clouds at near-infrared (NIR) wavelengths (Foster & Goodman 2006; Ysard et al. 2016), while coreshine is caused by scattered photons from deeper within the dense cores, which are visible in the mid-infrared (MIR; Steinacker et al. 2010; Pagani et al. 2010). Coreshine and cloudshine provide a way to study the growth of grains in the dense ISM (Ysard et al. 2018). Both are affected by changes in dust scattering efficiency, which may be related to the surface irregularity of grains, changing grain size or fluffiness, coagulation, and ice coating (Ossenkopf 1993; Stepnik et al. 2003; Ridderstad & Juvela 2010; Ormel et al. 2011; Ysard et al. 2013, 2016; Köhler et al. 2015; Min et al. 2016).

Dust is heated by stellar radiation from the ultraviolet (UV) to the visible, and the absorbed energy is radiated away in a range from MIR to far-infrared (FIR) and millimetre wavelengths. The observed variation of the MIR-to-FIR ratios is believed to be due to dust grain evolution, grain growth, and ice mantle formation with increasing density (Ormel et al. 2009; Boogert et al. 2015; Köhler et al. 2015). These changes are reflected in the dust spectral energy distribution (SED), which thus provides important clues on dust evolution processes such as grain coagulation and fragmentation (Compiègne et al. 2011).

Stars form from collapsing clouds of dense interstellar gas and dust, dust emission being an important tracer of the process. Filamentary structures are common in molecular clouds (MCs; e.g. Men’shchikov et al. 2010; André et al. 2010, 2014, 2019; Arzoumanian et al. 2011; Hennemann et al. 2012; Juvela et al. 2012a; Malinen et al. 2012; Palmeirim et al. 2013; Wang et al. 2015) and they fragment to subparsec-scale cores, which may subsequently lead to the formation of young stellar objects (YSOs; Kirk et al. 2013; Offner et al. 2014; Könyves et al. 2015). SF is studied with both dust and molecular line observations and the latter are essential for investigations of cloud stability, kinematics, and chemistry (Motte et al. 1998; Bergin & Tafalla 2007; Enoch et al. 2007; Pattle et al. 2017). However, dust is a central tool also in the study of SF processes. The far-IR dust emission, often approximated as modified blackbody (MBB) emission, traces not only the column density but through the dust temperature also the ISRF changes that are associated with the general SF activity (Sadavoy et al. 2013; Planck Collaboration XI 2014). At shorter wavelengths, under 100 μm, emission from hot dust is important for the detection and characterisation of young stellar objects (YSOs; Lada 1987; Benedettini et al. 2018).

High-latitude clouds (|b| > 30°) are a fairly rare class of interstellar clouds (Dutra & Bica 2002; McGehee 2008). They are typically nearby objects with low column densities and no star formation. There are only a handful of high-latitude clouds with more molecular material and some low-mass SF activity (Lynds 1962; Malinen et al. 2014). These are excellent targets for studying SF triggered by supernovae and stellar winds, as SF due to gravitational collapse is less likely (Elmegreen 1998; McGehee 2008). Because of the low levels of line-of-sight (LOS) confusion, they are good targets also for studies into the dust properties in interstellar clouds.

LDN 1642, also referred to as MBM 20 and G210.90–36.55 (Lynds 1962; Magnani et al. 1985; Juvela et al. 2012b), is one of the star-forming high-latitude clouds (Table 1). It is gravitationally bound and hosts three YSO systems (McGehee 2008; Malinen et al. 2014). LDN 1642 is part of a larger (> 4°) cometary HI cloud, with an over5° long tail towards the Galactic plane (Gir et al. 1994; Alcalá et al. 2008). The cloud is projected on the Orion-Eridanus bubble (Brown et al. 1995), with which it may be interacting (Lehtinen et al. 2004). It is also located ~10° from the reflection nebula IC 2118 (the Witch Head nebula) (Kun et al. 2001; Alcalá et al. 2008). The cloud structure and the large-scale magnetic field of LDN 1642 are linked, the magnetic field possibly affecting the cloud evolution (Malinen et al. 2016). There is a clear change from magnetic-field-aligned to perpendicular structures around a column density of N_H = 1.6 × 10²¹ cm⁻² (Malinen et al. 2016). The light scattering in LDN 1642 at optical wavelengths has already been studied in Mattila et al. (2007, 2018).

LDN 1642 contains several denser regions, named by Lehtinen et al. (2004) as A1, A2, B, and C. Malinen et al. (2014) divide region B into two subregions, B1 and B2, due to an intensity maximum separate from the main clump. Three pre-main-sequence objects are associated with LDN 1642. Two of these, IRAS 04327-1419 = L1642-1 (V* EW Eri, HBC 413) and IRAS 04325-1419 = L1642-2 (HBC 410) were discovered by Sandell et al. (1987). We refer to these sources as B1 and B2, respectively, and their locations are indicated in Fig. 1. Both are faint binary stars. The primary of L 1642-1 is classified as a Type II YSO T-Tauri star of spectral class K7IV (Sandell et al. 1987; Malinen et al. 2014) and the primary of L 1642-2 as a flat-spectrum YSO M0 class Hα-emission star (Liljestrom et al. 1989; Malinen et al. 2014). A weak, bipolar outflow has been found around B2, and a Herbig-Haro object (HH123) originates from it (Liljestrom et al. 1989; Reipurth & Heathcote 1990). The 2MASS point source 2MASS J04351455-1 414 468 was originally classified as a potential foreground dwarf star (Cruz et al. 2003) but is now considered to be a Type III YSO associated with LDN 1642 (Malinen et al. 2014). We refer to this object as B3 (Fig. 1).

In this paper, we study the cloud LDN 1642 by combining Herschel satellite data with new NIR observations. With the help of radiative transfer (RT) modelling, we test the ability of dust models to consistently predict in LDN 1642 all three aspects of dust observations: sub-millimetre emission, NIR scattering, and NIR extinction. Its location at a small distance and high above the Galactic plane makes LDN 1642 a good target for this study.

The contents of the paper are the following. The observations at sub-millimetre, NIR, and optical wavelengths are presented in Sect. 2 and the main observational results in Sect. 3. Section 4 describes the RT modelling, where results are shown for dust emission in Sect. 4.1 and for NIR scattering in Sects. 4.2 and 4.3. We discuss the results in Sect. 5 before presenting the final conclusions in Sect. 6.

2 Observations

2.1 Dust emission

We used the pipeline-reduced Herschel observations from the Herschel science archive¹. Of the observations made with the PACS instrument (Poglitsch et al. 2010), we use the 160 μm data (level 2.5 data products, observation ID numbers 1342225212 and 1342225213), the maps produced with the Scanamorphos algorithm (Roussel 2013). The PACS 100 μm data show little extended emission and even the embedded sources are associated with little extended emission. The observations with the SPIRE instrument (Griffin et al. 2010) cover the wavelengths 250, 350, and 500 μm (observation ID 1342216940). The data correspond to extended-source calibration and, with a comparison to Planck data, have already been zero-point corrected by the pipeline (see Bernard et al. 2010). However, we analysed the data using background subtraction, which makes the results independent of the zero-point accuracy. The background values were estimated as the average intensity within 3.7 arcmin of the position RA = 4:34:49.35, Dec = –14:26:21.70. For the background determination, the maps were also first convolved to a common 41′′ resolution. Here, and later in the analysis, we use for SPIRE the convolution kernels provided by Aniano et al. (2011).

2.2 Near-infrared observations

The central part of LDN 1642 has been imaged in the J, H, and K_S bands using the HAWK-I instrument. HAWK-I is a cryogenic wide-field NIR camera installed at the ESO VLT telescope (Kissler-Patig et al. 2008). The field of view is 7.5′ × 7.5′, with a cross-shaped gap of 15′′ between the four HAWAII 2RG 2048 × 2048 pixel detectors. The pixel scale is 0.106′′ pixel⁻¹. Further details of the instrument can be found in Kissler-Patig et al. (2008). The observations were performed as ON–OFF measurements to recover the faint surface brightness. The observations consisted of three pointings arranged around the source B1, which was too bright for direct observations. The NIR photometry was done using the APPHOT task in Image Reduction and Analysis Facilities (IRAF) software, and the final calibration was provided by the comparison with the magnitudes in the 2MASS catalogue (Skrutskie et al. 2006).

To study the surface brightness, we created another set of images where the stars were eliminated. To remove the stars, we first used the DAOPHOT task ALLSTAR, which in many cases leaves significant residuals at the location of bright stars. The effect of some bright stars extends beyond the area masked by DAOPHOT. In such cases, the masks were extended manually, removing areas where the surface brightness enhancement was visibly above the general background. Finally, faint stars that were not identified by DAOPHOT were removed with median-filtering. The size of the median filter was 5.0′′ (or 19 pixels) for J, H, K_S, and WISE 3.4 μm bands.

For the J- and K_S-bands we calculated estimates of the absolute sky brightness behind the LDN 1642 cloud by subtracting from DIRBE measurements the combined flux of 2MASS stars, weighted by the DIRBE beam. This procedure gave 71 and 31 kJy sr⁻¹ for the J and K_S bands, respectively. A linear interpolation gives 53 kJy sr⁻¹ for the H band. The values have considerable uncertainty since the estimates vary by ~30% when derived from neighbouring independent DIRBE pixels. When combining the DIRBE and 2MASS data, we do not explicitly include colour corrections that are small compared to this uncertainty (Levenson et al. 2007).

Fig. 1 Column density N(H2) (frame a), dust temperature Tdust (frame b), and J-band extinction AJ (frame c) of LDN 1642. The resolution of the maps is indicated by the circles in the lower left corner. The values of N(H2) and Tdust are based on MBB SED fits to Herschel SPIRE data, after background subtraction. The red circle in frame a shows the reference region used for the background estimation. Areas with low column density have been masked in the plots (white pixels). The labelled stars indicate the locations of embedded sources B1–B3.

Fig. 2 Large-scale environment of LDN 1642 shown by dust emission and optical scattered light. Left panel: Planck 857 GHz surface brightness, the white arrow indicating the direction towards the Galactic Plane. Centre panel: optical scattered light from DSS red channel. The white box shows the location of the LDN 1642 cloud and covers an area of 0.7° × 0.7°. Right panel: a close-up of the area marked with the white box in the centre panel, showing the extended R-band surface brightness. The black stars indicate the locations of the embedded sources B1, B2, and B3, and the red contours show the N(H2) column density. The lowest contour is at 15 % of the peak column density of 5.14 × 1021 cm−2, with the contour levels increasing in 15 % intervals.

2.3 Optical observations

The optical data were obtained from the Infrared Science Archive (IRSA)² and are part of the Space Telescope Science Institute (STScI) Digitized Sky Survey (DSS). LDN 1642 was imaged with the UK Schmidt telescope to a photographic plate (emulsion type IIIaF) using an OG590 filter, which roughly corresponds to the R band, covering the wavelength range 6300–6900 Å. The plate covers a ~6° × 6° area on the sky and the digitised images have a pixel size of 1′′. The DSS data are shown in Fig. 2. The conversion to units of MJy sr⁻¹ is described in Appendix A.

There are two bright stars, HD29613 (red giant) and HD29503 (1 Eri, a triple star), at a small angular distance from LDN 1642. With the parallax measurements in the Gaia (Gaia Collaboration 2016) Data Release 2 catalogue (DR2, Gaia Collaboration 2018) (see also Bailer-Jones 2015; Luri et al. 2018) we estimate ~ 65 pc and ~ 36 pc as the distances of the two sources, respectively. The cloud is at a more than 70 pc larger distance, and the stars are thus not likely to contribute significantly to its surface brightness.

3 Results

3.1 Column densities and extinction

Column densities and dust temperatures were estimated using the Herschel 160–500 μm surface brightness measurements. The maps were resampled onto common pixels and modified blackbody fits were performed pixel by pixel. The dust opacity spectral index was fixed to β = 1.8 (Planck Collaboration XXV 2011; Juvela et al. 2015a) and the conversion from optical depth to hydrogen column density assumes a dust opacity of κ_ν = 0.1 g cm (Beckwith et al. 1990). This value of κ_ν corresponds to dust properties in very dense regions and has been used in earlier analysis of LDN 1642 (e.g. Juvela et al. 2012b).

One set of calculations was done with 160–500 μm maps convolved to a common 41′′ resolution, providing dust colour temperature and column density maps at the same resolution (Fig. 1a). Another column density map was calculated at a higher resolution following the procedure described in Palmeirim et al. (2013), combining estimates computed with 160–250, 160–350, and 160–500 μm data. This procedure provides a N(H₂) map at the resolution of the 250 μm observations that was further smoothed to 25′′ resolution.

Extinction maps were calculated using the combination of 2MASS (Skrutskie et al. 2006) data and our HAWK-I photometry. For the latter, we included all stars with measured magnitudes at least in the H and K_S bands. We used a variation of the NICER method (Lombardi & Alves 2001) that takes into account the estimated ratio between the average extinction within a resolution element and extinction towards individual stars (Juvela & Montillaud 2016). These ratios were read from the Herschel column density map. The resolution of that map is higher than the resolution of the final extinction maps (41′′ vs. 2′), which helps to reduce the noise caused by column density variations on scales below 2′. Figure 1c shows the resulting extinction map of τ(J). The extinction peaks at A_J = 2.6 mag (A_V = 9.3 mag for R_V = 5.1)).

3.2 Comparison of sub-millimetre and near-infrared data

To quantify the ratio of sub-millimetre and NIR opacities, we calculated the ratio τ(250 μm)∕τ(J). The first estimates are based on the NICER extinction map at 2′ resolution and the Herschel column density map convolved to the same resolution. To establish a common zero point, the average value in region defined by τ(250 μm) < 10⁻⁴ was subtracted from both maps. With data remaining above τ(250 μm) = 10⁻⁴ and sampled at 1 arcmin steps, we obtained a ratio ⟨τ(250 μm)⟩∕⟨τ(J)⟩ = 1.07 × 10⁻³. Figure 3a also shows a linear total least squares fit to part of the data (blue points). The slope of the fit gives τ(250 μm)∕τ(J) = (1.22 ± 0.04) × 10⁻³, with the formal error estimates from the least squares fit.

The above values are based on extinction maps at low resolution. For comparison, we also correlated the extinction estimates of individual stars with the column densities read from a map with 40′′ resolution. This fit is shown in Fig. 3b. The result was τ(250 μm)∕τ(J) = (0.95 ± 0.04) × 10⁻³. Thus, with more conservative error estimates, the optical depth ratio in LDN 1642 is τ(250 μm)∕τ(J) = (1.0 ± 0.2) × 10⁻³. Even this uncertainty may not fully cover all systematic errors such as the possible bias in τ(250 μm) caused by LOS temperature variations.

The derivation of τ(J), as shown in Fig. 3, assumed the standard extinction curve of Cardelli et al. (1989). We also correlated the N(H₂) data at 25′′ resolution directly with the NIR J–H and H–K colours of individual stars (Fig. 4). The relationships remain linear up to the highest values and for both colour excesses. The N(H₂) uncertainties were assumed to be 20% and the uncertainties of the NIR colours are the squared sums of the photometric errors and a constant that represents the dispersion in the intrinsic colours and is set so that the final estimates of uncertainty are consistent with the observed scatter. The numerical values quoted in the figure depend on the chosen value of κ(250 μm), but the ratio of the slopes gives independently E(H − K)∕E(J − H) = 0.73 ± 0.35. This value is slightly higher than in the Cardelli et al. (1989) extinction curve (0.66), but the difference is not significant considering the uncertainties. The result also would change by less than 1% if N(H₂) was taken from the other column density map, which was derived from the 250–500 μm Herschel data at 41′′ resolution. The extrapolation of the linear relationships to zero column density shows the average NIR colours in the region used for the N(H₂) background subtraction, ⟨J − H⟩ = 0.57 mag and ⟨H − K⟩ = 0.17 mag.

Fig. 3 Correlations between the optical depths at 250 μm and in the J band. The left frame uses a NICER τ(J) map with 2′ resolution.The pixel values are plotted as black points. The blue points correspond to the same data but excluding pixels close to the background level. The dashed line is the least-squares line to the blue points and, for comparison, the solid line has a slope corresponding to the ratio of the averages of the background-subtracted values. In the right hand frame, we compare τ(J) estimates of individual stars to the Herschel column density estimates at 40′′ resolution. The parameters of the least-squares lines are given in the frame.

Fig. 4 Estimated column density N(H2) as function of NIR colours J–H (in blue) and H–KS (in red). The parameters of the linear fits are given in the plot.

Fig. 5 Observed surface brightness in J, H, KS, and WISE 3.4 μm bands (frames a–d, respectively). Herschel column density and dust temperature maps are shown in frames e and f. Bright stars identified by DAOPHOT were removed and replaced with interpolated surface brightness, and faint stars have been eliminated with median filtering. In the WISE 3.4 μm band image (frame d), the white regions correspond to areas around bright stars that were excluded from subsequent analysis. The YSOs B1–B3, are identified in frame a. The yellow circle in frames a–d indicate the reference region used for setting a common zero level. The white box in frame e corresponds to the extent of the J, H, KS, and 3.4 μm maps (the area shown in frames a–d).

3.3 Scattered light

In this section, we compare the column density estimated from sub-millimetre dust emission with the surface brightness maps for NIR dust scattering. The final surface-brightness maps of the central part of the LDN 1642 cloud in the J, H, K_S, and WISE 3.4 μm bands are shown in Fig. 5, which also includes for comparison the column density and dust temperaturemaps derived from Herschel observations. All NIR maps have been convolved to the same 25′′ resolution, to enable comparison between the NIR surface brightness and column density. To establish a common zero point, we subtracted the median values of a region marked with the yellow circles in Fig. 5. The surface brightness values for J, H, and K_S are similar within a factor of two, while for WISE 3.4 μm, the values are almost ten times lower. LDN 1642 contains a central dense region surrounded by diffuse material. The densest part of the cloud contains three YSOs, B-1, B-2, and B-3. These three YSOs are detectable in J, H, K_S, and 3.4 μm maps as compact sources. North of B1, an elongated structure is visible in surface brightness and Herschel column density maps. It is more prominent in J and H compared to the longer wavelengths, K_S and 3.4 μm.

Sub-millimetre dust emission shows good correlation with the morphology of NIR emission as shown in Fig. 5. The correlation between the observed scattered light in the J, H, K_S, and 3.4 μm bands with Herschel column density is shown in Fig. 6. To make a pixel-to-pixel comparison of observed J, H, K_S and 3.4 μm surface brightness with Herschel column density maps, we selected two regions (marked in Figs. 7e and 8f). Within the larger region, we masked the positions of the three YSOs, B-1, B-2 and B-3 to reduce their effect. The maps at 25′′ resolution are sampled at 7′′ steps. In Fig. 6, the blue points show the comparison for the larger area and the red points for the smaller area. We fitted robust linear least squares lines for both regions. The observed NIR surface brightness shows strong correlation with Herschel column density in all the bands. The slope decreases for both regions as the wavelength increases from J-band to 3.4 μm (Fig. 6). Table 2 lists the parameters estimated from the comparison between the Herschel column density map and the surface brightness maps. The correlation coefficients indicate strong correlation in the J, H and 3.4 μm bands and moderate correlation in the K_S band.

Fig. 6 Correlation between the observed scattered light in the J, H, KS, and 3.4 μm bands and Herschel column density. The data are at 25′′ resolution and sampled at 7′′ steps from the two areas marked in Fig. 8f, blue points for the larger and red points for the smaller area. The dashed lines and the dotted lines are robust linear least squares lines fitted to the red points and blue points, respectively.The fitted slopes and intercepts are quoted in the frames. The Spearman’s correlation coefficients are also given in the corresponding colours.

4 Radiative transfer modelling

We constructed radiative transfer (RT) models for the dust emission to derive 3D models of the density distribution in LDN1642. Next, the obtained density field and assumptions of the external and internal radiation sources were used to calculate predictions for the scattered light.

4.1 Modelling of Herschel emission

A 3D density model of LDN1642 was first optimised to match the Herschel 250–500 μm observations. The cloud volume was divided to 144³ cells, each with a linear size of ~0.0136 pc. This length scale corresponds to 20′′ at the distance of 140 pc. The model was further refined according to the local density by adding up to two levels in the octree hierarchy, with an approximately equal number of cells on each of the three refinement levels. The model thus reaches a resolution of 5.0′′ over most of the dense areas.

The LOS density profile corresponding to a map pixel (along the third dimension) was set equal to the narrowest column density profile that existed for any line (any position angle) crossing that pixel on the plane of thesky. This default LOS profile is in the following referred to as having a relative width of W = 1. It favours a cylindrical 3D geometry for features that appear elongated on the observed surface brightness maps. Because of the fundamental difference between the density and projected column density profiles, this setup corresponds only approximately to cylinder symmetry. Furthermore, because the actual LOS extent of the cloud and the inclination of the structures are unknown, we tested for comparison models with 70% larger LOS extent (W = 1.7).

The angular distribution of the incoming radiation was obtained from the DIRBE all-sky maps (Hauser et al. 1998), which include direct observations of the J and K_S bands. At shorter wavelengths, the angular distribution is assumed to be the same as in the J band. This assumption is not entirely accurate but includes the main effects of the anisotropic external illumination. The spectrum of the external radiation field was rescaled to match the Mathis et al. (1983) estimates of the radiation field intensity in the solar neighbourhood. In the NIR regime, the Mathis et al. (1983) estimates are some 40% below the sky-averaged DIRBE values.

For the first models, the dust properties were taken from Jones et al. (2013, in the following J13). In dense clouds, such as LDN 1642, the dust is expected to evolve towards larger grain sizes and larger sub-millimetre opacity. To quantify the potential effects of this dust evolution, as a second option we tested models that contained only ice-coated aggregate grains (AMMI). A third set of model clouds was also created using spatial variations in the relative abundance of J13 dust, core–mantle–mantle (CMM) grains, and AMMI grains. In the following, these models are referred to as THEMIS models. Ysard et al. (2016) have already used CMM and AMMI dusts to model enhanced NIR and MIR scattering (cloudshine and coreshine). In the THEMIS cloud models, the relative abundances of the three dust components were set according to the function (1)

The threshold density values (n₁, n₂) were set equal to (1, 10⁴), (10⁴, 4 × 10⁴), and (4 × 10⁴, 10¹⁰) cm⁻³ for J13, CMM, and AMMI dust, respectively. Thus, the low-density parts of the model cloud have only J13 dust. Its abundance drops to zero around 10⁴ cm⁻³, where CMMis briefly the most abundant component before the relative abundance of AMMI rises from zero to one at densities above 4 × 10⁴ cm⁻³.

The calculations were performed with the Monte Carlo RT programme SOC (Juvela 2019), assuming that sub-millimetre emission can be predicted with calculations where the grains remain at an equilibrium temperature. Because of the associated larger computational cost, full calculations with stochastically heated grains were performed, for reference, only in one case (J13 dust, W = 1).

The model predictions for surface brightness were saved as 288 × 288 pixel maps with a pixel size of 5′′. To match the model intensities with the Herschel observations, the models were optimised iteratively. The columndensities were updated using the ratio of the observed and the model-predicted 350 μm surface brightness values. This gave for each map pixel a scaling factor that was used to update the densities in all cells along the corresponding LOS. The external radiation field was adjusted with a single scalar factor equal to the average 250–500 μm intensity ratio in observations, divided by the same ratio in the model predictions. This approach converges the average colour temperature of the model towards the observed average 250–500 μm colour temperature. In addition to the external radiation field, we included point sources at the locations of the three known embedded sources (see Fig. 1). These were modelled as T = 6000 K black bodiesand their luminosities adjusted so that the model predictions for the average dust colour temperature in a small 5 × 5 pixel area around the sources matched the observations. The heating from the embedded sources reduces the predicted column densities but only in a very limited area. When the scattered light is later modelled in the neighbourhood of the embedded sources, the density field predicted by this emission modelling is not used (Sect. 4.3).

The comparison of the observed and model-predicted surface brightness maps is shown in Appendix B.1. The emission is fitted almost equally well using any of the dust models and with both cloud shapes, W = 1 and W = 1.7. Figure 9 shows the abundance variations in the THEMIS model, by plotting the column densities weighted by the relative amounts of the J13, CMM, and AMMI dust.

The emission models predict NIR optical depths τ(J) that are above NICER estimates, that is to say compared to the direct extinction measurements based on background stars. To quantify this difference, we convolved the τ(J) model maps to the 2′ resolution of the NICER map and subtracted from both the average value in the area where τ(J) of the J13 model was within the 1–5% percentile range. The ratios of the average model-predicted and observed τ(J) were then calculated for the pixels falling in the 25–90% percentile range of the J13 map. The lower limit was chosen so that the comparison avoids regions where the signal is close to zero, and the upper limit of 90% was chosen to downweight the contribution of the cloud centre, where the low stellar density renders the NICER estimates more uncertain. However, the obtained optical-depth maps were quite flat, without systematic variations correlated with the column density.

The results are listed in Table 3. The NICER estimates are always calculated using the extinction curve of the corresponding dust model, although these change the results only a little. The NIR optical depths of the optimised model clouds are 2.3–4.6 times higher than observed. The differences depend mainly on the dust model. The LOS cloud extent has effects only at a level of a few per cent. The full treatment of stochastic grain heating increased the τ(J) value by only 25% and the estimate of the radiation field intensity by some 40%. Qualitatively, the effect of stochastic heating is expected to be similar for the other dust models. The three-component THEMIS model is in this comparison close to the J13 model, because the J13 dust component is dominant at large scales. The relative abundance of the dust components is illustrated in Fig. 9.

The differences in the predicted volume density and dust temperature distributions is discussed in Sect. 5.3.1. We return to the optical depth discrepancy between the emission models and the NICER maps in Sect. 5.3.2.

Fig. 7 Predicted intensity of scattered light. The cloud model is J13 with W = 1. Frames a–f: bands R, J, H, KS, 3.4, and 4.6 μm, respectively. Only the scattering of the external radiation field is included and the effects of attenuated sky background are not included. Frame e indicates areas used for correlation analysis, the larger area drawn with thin line (excluding the sources B1-B3) and the smaller area with thick line (see also Fig. 8).

Fig. 8 Comparison of observed (first row) and modelled (other rows) 1.25–3.4 μm surface brightness. The second row corresponds to J13 dust, W = 1, and no sky background (Ibg = 0). The effects of sky background are included in the other cases (J13 dust with W = 1 and W = 1.7, and AMMI with W = 1, and the THEMIS model with three dust components). Frame f indicates the areas used for correlation analysis, the larger area shown with a thin cyan line (excluding the sources B1–B3) and the smaller area with a thick cyan line. The yellow circles indicate the reference region used for setting a common zero level. Black contours are drawn at 1.5 times the maximum of the colour scale. The colour scale is the same for all models but different for observations.

Fig. 9 Hydrogen column density N(H2) associated with each dust component in the THEMIS model with spatially varying dust abundances. The three frames correspond to the J13, CMM, and AMMI dust, respectively. The black circle shows the reference area used in scattered-light analysis (cf. Fig. 5).

4.2 Modelling of extended scattered light

We calculated predictions for the scattered light in the R, J, H, K_S, and 3.4 μm bands, using the density fields obtained from the emission modelling (Sect. 4.1). Figure 7 shows the results for one of the cases (J13, W = 1). The figure includes the 4.6 μm band although, for the lack of signal in the actual observations, this wavelength will not be analysed any further. The figure only shows the scattered light and thus does not include the effects of the LOS sky background, which reduces the observed ON–OFF signal by I_bg(e^−τ − 1).

The R-band results are compared to DSS data in Appendix A. We do not have an estimate for the absolute brightness of the sky background in the R band and therefore only show the comparison with the scattered-light component from the model. That emission was found to be only a fraction of the observed signal and a positive sky background I_bg would decrease the model predictions further. The R-band data are not analysed further in this paper.

Figure 8 compares the model predictions at J, H, K_S, and 3.4 μm with the corresponding HAWK-I and WISE observations. After median filtering, the observed maps were convolved to 25′′ resolution. The model predictions have a similar resolution because the density field was fitted to Herschel 350 μm data that have ~ 26′′ resolution (Sect. 4.1). This resolution applies to the scattering of the external ISRF. The direct radiation from the embedded sources and their scattered radiation were convolved with FWHM = 25′′ Gaussian, although the effective resolution of the latter is not well defined. Near the point sources, the surface brightness variations are dominated by the radial change of the radiation field, which is to some extent resolved at a resolution higher than that of the underlying density field. However, we examine the environment near the embedded sources separately in Sect. 4.3 and here concentrate on the larger scales. At distances larger than ~ 1′ from the sources B1-B3, the scattering is dominated by the external ISRF. To establish a common zero point for the observations and the models, we subtracted from each map the median value of the reference area indicated in Fig. 8.

Figure 8 includes model predictions for five cases. The first three are for J13 dust, the first one showing the scattered light without the effect of the background sky brightness I_bg on the ON–OFF measurement. The background term I_bg(e^−τ − 1) is included in all other cases, clouds with W = 1.0 and W = 1.7 and with J13 dust, and the W = 1.0 cloud with AMMI and THEMIS dust cases. The predicted surface brightness is significantly below the observed levels and is even negative for the J band. The negative values result from the significant sky brightness, combined with the significant NIR optical depth of the model clouds.

Because of the tension between the τ(J) values of the NICER measurements and the cloud models that fit the dust emission, we repeated the scattering calculations using modified model clouds where the average τ(J) was decreased to the level of the observed NICER values. This was done by dividing all densities by the factors / of Table 3. The results are shown in Fig. B.3 and the correlations for the model J13 are plotted in Fig. 10. The intensity of the J13 model is still far too low but the pure AMMI model rises almost to within a factor of two of the observed surface brightness values. The interpretation of these results is discussed in Sect. 5.3.2.

Fig. 10 Correlation between the observed and modelled scattered light in the J, H, KS, and 3.4 μm bands. The model results correspond to models where the densities have been scaled down to match the NICER extinction measurements (J13 dust, and W = 1). The data have a resolution of 25′′ resolution and sample at 7′′ steps the areas marked in Fig. 8f. The blue points correspond to the larger area and the red points to the smaller area. The dashed lines show robust least squares fits to the red points. The values of the slopes and correlation coefficients r are given in the frames.

4.3 Scattering around embedded sources B2 and B3

The large-scale model of Sect. 4.2 does not have the resolution to accurately describe scattering near the embedded sources. Therefore, we made separate, spherically symmetric models for the sources B2 and B3. The source B1 was not covered by H and K_S observations and is not considered here. For B2, the observations are partly saturated but only within the innermost couple of arcsec.

The spherical models have a spatial resolution of 0.5′′ (70 au at the distance of 140 pc) and extend to a distance of one arcmin. Section 4.2 showed that the surface brightness produced by the external radiation field is relatively uniform at this scale. Therefore, we computed predictions for the scattered light from the spherically symmetric models ignoring the external illumination. When the models are compared to observations, we subtract from both the average signal at 55–60′′ radial distances. The comparison will thus not be affected by the external radiation field if its contribution is significant only at distances larger than ~ 1′ or if it can be approximated as a flat background. The effects of the LOS sky background are also not considered. First, the intensity of the scattered light within the 1′ region is high compared to the sky background. Second, the effects of the I_bg are decreased by a factor of e^−τ, where τ refers to other LOS extinction, if that is uncorrelated with the structure inside the 1′ region.

The models were fitted to observations by modifying the radial density profile and the source luminosity. For the radial density profiles, we tested both truncated power laws and Plummer-like functions, with no significant difference in the fit quality. We show results for the Plummer functions with three free parameters, that is the centre density n₀, the characteristic radius R, and the asymptotic powerlaw exponent α, (2)

The point sources were modelled as 6000 K black bodies, with the total luminosity as a free parameter. Since the scattered light is linearly proportional to the intensity of the light sources, the results can be easily rescaled for any assumption of the spectrum of the central source.

Figure 11 shows the results for the B2 and B3 regions (model J13), including the attenuated direct radiation from the sources. Model maps are convolved with the point spread function (psf) estimated from the observations of unsaturated stars in the field. The model parameters are fitted using data at 5–40′′ radial distances.

The figures show a good match between the observations and the model, which in the 10–40′′ distance range mainly consists of scattered light. Towards the centre, the signal is dominated by the direct light from the point source and unresolved scattering. Also this part is well matched, the intensity profiles following the shape of the point spread function. The psf was estimated up to a radial distance of 15′′ but there it is already orders of magnitude below the peak value and cannot be measured reliably. In the outer part, beyond 40′′, the final drop is caused by the background subtraction.

The bestfit density profiles are practically flat with α≳ − 0.15 and the observed surface brightness is almost consistent with a radial decrease of the radiation field intensity in a homogeneous medium. The LOS optical depths from the source to the observer are for B2 below one and for B3 of the order of one. The contribution of thermal emission should be small for the assumed 6000 K sources (Sellgren et al. 1996), because of the low colour temperature and because of the large extinction that further removes short-wavelength photons from the radiation field. The fitted luminosities (see Fig. 11) suggest higher effective temperatures for the central sources, but these values should be sensitive to the assumed structure and opacity of the dust layers closest to the central source. These are not well constrained by the available data. Because the relative contribution of thermal emission should decrease with radial distance (as short-wavelength photons are absorbed), it should result in steeper density profiles in the models where only the scattered light is considered. Therefore, the fact that the derived density profiles were flat argues against any significant contribution from dust emission.

Appendix B shows corresponding plots for models with AMMI dust. Because of the higher dust albedo, the derived source luminosities are lower by almost a factor of two but the derived density profiles are still flat. For B2, the powerlaw exponent is smaller, α = −0.97, but the characteristic radius is large, R₀ = 0.016 pc, which corresponds to ~ 26′′. Thus, although the fit parameters are partly degenerate, all models show only little radial density variation. The degeneracy also applies to the parameters n₀ and R₀.

Fig. 11 NIR surface brightness profiles in the vicinity of the embedded sources B2 (frames a–c) and B3 (frames d–f). The black curvesshow the azimuthally averaged observed profile. The red dashed curves show the model predictions for J13 dust, consisting of the scattered light and the attenuated direct contribution of the central source. The average signal at 55–60′′ distance is subtracted from both curves. The cyan lines show the estimated direct contribution of the point source, based on the fitted source luminosity, the LOS extinction, and the psf shape.

5 Discussion

5.1 Sub-millimetre dust opacity

We derived with HAWK-I NIR observations and Herschel sub-millimetre data three estimates for the submm-to-NIR extinction ratio that were consistent with τ(250 μm)∕τ(J) = (1.0 ± 0.2) × 10⁻³. The comparison of data at 2′ resolution gave ⟨τ(250 μm)⟩∕⟨τ(J)⟩ = 1.07 × 10⁻³ while the slope of the least-squares fit gave τ(250 μm)∕τ(J) = (1.22 ± 0.04) × 10⁻³. The latter fit excluded low-column-density regions and the higher value could thus be more representative of the central regions of LDN 1642. As discussed in Juvela et al. (2015b), both optical depth estimates τ(250 μm) and τ(J) could be biased. However, because of the modest optical depths, the bias in τ(250 μm) should remain small. The maps of τ(250 μm)∕τ(J) were quite flat, also suggesting that the ratio could be measured reliably.

LDN 1642 was in the sample of cold clumps analysed in Juvela et al. (2015b) as source G210.90–36.55, where the combination of Herschel and 2MASS data resulted in an higher value of τ(250 μm)∕τ(J) = 1.6 × 10⁻³. This value was close to the median ratio for a sample of clouds extracted from the Planck survey of cold clumps (Planck Collaboration XXVIII 2016). The sub-millimetre opacity was obtained from MBB fits with β = 2.0 and with β = 1.8 (as used in the present paper), the value would decrease by some 20–30%, becoming marginally consistent with our new results.

The ratio τ(250 μm)∕τ(J) of LDN 1642 is significantly hi- gher than in the diffuse medium. Assuming the R_V = 3.1 extinction curve and the ratio N(H₂)∕A_V = 9.4 ×10²⁰ cm⁻² mag⁻¹ (Bohlin et al. 1978) our result corresponds to τ(250 μm)∕N_H = 1.5 ×10⁻²⁵ cm² H⁻¹. This value is three times higher than the Planck measurement τ(250 μm∕N_H) ~ 0.5 × 10⁻²⁵ cm² H⁻¹ obtained at high latitudes. We scaled the value to 250 μm using the opacity scaling ν^1.53 given in that paper (Planck Collaboration Int. XVII 2014). Fukui et al. (2014) examined atomic regions around selected high-latitude clouds and obtained a value τ(850 μm)∕N_H = 1.5 × 10⁻²⁶ cm² H⁻¹. For spectral indices β = 1.5–1.8, this results in similarly low 250 μm dust opacities, τ(250 μm∕N_H) = (0.4 − 0.6) × 10⁻²⁵ cm² H⁻¹.

The LDN 1642 value is similar to previous measurements of dense molecular clouds, indicating clear dust evolution relative to the diffuse medium (Stepnik et al. 2003; Martin et al. 2012; Roy et al. 2013). Below we mention some recent studies. When the original results were reported for different wavelengths, we assume β = 1.8 for the long wavelengths and the standard extinction curve for NIR (R_V = 3.1), to scale the results to τ(250 μm)∕τ(J).

Suutarinen et al. (2013) used Herschel data and dedicated NIR observations to derive directly a value τ(250 μm)∕τ(J) = 1.4 × 10⁻⁴ for the Corona Australis cloud. Lombardi et al. (2014) used both Planck and Herschel data to derive ratios of 850 μm and K_S-band opacity in Orion. The results correspond to τ(250 μm)∕τ(J) = 1.5 × 10⁻³ and 1.1 × 10⁻³ for the Orion clouds A and B, respectively. With similar analysis, Zari et al. (2016) found a value of 1.0 × 10⁻³ in the Perseus molecular cloud and Lada et al. (2017) found 1.1 × 10⁻³ in the California Nebula. Even larger relative increases of dust opacity have been reported at longer wavelengths (Mason et al. 2020). On the other hand, Forbrich et al. (2015) found for FeSt 1–457 (a core in the Pipe nebula) a dust opacity that was only slightly higher than in the diffuse medium, τ(250 μm)∕τ(J) = 0.65 × 10⁻³, in spite of the data covering extinctions up to A_K = 5 mag.

The high τ(250 μm)∕τ(J) values are more consistent with models of evolved dust, with increased grain sizes and the formation of aggregates, possibly covered by ice (Ossenkopf & Henning 1994; Ormel et al. 2011; Köhler et al. 2012; Ysard et al. 2016, 2019). Corresponding changes should exist also in the shorter-wavelength scattering properties of the grains.

5.2 Scattered light

The light scattered by dust grains has been observed in the NIR, referred as cloudshine, towards many molecular clouds (Juvela et al. 2008, 2012c; Lefèvre et al. 2014). However, recently discovered MIR dust scattering through Spitzer 3.6 μm and WISE 3.4 μm data, referred as coreshine, was the first direct indication of a significant population of ~1 μm grains in pre-stellar cores (Pagani et al. 2010, Steinacker et al. 2010, Juvela et al. 2012c). Juvela et al. (2008, 2009, 2012c) studied the Corona Australis cloud in detail and concluded that NIR scattered light can be used to estimate better resolution column density maps at low visual extinction (A_V≲10). They also found a linear relationship between Herschel column density estimates and NIR scattering.

Our data included HAWK-I NIR and WISE 3.4 μm measurements. We have assumed that the surface brightness is due to scattered light only. Figure 5 clearly shows extended emission from the densest part of the LDN 1642 cloud. Malinen et al. (2014) studied the LDN 1642 cloud at multiple wavelengths and suggested that the extended emission from the densest part of LDN 1642 in the WISE 3.4 μm image could be due to scattered MIR light (the coreshine phenomenon) associated with grain growth (Steinacker et al. 2010). We masked the bright sources in the WISE3.4 μm image and showed that the extended emission is similar to the NIR surface brightness and the column density maps. It can be considered an upper limit for dust scattering. No extended emission is found in the WISE 4.5 μm map. In thecomparison with the Herschel column density, we found a strong linear correlation in the J, H, and 3.4 μm bands. The correlation was weaker in the K_S band because of the smaller optical depth and problems with the data quality.

Lefèvre et al. (2014) studied dust grain properties inside molecular clouds using coreshine modelling and showed how the intensity of the coreshine depends on the incident radiation, the extinction of the background radiation, the grain properties, and the core properties. They found a higher NIR/MIR (K∕3.6μm) ratio for the Taurus–Perseus region, which could be explained by the presence of ice mantles. Alternatively, the grain size distribution having bigger silicates than carbonaceous grains could explain the higher NIR/MIR ratio (Lefèvre et al. 2014). For LDN 1642, the NIR/MIR (K/3.6 μm) ratio is ~8 (11) for the larger (smaller) region shown in Fig. 8f. These high values are consistent with that of the Taurus–Perseus region, in that interpretation suggesting bigger silicate dust grains (Lefèvre et al. 2014). The J/K ratio for LDN 1642 is ~1.6 (0.9) for the larger (smaller) region. These values are also similar to the typically observed range of 0.3–3 found by Lefèvre et al. (2014).

Fig. 12 Comparison of model column density N(H2) (first row), NIR optical depth τ(J) (second row), and a cross section of the physical dust temperature T(dust) (third row) in alternative models. The first column shows the reference model with J13 dust and W = 1 LOS cloud extent. The other frames show differences to this model, for the dust and W values indicated on top. The last column (J13, SHG) is the same as the reference model but with calculations including the full treatment of stochastic heating. Temperature maps are not shown for the THEMIS and stochastically-heated-grain cases, wherethere is no single temperature per cell.

5.3 Radiative transfer models

In this section, we discuss the interpretation of the results based on the RT modelling. We start with models fitted to the sub-millimetre emission and their consistency with NIR extinction. We then discuss the NIR surface brightness and the tension between the models matching observations of either the sub-millimetre dust emission or the combination of NIR scattering and extinction.

5.3.1 Models of dust emission

The fit residuals were mostly below 10% for all sub-millimetre bands and comparable to the observational uncertainties (Fig. B.1) and even the 160 μm extrapolated values were almost at the correct level. The fit quality was mostly independent of the tested dust models (Fig. B.2) and large residuals are not expected because of the large number of free parameters. However, because the radiation field scaling affects the whole map, differences in dust opacities (sub-millimetre vs. optical) could have caused systematic differences in the relative colour temperature of between dense and diffuse regions. There tended to be positive residuals at short wavelengths, partly correlated with the column density. These residuals suggest that the dust temperature is too low in the dense part of the cloud. In agreement with this possibility, the residuals were smaller in the W = 1.7 case when the cloud structure allowed more radiation to reach the cloud centre. The residuals could also reflect dust property variations that reduce the short-wavelength cloud opacity. However, the sub-millimetre data alone do not give strong constraints on dust models.

Although the sub-millimetre data can be fitted with any of the dust models, they do lead to significant differences in other parameters. Figure 12 shows the changes in the column density, J-band optical depth, and dust temperature relative to the J13 model.

In the case of AMMI, the central temperature is almost unchanged while in the outer parts there is an over 2 K drop relative to the J13 case. Although ice-coated grains are not be expected in diffuse regions, they can affect the temperature profiles deeper in the cloud. The dust temperature is an important parameter for cloud and core chemistry and even for their gravitational stability (Bergin et al. 2006; Sipilä et al. 2017). In these respects, even differences of 1 K can be significant.

The dust differences also affect the estimated mass distributions. Assumption of higher sub-millimetre emissivity would result in lower optical depths in the UV-optical-NIR regime, thus leading to a more uniform temperature distribution and smaller column-density variations. Such muted differences would clearly be noted if the comparison was made between models with the same radiation field. In our case, these effects are partly compensated by changes in the ISRF level. For AMMI dust, the column density is up to 75% lower than for the J13 model. The reduction is some 10% larger towards the cloud centre than in the outer cloud regions. Such trends are significant for the accuracy of mass estimates but also for estimates of the core density profiles. The spatial dust-property variations in real clouds naturally introduces additional uncertainties.

In spite of the lower model column densities, AMMI and THEMIS led to 40–100% increases in the predicted NIR opacity because both dusts have significantly higher opacities also at NIR wavelengths. This is illustrated by Fig. 13 that compares the J13 and AMMI opacity curves. Although the average increase in τ(J) is similar for AMMI and THEMIS (Fig. 12 and Table 3), the spatial distributions of the τ(J) increase are different. For AMMI the increase is smallest towards the cloud centre while for THEMIS it follows the dust abundance variations and peaks at the cloud centre. Such differences will be reflected in the contrast of the scattered-light images.

Figure 12 highlights two further factors affecting the emission modelling. If the cloud is more extended in the LOS direction, the optical depths become lower in the perpendicular direction, the central temperature rises, and the observed surface brightness is reproduced with a lower column density. When the LOS extent was increased by 70% (W = 1.7), the column density was up to ~20% lower and the core temperaturehigher by ~1 K. These effects are thus almost of similar magnitude as the differences between the dust models.

Because of the higher computational cost of the full treatment of stochastically heated grains, the long-wavelength emission was calculated assuming grains at an equilibrium temperature, which is true only for large grains. The absorption and MIR emission by smaller grains is energetically significant and when part of the emitted energy is transferred to shorter wavelengths, a higher column density is needed to produce a given sub-millimetre intensity. Figure 12 includes results for the J13 model (W = 1.0) when the full grain size distributions and the stochastic heating are taken into account. The column densities of the fitted model are higher by up to ~10% in the outercloud regions, where the MIR emission is strong. The difference decreases with column density and is only a couple of per cent towards the cloud centre. For the most accurate results, the full treatment of grain heating would still be preferred, if computationally feasible. For our models consisting of ~ 4 × 10⁶ cells, the full treatment slowed down the calculations by more than a factor of ten, to about one hour per iteration.

We assumed the Mathis et al. (1983) ISRF model as the reference and included the scaling k_ISRF as a free parameter. The spectral shape of the incoming radiation has additional second-order effects. The cloud models cover a limited volume, which can be assumed to be surrounded by outer cloud layers that attenuate the incoming radiation. By preferentially removing short-wavelength radiation, an external layer would increase the mean free path for the remaining radiation, resulting in smaller temperature gradients. The observations used in the modelling were similarly background-subtracted, to eliminate the extended foreground and background emission. Using these original surface brightness data with absolute zero points (Sect. 2.1) and the J13 dust properties, we estimate τ(J) = 0.52 for the average optical depth in the area used for background subtraction. If this corresponds to a layer between LDN 1642 (the modelled volume) and the stars providing the ISRF, the incoming radiation should be attenuated by e ^{−τ_Ext(ν)}, where the optical depth of the external layer could be up to τ(J) ≈ 0.26, half of the LOS value. We recomputed the J13 model with this change in the shape of the incoming ISRF radiation. The optical depth of the optimised model decreased only by 5%. The effect of τ_Ext is likely to be even smaller, since at least part of the LOS material is mixed with the stellar distribution. Therefore, the discrepancy in the τ(J) values (models vs. NICER) cannot be resolved by an external cloud layer or other similar changes in the shape of the ISRF spectrum.

Fig. 13 Dust optical depth for J13 and AMMI models relative to hydrogen column density N(H) = 1021 cm−2. The solid lines are for absorption and the dashed lines for scattering. The plot lists τ(250 μm)∕τ(J) ratios where the first number corresponds to absorption (points marked with crosses) and the second to extinction (marked with circles).

Fig. 14 Comparison of observed net surface brightness (dashed lines) and models with J13, AMMI, and THEMIS cases. The frames correspond to different model assumptions: the default parameters (frame a), cloud model with larger LOS extent (W = 1.7, frame b), four times larger intensity of the local radiation field (frame c), and the default model assuming no LOS sky background (Ibg = 0, frame d). Estimates Fig. 8 and include estimates for the areas marked in Fig. 8. The solid lines correspond to the larger area and the partly transparent lines to the smaller area. THEMIS results are shown only for the default case.

5.3.2 Models of scattered light

The morphology of the model predictions was consistent with the observations of the large-scale scattering but only if the background component I_bg(e^−τ − 1) was ignored. Better results were obtained by using model clouds with lower column density, in agreement with the direct NIR extinction measurements.

When we used model clouds obtained from the emission modelling, the NIR observations were underestimated most severely at the shortest wavelengths. The J-band surface brightness was mostly negative and the K_S-band and 3.4 μm signals were half those observed (Fig. 8). Similar wavelength dependencies existed for all dust models. This behaviour was suggestive of excessive NIR optical depths in the cloud models. In addition to the column density, the surface brightness depends on the dust model, the 3D shape of the cloud, and the sky brightness behind the cloud. We summarise these effects in Fig. 14, for the cloud models obtained from the fitting of dust emission.

AMMI dust increases the levels of scattered light but while the observed signal is underestimated at the shortest wavelengths, the 3.4 μm prediction istoo high. Apart from the level of the NIR signal, the change of the dust model does not improve the match to the observed shape of the SED. The larger LOS cloud extent clearly increases the predicted surface brightness and changes the SED shape in thecorrect direction. However, the effect remains too small.

The NIR modelling adopted the Mathis et al. (1983) ISRF values but, because the scattered light is directly proportional to the incoming radiation, the results could be easily rescaled. For example, the DIRBE measurements suggest some 40% higher NIR intensities (Lehtinen & Mattila 1996). As shown in Fig. 14c, we would need a much larger factor to match the general level of the NIR observations and, if the same scaling was applied to all bands, the SED shape would still not match the observations. The 1.25–3.4 μm radiation field had an effective colour temperature of 3500 K. Even if the field were dominated by 10 000 K sources, the J-band to 3.4 μm ratio would increase by less than a factor of three. In the absence of nearby massive stars, the discrepancy in the level and spectral shape of the NIR excess cannot be solved by modifications of the radiation field either.

If sky brightness behind the cloud were severely overestimated and we set I_bg to zero (Fig. 14d), the NIR signal would increase but this would still not fix the incorrect SED shape. If the ISRF is boosted to increase the scattering, the background term becomes relatively small and has only little effect on the NIR spectrum. Therefore, none of the above modifications can resolve the discrepancy between the observed and modelled NIR surface brightness. Of course, the independent NIR extinction measurements already suggested that the main problem resides in the high NIR optical depth.

A higher optical depth of an already optically thick cloud can directly decrease the intensity of the scattered light (Juvela et al. 2006), at the same time makingthe term I_bg(e^−τ − 1) more negative. Section 5.3.1 noted that the discrepancy between the NIR optical depths deduced from the dust emission models and the direct extinction measurements was between a factor of 2.3 and 4.6 (Table 3). When the calculations were repeated with cloud models with τ(J) decreased tomatch the NICER estimates, the results were clearly improved (Fig. 10).

Figure 15a summarises the results for lower-density cloud models when other parameters (I, W, and ISRF) are kept at their default values. The predicted spectra are now closer to the observations. The J13 model remains in the J band a factor of several below the observations. For the AMMI model both the intensity level and SED shape are much closer to the observations. THEMIS provides a higher average NIR brightness but its SED shape is more inconsistent with the observations. We emphasise that Fig. 15 shows the surface brightness excess relative to the reference area shown in Fig. 8. Therefore, the THEMIS results are sensitive to the density thresholds used in setting the relative abundances of the dust components. The reference area is mainly below the densities where the final transition to AMMI dust takes place. If this transition took place at a lower density, the THEMIS result would become more similar to the AMMI one.

In Fig. 15b, the ISRF is further assumed to be 50% higher and the LOS background is 50% lower. The first change would be in agreement with the DIRBE measurements, as discussed above. A 50% error in the I_bg is unlikely but the effect of I_bg is already relatively small, because of the lower optical depth and the higher intensity of the scattered light. With these changes, the J13 predictions increase significantly but remain below the observed values at the shortest wavelengths. AMMI matches the NIR data well, while the 3.4 μm value is slightly overestimated. For THEMIS, the MIR signal is clearly above the observations.

A further reduction in model column density should improve the match with the observed SED shape, increasing the short-wavelength signal relative to the longer wavelengths. However, the SED shape also depends on other factors, such as the details of the spatial variations of dust properties. We conclude that the LDN 1642 NIR observations can be explained by using models with evolved dust components, such as AMMI. The types of dust grains found in the diffuse ISM are not able to provide sufficient surface brightness or the correct spectral shape, not without improbable modifications to several parameters, including the intensity and the spectral shape of the external radiation field.

The modelling of the immediate environment of the sources B2 and B3 showed that the surface brightness is compatible with pure scattering (Sect. 4.3). The observations should have been sensitive to the potential additional NIR component of dust emission, unless that is restricted to within 10′′ (~ 0.01 pc) of the sources. The similarity of the surface brightness profiles between B2 and B3, and the good match with the scattered-light models also confirm that B3 is indeed part of the LDN 1642 cloud.

Fig. 15 Comparison of observed net surface brightness (dashed lines) and model predictions for J13 and AMMI dust. The density of the cloud models has been rescaled to match the NIR extinction measurements. In frame b, the general intensity of the ISRF has been scaled further by a factor of 1.5 and the intensity of the sky background Ibg by 0.5.

5.3.3 Discrepancy between emission and scattering

A major discrepancy existed between the higher NIR optical depths predicted by the dust emission modelling and, on the other hand, the lower values suggested both by the simulations of the NIR scattering and the direct extinction measurements with background stars. If the NIR optical depth could be reduced, the NIR surface brightness could be explained with dust with high albedos, such as AMMI.

The W = 1 and W = 1.7 models showed that the cloud shape does not provide a solution, possibly with the exception of the unlikely scenario of a very long filament viewed along the main axis. As an alternative to limit the dust temperature variations (and thus to decrease the column density), we also briefly tested the effects of a clumpy cloud structure by scaling the density values with Gaussian random fields with different powerlaw indices, with σ = 1 on logarithmic scale or with direct multiplication with N(μ = 1, σ = 0.35). However, the effects on the predicted τ(J) remained smaller thanthe difference between the W = 1 and W = 1.7 models.

Of the dust properties, the albedo and the scattering function are important for the NIR surface brightness but do not affect the main problem of the NIR extinction. The problem does not concern only the modelling but was already shown by the empirical result τ(250 μm)∕τ(J) ~ 10⁻³. This ratio is a factor of two lower for J13, τ(250 μm)∕τ_J = 0.49 × 10⁻³, and even lower for AMMI, 0.40 × 10⁻³.

The difference in the NIR extinction curve of the three dust models had a negligible effect on the τ(J) values (Table 3). Previous studies also have concluded the NIR extinction curve to be relatively constant, with variations at most at a 5% level (Indebetouw et al. 2005; Lombardi et al. 2006; Román-Zúñiga et al. 2007; Stead & Hoare 2009; Fritz et al. 2011; Ascenso et al. 2013; Wang & Jiang 2014). Some of these studies have targeted clouds with column densities higher than in LDN 1642. Our estimate E(H-K)∕E(J-H) = 0.73 ± 0.35 was fully consistent with the standard extinction curve, given its large uncertainty. The uncertainty caused by the shape of the extinction curve is thus likely to be below ~10%. Small-scale cloud structure could bias τ(J) values but in the other direction, reducing the extinction estimates (Lombardi 2009). Down to 40′′ scales, the small-scale structure was already taken into account with the help of Herschel observations. The ratio between the τ(J) values from the emission models and from the NICER calculations was spatially constant, which also suggests that errors related to cloud gradients or variations in the local stellar density are not significant.

Assuming that the observed τ(J) values are accurate and taking the observed τ(250 μm)∕τ(J) ratio as thestarting point, we tested ad hoc dust models where the opacities at λ > 60 μm were scaled with a constant to match τ(250 μm)∕τ(J) = 10⁻³. When thesemodified dusts were used in the emission modelling, the τ(J) values were reduced almost proportionally to the increase of the sub-millimetre opacity (Table 4). For J13, the τ(J) value of the emission model is nearly consistent with the NICER estimate while for AMMI there still remains a factor of two discrepancy. With the modified dust models, the radiation field estimates were increased to k_ISRF ~1.2–1.5, thus mainly between the lower Mathis et al. (1983) estimates and the higher values obtained from DIRBE observations (Lehtinen & Mattila 1996).

Figure 16 shows the NIR surface brightness predictions for these ad hoc dust models. Figure 16a can be compared to Fig. 15a, where, considering the NIR data, the only difference is in the NIR optical depth (with the factors in the fourth column of Table 4). The decreased NIR optical depth provided by the ad hoc dust models is sufficient to bring the NIR signal close to the observed level for the AMMI dust. The further reduction of the optical depth by a factor of ~2 in Fig. 15a does not lead to significant additional improvement, apart from the higher values in the J band. The J13 dust model is still excluded, however, even after the modifications. It would match NIR observations only with a much stronger radiation field, which would be in contradiction with the sub-millimetre emission modelling. For the THEMIS model, the short-wavelength intensities have increased, bringing them close to the observations. Comparisons with Fig. 14 and Fig. 15 show, however, that the increased brightness cannot be explained simply by the average cloud optical depth. Instead, it is partly caused by a change in the column density structure that has increased the NIR surface brightness contrast relative to the reference area. In the THEMIS case, the results depend on the densities at which the dust properties are assumed to change. If the density thresholds were lower by a factor of two, the contrast between the cloud centre and the reference area would increase, and the THEMIS spectrum would rise above the AMMI curve.

The modified THEMIS model had τ(J) close to the observed value (Table 4) and, together with the pure AMMI case, was closest to the observed NIR spectrum (Fig. 16). With further tuning of the relative abundances of the dust components, the NIR signal might be brought to an even better agreement with the observations. The modified THEMIS model was made explicitly consistent with the measured NIR vs. sub-millimetre optical depth ratio, τ(250 μm)∕τ(J) ~ 10⁻³. The high values of τ(250 μm)∕τ(J) (shown by direct observations and required by the modelling) also are qualitatively consistent with recent laboratory work. In these studies, interstellar dust analogues are shown to have much higher dust opacities at longer wavelengths (λ > 20 μm) compared to the current silicate dust models (Demyk et al. 2017). The laboratory measurements also reveal a significant temperature-dependence in the sub-millimetre opacities and opacity spectral indices that should be taken into account in future cloud models.

Fig. 16 Comparison of observed net surface brightness (dashed lines) and model predictions for dust models with ad hoc increase in the sub-millimetre emissivity. In frame b, ISRF has been scaled with a factor of 1.5.

5.4 Comparison to other studies

Juvela et al. (2012c) modelled the light scattering and dust emission of the northern filament of the Corona Australis cloud. Similar to the present study, the models fitting observations of sub-millimetre dust emission predicted NIR cloud opacities that were clearly higher than the direct extinction measurements with background stars. The observed level of NIR scattering could be matched only by assuming a significant increase in the intensity of the NIR radiation field. In that paper, models were calculated for two dust models, the R_V = 5.5 dust from Draine (2003) and the Ossenkopf & Henning (1994) dust with thin ice mantles. The latter was found to provide a better description of the Corona Australis filament, which contains a couple of dense clumps with column densities higher than in LDN 1642.

Ysard et al. (2013) modelled dust emission of the Taurus L1506 cloud, including density-dependent dust evolution. The strong sub-millimetre emission suggested the presence of dust aggregates at densities above 1500 cm⁻³. Ysard et al. (2016) carried out corresponding modelling of dust scattering using the THEMIS evolutionary dust model, comparing the results to NIR observations and the MIR Spitzer IRAC data on 21 starless cores. The coreshine observations required the presence of evolved dust, such as a combination of CMM and AMM (aggregates without ice mantles) or a combination of CMM and AMMI (aggregates with ice mantles), the cloudshine data being more compatible with the latter. The intensity and the balance between NIR and MIR brightness could be adjusted by changing the density thresholds for the transition between different dust populations. However, in that study the scattering was not modelled simultaneously with dust emission, relying on generic spherical density distributions instead. One of the main conclusions was that, thanks to the H-rich carbon mantles, the NIR–MIR scattering could be explained with a smaller increase in the grain volumes. In other earlier studies, coreshine was associated with grains larger than 1 μm (Andersen et al. 2013; Steinacker et al. 2015), but also was seen to be directly linked with the appearance of ice features (Andersen et al. 2014).

Togi et al. (2017) studied dust emission and NIR–MIR scattering in the cloud B207. The cloud hosts a single protostar and has a peak column density of N_H2 ~ 3.5 × 10²² cm⁻², which is three times higher than in LDN 1642. The analysis pointed to high dust albedo values that peak at A = 0.84 in the I band. The comparison with Ysard et al. (2016) models showed that observations could be explained best with dust properties similar to CMM+AMM or CMM+AMMI. While the models matched the NIR–MIR signal in the cloud core, they underestimated the sub-millimetre emission by a factor of two. The comparison is complicated because the peak column density is 2.5 times higher in B207 than in the model that it was compared with. If the model column density were scaled upwards, the sub-millimetre surface brightness would not increase proportionally, because of the simultaneous drop in dust temperatures. Furthermore, as seen in the LDN 1642 modelling, the predicted NIR intensities would be reduced, because of the reduction in the number of scattered photons and because of the larger negative contribution of the I_bg (e^−τ − 1) term. Therefore, also the B207 data seem to point towards the dust having a high opacity ratio between the sub-millimetre and NIR wavelengths.

In the present paper, we examined signs of dust evolution mostly by comparing the results for a diffuse-medium dust model, the model J13, and for a dense-medium dust model, the AMMI model for aggregates with ice mantles. In any realisticscenario, the dust properties should vary inside the cloud in a continuous fashion. Our NIR data covered only the central part of the LDN 1642 cloud and thus do not trace the full transition from diffuse medium to cloud cores. With optical data over a more extended area, Mattila et al. (2018) estimated in LDN 1642 an i-band albedo of A ~ 0.72 and showed those observations to be consistent with pure CMM dust. Saajasto et al. (2020) studied the thermal emission and NIR scattering in the cloud LDN 1521, also attempting self-consistent modelling of both FIR emission and NIR scattering. The best fitting models included two or three dust components and the dust evolution was modelled by modifying their relative abundances as a function of density. Compared to our results, the ratios (NIR optical depths in emission models vs. direct NIR extinction measurements) reported by Saajasto et al. were closer to unity, 1.56 and 0.74 in tests with the SIGMA (Lefèvre et al. 2019) and THEMIS dust models, respectively. The SIGMA model clearly overestimated the NIR surface brightness, while the THEMIS model predicted better the surface brightness in the dense parts of the cloud.

6 Conclusions

We have examined dust emission, scattering, and extinction in the high-latitude, star-forming molecular cloud LDN 1642. The new HAWK-I data provided estimates for the NIR extinction and net surface brightness, which is the sum of scattered light and attenuated LOS sky background. Together with the Herschel sub-millimetre maps, these data provided a good starting point for the testing of different dust models. The study led to the following conclusions:

• The maximum extinction in LDN 1642 is A_J = 2.5 mag, which corresponds to A_V = 9.3 mag (R_V = 3.1), at a resolution of 2′.

• There are no indications of NIR extinction-curve variations; the NIR colour excesses increase linearly with N(H₂) up to the highest column densities, and the observed ratio E(H-K)∕E(J-H) = 0.73 ± 0.35 is consistent with the standard extinction curve.

• We find an optical depth ratio of τ(250 μm)∕τ(J) ≈ 10⁻³. This result is similar to previous ratios found for cold clumps and a few times higher than in the diffuse medium, thus confirming the increase of the dust sub-millimetre emissivity.

• The sub-millimetre observations could be fitted well with radiative transfer models, irrespective of the assumed dust model. However, these result in tens of percent differences in the absolute N(H₂) values and the relative values between regions of low and high column density.

• Compared to the diffuse-medium dust model J13, the evolved dust model AMMI results in up to 2 K lower temperatures. This difference is a combined effect of changes in the sub-millimetre vs. optical opacity ratios and changes in the absolute opacity values.

• The models fitting the sub-millimetre dust emission predict NIR extinctions that are 2.3–4.6 times higher than the direct extinction measurements. The discrepancy affects all of the tested dust models and, in the modelling of the NIR surface brightness, results in SEDs with too low intensities and wrong spectral shapes.

• With dust properties appropriate for the diffuse medium (dust model J13), the modelled intensity of the scattered light was at least a factor of three below that of the observations. This difference remained true even if one assumed a 50% higher radiation field and a 50% lower sky brightness behind the cloud. This excludes J13 as a viable dust model for LDN 1642.

• We tested ad hoc variations of the dust models where the τ(250 μm)∕τ(J) ratio was increased to the empirically found value. The modified J13 dust model was still excluded because of weak NIR scattering. The modified AMMI and THEMIS models resulted in NIR–MIR signal almost at the observed level, with approximate agreement also in the ISRF scaling between the NIR and sub-millimetre models (Fig. 16b).

• The study shows that LDN 1642 contains evolved dust with high sub-millimetre opacity and NIR scattering cross section. The direct observations of the τ(250 μm)∕τ(J) ≈ 10⁻³ ratio and the modelling of dust emission and scattering show that the ratio of dust sub-millimetre and NIR dust opacitiesis higher than in the current dust models.

In this paper, we have examined dust properties using data on NIR extinction, NIR scattering, and sub-millimetre emission. Further crucial and complementary pieces of information may be provided by future observations, such as James Webb Space Observatory (Gardner et al. 2006) measurements of the MIR ice and silicate features, or observations with the planned SPICA satellite (Roelfsema et al. 2018) of dust polarisation and MIR-to-FIR dust spectra, clarifying the picture of both the large-grain properties and the populations of very small grains.

Acknowledgements

N.S. acknowledges the financial support from the visitor and mobility programme of the Finnish Centre for Astronomy with ESO (FINCA), funded by the Academy of Finland grant number 306531. E.M. is funded by the University of Helsinki doctoral school in particle physics and Universe sciences (PAPU). VMP acknowledges support by the Spanish MINECO under project AYA2017-88754-P.

Appendix A Comparison with DSS data

The digitised DSS images are in units of photon density. We use the Gaia observations to calibrate them to units of Jy sr⁻¹, using photometry from ~ 1000 stars from the DR2 catalogue. The selectedstars are within 1.5° of LDN 1642, but because the Gaia observations are considerably deeper than the DSS images and because of the non-linear response of the photographic plates (see right panel of Fig. A.1), we restrict the analysis to stars with intensities in the range 6500–46 000 e s⁻¹ in the red part of the G band filter (G_r). We convert the Gaia fluxes to physical units following the Gaia documentation (Busso et al. 2018, chapter 5.3.6), by first converting the flux to an AB magnitude (Oke & Gunn 1983). The Gaia instrumental magnitude is defined as G = −2.5logI + G_0,AB, where I is the weighted mean flux of the source and G_0,AB is the zero point in the AB system. The AB system can be generalised (Bessell & Murphy 2012) to be used with broad photometric bands so that AB = −2.5log⟨f_ν⟩− 56.10, where ⟨f_ν⟩ is the source mean flux per frequency and the constant 56.10 takes into account the fact that Gaia fluxes are in units of W m⁻² Hz⁻¹. Combining the two equations we have ⟨f_ν⟩ = I × 10^{−0.4 (G_0,AB +56.10)}. We use a value of 25.1161 for the G_0,AB, which is the zero point of the G_r band (Evans et al. 2018). We estimate the DSS photon density flux of the stars using aperture photometry with a fixed aperture size of 7.5′′. The conversion factors are estimated with a linear fit to the relation between the Gaia fluxes and the DSS photon densities, as shown in the left panel of Fig. A.1a. The right panel in Fig. A.1 shows the sensitivity of the DSS photographic plate, computed as averages over the sensitometric spots on the plate. The conversion factors are then used to convert the DSS image to Jy sr⁻¹. We assume an uncertainty of ± 20% in the regions where the photon density is in the range 4000–17 500 e s⁻¹. The conversion becomes uncertain outside of this range.

Fig. A.1 Fits to the DSS and Gaia observations. Left panel: a linear fit to the relation between the DSS and Gaia observationsand the best fit parameters of the fit. Right panel: response of the photographic plate from the sensitometer spots of the photographic plate. The red curve is a fifth order polynomial fit to the sensitometric data, with the coefficients listed in the figure.

We computed R-band predictions only for scattered light because we do not have estimates for the absolute sky brightness I_bg at this wavelength. Figure A.2 shows that the RT model predictions are one fourth of the observed sky brightness. The ratio shows a gradient, which could be an artefact from the DSS plate (Fig. 2). Figure A.2b shows that, after removing the mean gradient (some 1.5% per arcmin), the least-squares slope is even lower, . It may be biased towards lower values by residual contribution from point sources (high values of observed intensity). On the other hand, any surface brightness in the background sky would decrease the slope further. Therefore, the default models definitely fail to produce enough surface brightness in the R band.

Fig. A.2 Comparison of observed and modelled (J13, W = 1) R-band surface brightness. Frame a: map of the ratio with contours at equal to 20 kJy sr−1 and 40 kJy sr−1 (cf. Fig. 7a). The red circle indicates the reference area used for establishing a common zero point. Frame b: correlation as 2d histogram, with a logarithmic colour scale for the point density. stands for observations corrected for the main gradient. The dashed line shows the least-squares fit to data with kJy sr−1. The effect of background sky brightness is not included in the model.

Appendix B Further model calculations

B.1 Dust emission models

Figure B.1 compares the sub-millimetre observations to the modelling with J13 dust. The RT models were optimised to match the 250–500 μm data but the figure also shows a comparison with the 160 μm data.

Figure B.2 shows the fit residuals for alternative models with different LOS cloud extents (W = 1.0 and 1.7) and dust models J13, AMMI, and THEMIS. Except for the THEMIS model, the dust properties are constant throughout the model volume. The fits are found to be of similar quality, although with more variation in the 160 μm residuals.

B.2 Predictions of extended scattering

Figure 8 showed predictions of NIR surface brightness for the ISRF intensity of Mathis et al. (1983). Figure B.3 shows the same comparison when the cloud optical depths have been scaled down to match the average NICER extinction.

B.3 Scattering near embedded sources

Section 4.3 showed results for scattered light from spherically symmetric models of the source B2 and B3 environment, based on the use of the J13 dust model (Compiègne et al. 2011). Figure B.4 show the corresponding results for the AMMI dust models.

Fig. B.1 Sub-millimetre surface brightness in L1642. The observed 250, 350, and 500 μm maps (frames a–d) and the surface brightness predicted by the models with J13 dust with W = 1 (frames e–h). Frames i–l show the relative errors of the fits, .

Fig. B.2 Residuals of the modelled 250–500 μm surface brightness, . The columnscorrespond to 160, 250, 350, and 500 μm, respectively. The 160 μm data were not used in the model optimisation. The rows correspond to models with different dust (J13, AMMI, and THEMIS) and LOS cloud extent (W = 1.0 or W = 1.7). The averagevalue and the standard deviation of the relative residual is quoted at the bottom of each frame.

Fig. B.3 Comparison of 1.25–3.4 μm surface brightness between observations and models. The figure is the same as Fig. 8 but the model clouds have lower column densities that correspond to the NICER NIR extinction measurements. The colour scales are the same for all model plots of the same band. The yellow circles indicate the reference region used for background subtraction. The black contours are drawn at 1.5 times the maximum of the colour scale.

Fig. B.4 Surface brightness profiles in the vicinity of the embedded sources B2 (frames a–c) and B3 (frames d–f), for the J, H, and KS bands. This plot is the same as Fig. 11 but with the AMMI dust model (Ysard et al. 2016).

References

All Tables

All Figures

Fig. 1 Column density N(H2) (frame a), dust temperature Tdust (frame b), and J-band extinction AJ (frame c) of LDN 1642. The resolution of the maps is indicated by the circles in the lower left corner. The values of N(H2) and Tdust are based on MBB SED fits to Herschel SPIRE data, after background subtraction. The red circle in frame a shows the reference region used for the background estimation. Areas with low column density have been masked in the plots (white pixels). The labelled stars indicate the locations of embedded sources B1–B3.

In the text

Fig. 2 Large-scale environment of LDN 1642 shown by dust emission and optical scattered light. Left panel: Planck 857 GHz surface brightness, the white arrow indicating the direction towards the Galactic Plane. Centre panel: optical scattered light from DSS red channel. The white box shows the location of the LDN 1642 cloud and covers an area of 0.7° × 0.7°. Right panel: a close-up of the area marked with the white box in the centre panel, showing the extended R-band surface brightness. The black stars indicate the locations of the embedded sources B1, B2, and B3, and the red contours show the N(H2) column density. The lowest contour is at 15 % of the peak column density of 5.14 × 1021 cm−2, with the contour levels increasing in 15 % intervals.

In the text

Fig. 3 Correlations between the optical depths at 250 μm and in the J band. The left frame uses a NICER τ(J) map with 2′ resolution.The pixel values are plotted as black points. The blue points correspond to the same data but excluding pixels close to the background level. The dashed line is the least-squares line to the blue points and, for comparison, the solid line has a slope corresponding to the ratio of the averages of the background-subtracted values. In the right hand frame, we compare τ(J) estimates of individual stars to the Herschel column density estimates at 40′′ resolution. The parameters of the least-squares lines are given in the frame.

In the text

	Fig. 4 Estimated column density N(H2) as function of NIR colours J–H (in blue) and H–KS (in red). The parameters of the linear fits are given in the plot.
In the text

Fig. 5 Observed surface brightness in J, H, KS, and WISE 3.4 μm bands (frames a–d, respectively). Herschel column density and dust temperature maps are shown in frames e and f. Bright stars identified by DAOPHOT were removed and replaced with interpolated surface brightness, and faint stars have been eliminated with median filtering. In the WISE 3.4 μm band image (frame d), the white regions correspond to areas around bright stars that were excluded from subsequent analysis. The YSOs B1–B3, are identified in frame a. The yellow circle in frames a–d indicate the reference region used for setting a common zero level. The white box in frame e corresponds to the extent of the J, H, KS, and 3.4 μm maps (the area shown in frames a–d).

In the text

Fig. 6 Correlation between the observed scattered light in the J, H, KS, and 3.4 μm bands and Herschel column density. The data are at 25′′ resolution and sampled at 7′′ steps from the two areas marked in Fig. 8f, blue points for the larger and red points for the smaller area. The dashed lines and the dotted lines are robust linear least squares lines fitted to the red points and blue points, respectively.The fitted slopes and intercepts are quoted in the frames. The Spearman’s correlation coefficients are also given in the corresponding colours.

In the text

Fig. 7 Predicted intensity of scattered light. The cloud model is J13 with W = 1. Frames a–f: bands R, J, H, KS, 3.4, and 4.6 μm, respectively. Only the scattering of the external radiation field is included and the effects of attenuated sky background are not included. Frame e indicates areas used for correlation analysis, the larger area drawn with thin line (excluding the sources B1-B3) and the smaller area with thick line (see also Fig. 8).

In the text

Fig. 8 Comparison of observed (first row) and modelled (other rows) 1.25–3.4 μm surface brightness. The second row corresponds to J13 dust, W = 1, and no sky background (Ibg = 0). The effects of sky background are included in the other cases (J13 dust with W = 1 and W = 1.7, and AMMI with W = 1, and the THEMIS model with three dust components). Frame f indicates the areas used for correlation analysis, the larger area shown with a thin cyan line (excluding the sources B1–B3) and the smaller area with a thick cyan line. The yellow circles indicate the reference region used for setting a common zero level. Black contours are drawn at 1.5 times the maximum of the colour scale. The colour scale is the same for all models but different for observations.

In the text

	Fig. 9 Hydrogen column density N(H2) associated with each dust component in the THEMIS model with spatially varying dust abundances. The three frames correspond to the J13, CMM, and AMMI dust, respectively. The black circle shows the reference area used in scattered-light analysis (cf. Fig. 5).
In the text

Fig. 10 Correlation between the observed and modelled scattered light in the J, H, KS, and 3.4 μm bands. The model results correspond to models where the densities have been scaled down to match the NICER extinction measurements (J13 dust, and W = 1). The data have a resolution of 25′′ resolution and sample at 7′′ steps the areas marked in Fig. 8f. The blue points correspond to the larger area and the red points to the smaller area. The dashed lines show robust least squares fits to the red points. The values of the slopes and correlation coefficients r are given in the frames.

In the text

Fig. 11 NIR surface brightness profiles in the vicinity of the embedded sources B2 (frames a–c) and B3 (frames d–f). The black curvesshow the azimuthally averaged observed profile. The red dashed curves show the model predictions for J13 dust, consisting of the scattered light and the attenuated direct contribution of the central source. The average signal at 55–60′′ distance is subtracted from both curves. The cyan lines show the estimated direct contribution of the point source, based on the fitted source luminosity, the LOS extinction, and the psf shape.

In the text

Fig. 12 Comparison of model column density N(H2) (first row), NIR optical depth τ(J) (second row), and a cross section of the physical dust temperature T(dust) (third row) in alternative models. The first column shows the reference model with J13 dust and W = 1 LOS cloud extent. The other frames show differences to this model, for the dust and W values indicated on top. The last column (J13, SHG) is the same as the reference model but with calculations including the full treatment of stochastic heating. Temperature maps are not shown for the THEMIS and stochastically-heated-grain cases, wherethere is no single temperature per cell.

In the text

	Fig. 13 Dust optical depth for J13 and AMMI models relative to hydrogen column density N(H) = 1021 cm−2. The solid lines are for absorption and the dashed lines for scattering. The plot lists τ(250 μm)∕τ(J) ratios where the first number corresponds to absorption (points marked with crosses) and the second to extinction (marked with circles).
In the text

Fig. 14 Comparison of observed net surface brightness (dashed lines) and models with J13, AMMI, and THEMIS cases. The frames correspond to different model assumptions: the default parameters (frame a), cloud model with larger LOS extent (W = 1.7, frame b), four times larger intensity of the local radiation field (frame c), and the default model assuming no LOS sky background (Ibg = 0, frame d). Estimates Fig. 8 and include estimates for the areas marked in Fig. 8. The solid lines correspond to the larger area and the partly transparent lines to the smaller area. THEMIS results are shown only for the default case.

In the text

	Fig. 15 Comparison of observed net surface brightness (dashed lines) and model predictions for J13 and AMMI dust. The density of the cloud models has been rescaled to match the NIR extinction measurements. In frame b, the general intensity of the ISRF has been scaled further by a factor of 1.5 and the intensity of the sky background Ibg by 0.5.
In the text

	Fig. 16 Comparison of observed net surface brightness (dashed lines) and model predictions for dust models with ad hoc increase in the sub-millimetre emissivity. In frame b, ISRF has been scaled with a factor of 1.5.
In the text

	Fig. A.1 Fits to the DSS and Gaia observations. Left panel: a linear fit to the relation between the DSS and Gaia observationsand the best fit parameters of the fit. Right panel: response of the photographic plate from the sensitometer spots of the photographic plate. The red curve is a fifth order polynomial fit to the sensitometric data, with the coefficients listed in the figure.
In the text

Fig. A.2 Comparison of observed and modelled (J13, W = 1) R-band surface brightness. Frame a: map of the ratio with contours at equal to 20 kJy sr−1 and 40 kJy sr−1 (cf. Fig. 7a). The red circle indicates the reference area used for establishing a common zero point. Frame b: correlation as 2d histogram, with a logarithmic colour scale for the point density. stands for observations corrected for the main gradient. The dashed line shows the least-squares fit to data with kJy sr−1. The effect of background sky brightness is not included in the model.

In the text

	Fig. B.1 Sub-millimetre surface brightness in L1642. The observed 250, 350, and 500 μm maps (frames a–d) and the surface brightness predicted by the models with J13 dust with W = 1 (frames e–h). Frames i–l show the relative errors of the fits, .
In the text

Fig. B.2 Residuals of the modelled 250–500 μm surface brightness, . The columnscorrespond to 160, 250, 350, and 500 μm, respectively. The 160 μm data were not used in the model optimisation. The rows correspond to models with different dust (J13, AMMI, and THEMIS) and LOS cloud extent (W = 1.0 or W = 1.7). The averagevalue and the standard deviation of the relative residual is quoted at the bottom of each frame.

In the text

Fig. B.3 Comparison of 1.25–3.4 μm surface brightness between observations and models. The figure is the same as Fig. 8 but the model clouds have lower column densities that correspond to the NICER NIR extinction measurements. The colour scales are the same for all model plots of the same band. The yellow circles indicate the reference region used for background subtraction. The black contours are drawn at 1.5 times the maximum of the colour scale.

In the text

	Fig. B.4 Surface brightness profiles in the vicinity of the embedded sources B2 (frames a–c) and B3 (frames d–f), for the J, H, and KS bands. This plot is the same as Fig. 11 but with the AMMI dust model (Ysard et al. 2016).
In the text