© ESO 2021

1 Introduction

Understanding how stars are formed is one of the crucial questions in astronomy. The Herschel space observatory has provided us with detailed observations of nearby molecular clouds and shown that star forming regions have vastly diverse morphologies, from dynamically active filamentary fields to more quiescent clouds with simple geometries (Molinari et al. 2010; Men’shchikov et al. 2010; Juvela et al. 2012). These far-infrared (FIR) observations can be used to derive column density estimates and to study possible variations in dust properties. However, these studies are limited by our understanding of the emission properties of the grains, in particular the dust opacity and to a lesser degree the dust opacity spectral index.

The light scattered by dust grains at near-infrared (NIR) and mid-infrared (MIR) wavelengths can be studied and analysed without explicit assumptions on the FIR thermal emission properties of the grains, thus the scattering observations can be used to place additional constraints on dust properties and the density distribution. Lehtinen & Mattila (1996) were the first to study the extended surface brightness of dense interstellar cores at NIR wavelengths and later Foster & Goodman (2006) showed that it is possible to make large-scale maps of this extended surface brightness, naming it ‘Cloudshine’, and they attributed it to scattered light. Padoan et al. (2006) used radiative transfer modelling to show that the observed Cloudshine was consistent with the hypothesis of light scattering and could be used for studies of cloud structure. Star-forming clouds, but not the cores, are usually only moderately optically thick in the NIR and therefore there is a near-linear dependence between the surface brightness and the column density, assuming that the dust properties do not change with column density.

At MIR wavelengths, especially in the Spitzer 3.6 and 4.5 μm bands, a surprisingly high surface brightness was detected towards several cloud cores (Steinacker et al. 2010; Pagani et al. 2010; Juvela et al. 2012). Explaining the higher-than-expected surface brightnesses with the classical grain size distributions (Mathis et al. 1977) proved difficult, implying the presence of larger grains with sizes of the order of ~ 1 μm (Steinacker et al. 2010, 2015; Pagani et al. 2010; Andersen et al. 2013; Lefèvre et al. 2014). The high surface brightness towards cloud centres was named ‘Coreshine’ and is considered to be a direct evidence of dust growth in dense clouds.

The comparison of thermal dust emission and light scattering provides crucial information on the dust properties. On the other hand, because of limited resolution, or issues caused by anisotropic illumination and line-of-sight confusion, careful radiative transfer modelling is often needed to place constraints on theoretical models of dust.

We have chosen the cloud LDN1512 (hereafter L1512) to study and model NIR light scattering and thermal dust emission at FIR wavelengths. The cloud has a simple cometary morphology, is nearby (140 ± 20 pc, Launhardt et al. 2013), and based on the reported small line width of N₂ H⁺ (Caselli et al. 1995, 2002; Lin et al. 2020), appears to be quiescent.

The cloud has been previously mapped with the Herschel space observatory using both the photodetecting array camera and spectrometer (PACS; Poglitsch et al. 2010) and the spectral and photometric imaging receiver (SPIRE; Griffin et al. 2010) instruments. As discussed by Launhardt et al. (2013) and Lippok et al. (2013) the Herschel observations show a single starless core surrounded by a diffuse envelope. Based on the fitting of the ¹³CO line, Lippok et al. (2013) showed that the envelope of the core has a higher gas temperature compared to the central regions and their stability analysis shows that the core is thermally supercritical. The N₂ H⁺ observations of Lin et al. (2020) show a low kinetic temperature of ~8 ± 1 K within the innermost ~0.017 pc of the core, and based on their chemical modelling, the cloud is sufficiently evolved that the N₂ chemistry has reached a steady state. Their results suggest that L1512 is likely older than 1.4 Myr and that ambipolar diffusion has led to the formation of the core.

We use J, H, and K_S bands images, shown as a three-colour image in Fig. 1, from the Wide InfraRed CAMera (WIRCAM) on the Canada-France-Hawaii Telescope (CFHT). The observations are presented in Lin et al. (2020), where only the stellar content of the images has been exploited. in this paper we concentrate on the extended emission. After background subtraction, all three bands show a clear extended surface brightness component.

In this study, we simultaneously model the cloud at NIR and FIR wavelengths using an anisotropic radiation field and a cloud model derived from the Herschel observations together with dust models for the diffuse interstellar medium (ISM) dust, based on the model described by Compiègne et al. (2011) and three models that take into account the evolution of dust grains in shape and size. The aim of our study is to find a solution, or place constraints, for the cloud, the radiation field and the properties of the dust grains, that would give a consistent explanation for the observed extinction, NIR scattering, and dust emission.

This paper is organised as follows. In Sect. 2, we give an overview of the archival observations used in the paper and describe our NIR observations. In Sect. 3, we describe our radiative transfer methodology and explain our radiation field and cloud models. In Sect. 4 we present our results and discuss our findings in Sect. 5. Finally, in Sect. 6, we summarise our findings and provide our conclusions.

2 Observations

The Herschel observations were downloaded from the Herschel science archive and are a part of the Herschel key project The Earliest Phases of Star Formation (EPoS, PI: O. Krause). The EPoS project used both the PACS and the SPIRE instruments to cover a total of 12 different fields between wavelengths from 100 to 500 μm. In the projects source list, the source CB 27 corresponds to the cloud L1512. The nominal scanning speeds for the PACS and SPIRE instruments were set to 20′′ s⁻¹ and 30′′ s⁻¹, respectively. The general noise level in the EPoS projects maps is ~6 mJy beam⁻¹, but for the L1512 the noise levels are slightly higher ~11 mJy beam⁻¹. A more detailed description on the data reduction is given by Launhardt et al. (2013).

The NIR observations cover the wide filters J (1.25 μm), H (1.6 μm), and K_S (2.15 μm) (see Fig. D.1) and were obtained in 2013 using the WIRCAM instrument on the CFHT. The observations were carried out using the Sky-Target-Sky (STS) dithering mode, to subtract the atmospheric IR-emission and to preserve any extended scattered light from the source. The seeing conditions during the observations were good, with typical values less than 1′′. Data reduction was performed at the TERAPIX center using a specific reduction method to recover the extended emission.

3 Methods

Our aim is to simultaneously model the cloud in NIR and FIR and to compare our modelling results with observations. In this section we describe our model of the density distribution within the cloud and our radiation field model.

3.1 Cloud model

We use a three-dimensional model cloud discretised onto a grid of 225³ cells with a Gaussian density distribution along the line-of-sight (LOS). The angular cell size of our model cloud is 6′′ which corresponds to the pixel size of the Herschel SPIRE 250 μm map. Thus, the total angular extent of our model is 22.5′ × 22.5′, which at a distance of 140 pc corresponds to a physical size of ~0.4 × 0.4 pc. The density distribution is based on the Herschel SPIRE 350 μm surface brightness (determining the column density distribution on the plane of the sky) and scaled along the LOS (z-coordinate) by a Gaussian distribution $$\rho (z)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{{(z-\mu )}^2}{2{\sigma }^2}\right),$$(1)

where μ = 0.0 pc and with either σ = 0.07 pc, labelled Narrow, or with σ = 0.11 pc, labelled Wide. The column density distribution is optimised during the model fitting, taking into account the temperature structure of the model cloud.

3.2 Radiation field

The radiation field used in our modelling is based on the Diffuse Infrared Background Experiment (DIRBE¹) Zodi-Subtracted Mission Average (ZSMA) maps (Hauser et al. 1998) at 1.2–240 μm. Outside the DIRBE wavelengths we followed the Mathis et al. (1983) model, scaling the values below 1 μm by 1.4 to match the level of the DIRBE data. Similarly, we used linear interpolation in the range from 240 to 650 μm, to avoid a discontinuity in our model (Fig. 2). The intensity was distributed over the map pixels following the surface brightness distribution of the DIRBE maps. The spatial distribution at wavelengths shorter than 1 μm follows the distribution of the DIRBE J band (Fig. 3). At wavelengths longer than 240 μm, we assume the spatial distribution of the DIRBE 240 μm band. Compared to the Mathis et al. (1983) model, our radiation field has significantly higher intensity in the range [10, 240] μm, but the contribution from these wavelengths to the heating of the dust particles is not significant. However, the increased radiation field strength at NIR and MIR wavelengths will affect our light scattering modelling.

Fig. 1 Three colour image of the field L1512. The colours correspond to the NIR surface brightness at J (blue colour), H (green), and KS bands (red).All point sources have been masked.

Fig. 2 Intensity of the radiation field, green line, as a function of wavelength. The purple line shows the Mathis et al. (1983) model. The black dots show the intensity of the DIRBE observations averaged over the HEALPIX map.

Fig. 3 Model for the all-sky sky brightness distribution in the J band. The map has been rotated so that the direction towards L1512 is in the centre. The intensity of the model has been truncated at 20 MJy sr−1 for plotting and the image has been smoothed by a Gaussian with FWHM = 2°.

3.3 Radiative transfer

To solve the dust emission and the surface brightness of the scattered radiation, we used the Monte Carlo radiative transfer programme SOC (Juvela 2019). A schematic overview of our modelling is shown in Fig. E.12.

The dust grain properties are defined by the absorption and scattering efficiencies, Q_abs and Q_sca, and a scattering function described using the asymmetry parameters of the Henyey-Greenstein scattering functions, g, summed over the individual dust components. The default parameters are taken from the Compiègne et al. (2011) model for the diffuse ISM. Variations in these parameters, as well as three models including dust evolution, are also considered. A more detailed description of the dust models is provided in Appendix A. The absorption and scattering properties of the dust grains vary significantly between the models, an example of the optical properties of selected models is shown in Fig. 4.

As a first step, the dust properties, the radiation field, and cloud model are used to compute the dust emission at the Herschel SPIRE wavelengths of 250, 350, and 500 μm. We fit our model to the observed dust emission by an iterative process where the density of the model cloud and the strength of the radiation field are fitted to the observations. For each iteration step, we compare the difference between the observed and simulated surface brightness at 350 μm, I_obs (350) − I_sim(350) and, for each map pixel, scale the density of each model cloud cell along the corresponding LOS by $$K_{\rho }=\frac{I_{\text{obs}}(350)}{I_{\text{sim}}(350)}.$$(2)

Similarly, we scale the radiation field by an average scalar factor $$K_{\text{ISRF}}=\langle \frac{I_{\text{obs}}(250\mu \text{m})\text{/}{I}_{\text{obs}}(\text{500}\mu \text{m})}{I_{\text{sim}}(250\mu \text{m})\text{/}{I}_{\text{sim}}(\text{500}\mu \text{m})}\rangle .$$(3)

In the K_ISRF computation, we use only the pixels in the central part of the map, within the red circle shown in Fig. 5. The scaled radiation field and cloud model are then used to compute a new estimate for the dust emission.

Once the emission fit has converged, we use the final scaled cloud model and radiation field in a separate radiative transfer calculation to predict the scattered light at NIR wavelengths.

The observed surface brightness excess is relative to the background sky $$I_{\nu }^{\Delta }=I_{\text{sca}}+I_{\text{bg}}\times (e^{-\tau }-1),$$(4)

where I_sca is the intensity of the scattered light, I_bg is the absolute value of the background radiation seen towards the cloud and τ is the optical depth of the cloud. To compare the observations and simulations a background needs to be subtracted from the simulated maps. For the J and K_S bands, the background is estimated from the DIRBE ZSMA maps as an average over the four closest pixels to L1512 and for the H band we estimate the value from the J and K_S band values following the Mathis et al. (1983) model. However, the DIRBE observations include the contribution from point sources which for Eq. (4) needs to be removed. We estimate this contribution from the 2MASS point source catalogue by integrating the intensity of all point sources that are within the DIRBE pixels (see Lefèvre et al. (2014) for details). The integrated intensity values are then subtracted from the average DIRBE values resulting in background levels of 0.059, 0.061, and 0.040 MJy sr⁻¹ for the J, H, and K_S bands, respectively. The pixel-to-pixel variations in the DIRBE observations are at most ~ 15% for the J and K_S bands. Because of the relative proximity of the cloud, D = 140 pc, we assume that all of the background radiation is emitted from behind the cloud (i.e. there is no foreground emission).

Fig. 4 Comparison between the optical depth per H as a function of frequency for the Default, LGM, SIGMA, and THEMIS models. The yellow shaded regions show the frequency range from 250 to 500 μm and from 1.2 to 2.5 μm.

4 Results

In this section we report the main results of our study. In Sect. 4.1, we analyse the NIR and Herschel observations and describe the results of our modelling in Sect. 4.2.

4.1 Observations

Figure 5 shows the column density map obtained with modified blackbody (MBB) fits of the Herschel SPIRE data. We assume a spectral index β = 1.8, which corresponds to the average in nearby molecular clouds (Planck Collaboration XXV 2011). The SPIRE observations were convolved to a common resolution of40′′ and colour corrected using the factors of Sadavoy et al. (2013). The J band surface brightness has a local minimum at the location of the column density maximum, showing that the J-band scattering has saturated because of high optical depth. The dust emission peaks south of the column density peak in all the SPIRE channels, indicating stronger heating from that direction.

The Gaia data release 2 (DR2, Gaia Collaboration 2018) provides the parallaxes and magnitudes of over 1.6 billion sources, a catalogue that can be used to estimate distances within the Galaxy. To locate stars that could produce additional illumination, we select all stars that are within two degrees of the cloud centre and brighter than 15^mag in the Gaia G band. We convert the parallaxes to distance estimates assuming $r=1/\overline{\omega }$, where $\overline{\omega }$ is given in seconds of arc (see discussion by Bailer-Jones 2015 and Luri et al. 2018 for the accuracy of this method). We find nine stars that are likely to be closer than 100 pc to the cloud and have Gaia G band magnitude brighter than 10 magnitudes. Only three of these stars are on the southern side of the cloud (Fig. 5). However, their estimated G band illumination compared to the observed NIR surface brightness is I_stars∕I_SB = 0.001, with I_stars = 5.6 × 10⁻⁵ MJy sr⁻¹ and assuming an average surface brightness of 0.07 MJy sr⁻¹. Thus, compared to the illumination from the ISRF, the contribution from these stars to the illumination of the cloud is minimal.

We derive the J band optical depth from the WIRCam observations using the Near-Infrared Color Excess Revisited (NICER; Lombardi & Alves 2001) method and assume a standard extinction curve (Cardelli et al. 1989). The submillimetre optical depth was obtained by fitting the spectral energy distribution (SED) of the Herschel SPIRE observations with MBB curves, assuming an opacity spectral index β = 1.8 and a dust opacity of ${\kappa }_{\nu }=0.1{(\nu /1000\text{GHz})}^{\beta }$ cm² g⁻¹ (Beckwith et al. 1990), thus $${\tau }_{250}=\frac{I_{250}}{B_{250}(T)},$$(5)

where I₂₅₀ is the fitted 250 μm intensity, B₂₅₀ is the Planck function at 250 μm, and T is the colour temperature from the SED fit.

The maps for J-band extinction and τ₂₅₀ are shown in Fig. 6 and the correlation between the two quantities in Fig. 7. In Fig. 7, we have excluded the central region where the extinction estimates are uncertain due to the low number of background stars (the region within the contours in Fig. 6). A linear least squares fit gives τ₂₅₀ ∕A_J = (4.9 ± 0.3) × 10⁻⁴ mag⁻¹. Masking the high values above τ₂₅₀ > 0.002, reduces the ratio to τ₂₅₀ ∕A_J = (2.3 ± 0.3) × 10⁻⁴ mag⁻¹. The ratio of τ₂₅₀ ∕A_J = 4.9 × 10⁻⁴ mag⁻¹ is slightly higher than what one would expect for diffuse clouds (Bohlin et al. 1978; Planck Collaboration XI 2014), but is lower than the average value in dense clumps found in Juvela et al. (2015).

Fig. 5 Contours of J-band surface brightness on Herschel column density map. The plus signs show the locations of the 250, 350, and 500 μm emission maxima (red, green, and black, respectively). The white arrows indicate the projected direction of the three brightest nearby stars and the numbers next to the arrows indicate the 3D distance to the star and the Gaia G band average magnitude. The red circle shows the area that is used to compute the KISRF scaling factor.

4.2 Radiative transfer models

In this section we describe the results of our emission and scattering modelling. We use the Compiègne et al. (2011) dust model as the default model, but will also test variations, for example the effects of different size distributions and grain optical properties. We refer to these as ‘COM models’.

Table 1 lists the dust models used, and a more detailed explanations are provided in Appendix A. We computed for each model the column density, background-subtracted intensity in the J, H, and K_S bands, the J band and 250 μm optical depths (τ_J,em and τ₂₅₀), and the scaling factor of the radiation field K_ISRF. The values in Table 1 are averages over 5 × 5 map pixels centred on point 1. The temperature of the core is computed as an average value over 5³ cells, centred at the core. The values derived from observations are all averages over 5 × 5 pixels centred on point 1. The core temperature is based on the N₂H⁺ observationsby Lin et al. (2020). In addition to the models in Table 1, we tested three further modifications to the dust properties. The differences to the Default model were minor, and these results are presented only in Appendix C.

Our radiative transfer computations do not include stochastically heated grains (SHG), as the emission from these grains is minimal in the SPIRE bands and including the stochastic heating would increase the computations times significantly. We computed stochastic heating only for two test cases, for models that were previously fitted assuming dust at an equilibrium temperature. The results are shown in Figs. B.1and B.2. Including the stochastic heating decreases the emission in the SPIRE bands by ~ 10−15%, because some of the energy is now emitted in the MIR wavelengths. The emission is at 100 μm 30–45% and at 160 μm 15% above the observed values. For the Default model, the morphology of the simulated 100 and 160 μm maps agrees with the observed maps. For the THEMIS model, the bright rim in both maps is up to 40 MJy sr⁻¹ brighter than observed. For these two cases, we also computed the predicted SEDs for the full wavelength range 3.6–500 μm (Fig. B.3).

Fig. 6 J band extinction map and the optical depth at 250 μm. The white contours in the J band extinction map show the estimated uncertainties of 0.21, 0.3, 0.41 mag, from lowest to highest contour level.

Fig. 7 Correlation between the optical depth at 250 μm and the J band extinction. The purple line is a least squares fit to the data and the Pearson correlation coefficient is 0.65 (p ≪ 0.01).

4.2.1 COM models: dust emission

Figure 8 compares the dust emission in the Default model with the observations. The fit residuals (observed value minus the model prediction) of the simulated 350 μm emission are ± 1.5%. For the 250 and 500 μm bands, the residuals are on average ± 3%, but increase up to − 10% in the densest region.

The fitted radiation field strength is lower than our reference model for all test cases, but on average, the observations can be fitted with a scaling factor of K_ISRF = 0.5. The lowest value is K_ISRF ~ 0.44 for the model LGM, while the highest value is K_ISRF ~ 0.66 for the Scaled4 model. Thus, increasing the emissivity of the grains, models Scaled2 and Scaled4, increases the required K_ISRF to compensate for the grains being colder for a given radiation field. The τ₂₅₀ values of the models are within higher but within 15% of the values derived from the observations with SED fits. In models with larger grain sizes, such as LG and LGM, the difference increases to ~ 30%.

The column density of the simulated Default case has a maximum of ~ 2.6 × 10²² cm⁻², whereas themaximum value derived from the Herschel observations via MBB analysis is 1.2 × 10²² cm⁻². Changes in the assumed emissivity of the dust grains naturally affect the derived column density, see also the discussion by Malinen et al. (2011) and Ysard et al. (2012) on the uncertainties of the MBB analysis.

A low kinetic temperature of ~8 ± 1 K within the innermost ~0.017 pc of the core was derived by Lin et al. (2020) using N₂H⁺ line observations. We computed for each model an average temperature T_core over a cube of 5³ cells centred on the core. For the Default model T_core = 9.34 K, which is ~1 K higher than the N₂H⁺ estimate. The other models also show temperatures 1–2.5 K above the N₂ H⁺ estimate. The models with higher emissivity, Scaled2 and Scaled4, and with larger grains, LG and LGM, are 0.7−1.5 K warmer than the Default model. This is caused by the higher submillimetre emissivity leading to lower column density and lower cloud optical depth at the short wavelengths responsible for dust heating. This, in turn, results in higher dust temperatures in the core. On theother hand, the increased LOS illumination, model Wide, only increases the core temperature marginally to 9.7 K.

4.2.2 COM models: light scattering

We calculated the scattered surface brightness for the J, H, and K_S bands using the density distribution and the radiation field obtained from the fits to the dust emission and by adding net effect of the background I_BG × (e^−τ − 1), resulting in simulated surface brightness maps.

The computed surface brightness maps are shown in Fig. 9 and the individual components of the signal in Fig. E.5. The modelled surface brightnesses are up to a factor of four lower compared to the observed surface brightness. Furthermore, the general morphology of the simulated maps does not match the observed maps. The faint striations (see Fig. 9E) are more pronounced and the morphology of the bright rim is narrower in the simulated maps. The central part of the cloud has a low intensity, ~ 0.01 MJy sr⁻¹, but the low-surface-brightness area is more extended. In the simulated H and K_S band maps, the rim and the striations are clearly visible, unlike in the observations. As in the case of the J band, the morphology of the core does not match the observations. In all simulated maps, the Veil is considerably more prominent than in the observations. Since the Veil cannot be seen in the WIRCam observations but is clearly seen in the Herschel observations, it is possibly not connected to the L1512 but is further away on the same line-of-sight.

The increase in emissivity, and thus lower column density, of the Scaled2 and Scaled4 models has increased the surface brightness by up to 40% in the J band. In the H and K_S bands the intensity has decreased by 10% in the diffuse regions and increased by up to 30% in the dense regions. This produces a more compact core and the Veil is less prominent. However, the intensities are still a factor of two to three below the observed values. For the models that include larger grains, LG and LGM, the surface brightness values are only ~ 30% lower in the H and K_s bands and ~20% higher in the J band than the observed values. However, the morphology of the maps does not match the observations: the core is less compact and the Striations and the Veil more prominent.

The values of τ_J,em are on average three times higher than the NICER estimates (see Table 1). The model Scaled4 is an exception, the optical depth τ_J,em = 1.55 being within 5% of the value derived from observations.

Fig. 8 Observed (first row) and simulated (second row) emission maps using the Default model at 250, 350, and 500 μm. Third row: difference between the observations and simulations.

Fig. 9 Observed (first row) and estimated surface brightness maps from the Default model at the J, H, and KS bands (second row). All maps have been background substracted.

4.2.3 Models with dust evolution

The Defaultmodel is designed for diffuse lines of sight, but in dense cores dust grains are expected to grow due tocoagulation and mantle formation. Thus, we also test dust models that take into account the evolution of dust grains. A comparison between the optical properties of the models is shown in Fig. 4.

We test three cases where the dust evolution is modelled by changing the properties of the dust grains with increasing density. A more detailed explanation of the models is provided in Appendix A. As a summary: the model DDust is a combination of models Default and LG, where the relative abundance of the LG component increases with increasing density. The SIGMA model uses the Default model for diffuse gas and a combination of aggregate silicate and carbon grains with ice mantles in the dense regions. Finally, the THEMIS model uses a diffuse component with core-mantle (CM) grains, an intermediate component with core-mantle-mantle grains (CMM), and a dense dust component with amorphous core-mantle-mantle aggregates that have ice mantles (AMMI) (Köhler et al. 2015; Ysard et al. 2016). The evolution of the dust grains, with increasing density, will affect all dust parameters. Thus, one cannot identify a single parameter or parameters that change, as not only the physicalparameters, emissivity absorption and scattering properties, and size distribution, but also chemical properties evolve due to formation of aggregate grains and formation of (ice)mantles. Furthermore, the changes in the chemical composition are also expected to effect the optical properties, thus the affect of dust evolution is intricate.

In the fits of the dust emission, the 250 and 500 μm residuals in the central region are smaller for SIGMA (Fig. E.10) than for the default dust model (e.g. Fig. 8). The residuals are slightly larger in the case of the THEMIS model (Fig. E.10) and in the case of the DDust (Fig. E.10) model are similar to the Default model.

Compared to the Default case, the SIGMA and DDust models have higher K_ISRF factor, 0.67 and 0.75 respectively,while the K_ISRF factor of the THEMIS model is 0.47. As with the models that include larger grains, LG and LGM, τ₂₅₀ of the SIGMA, THEMIS, and DDust models are higher by ~30% compared to the Default model. The NIR optical depths are τ_J = 3.4, 1.05, and 14.4 for the SIGMA, THEMIS, and DDust models, respectively. The corresponding core temperatures read from the 3D models are 10.6, 7.6, and 9.3 K. Thus, in spite of the lowest τ_J value (even below the NICER measurement), THEMIS results in the lowest dust temperatures, some 0.5 K below the estimate derived from N₂ H⁺ line observations.

The THEMIS, and DDust models produce column densities close to the Default model with N(H₂) = 1.4 × 10²² cm⁻² and 1.8 × 10²² cm⁻², respectively. The column density of the SIGMA model is lower, with N(H₂) = 9.61 × 10²¹ cm⁻². The model column density profiles (Fig. 10) agree with the profile derived from the Herschel observations via MBB fits, but are narrower than the FWHM = 2.0′ estimated from the Herschel estimate. The profile of the Default model, FWHM = 1.4′, is close to the modelled N₂ H⁺ profile with FWHM = 1.32′ (Lin et al. 2020).

THEMIS produces two to three times higher NIR surface brightness than the Default model (Fig. E.4), but the SED shape and the morphology still do not match the observations. Although the central region of the cloud in H and K_S bands has surface brightness values within 0.02 MJy sr⁻¹ of the observed value, the more diffuse regions are brighter than observed regions. The striations are clearly visible in the THEMIS maps (H and K_S bands), whereas in observations no extended surface brightness is seen. The model SIGMA has the highest surface brightness of all our test cases, the intensities in the diffuse regions are in excess of 0.1 MJy sr⁻¹ for the H and K_S bands and even for the J band the values are in range [0.04,0.12] MJy sr⁻¹. The morphology is close to the observed morphology, with a compact core in the H and K_S bands and the surface brightness of the Veil is lower. The model DDust shows a clear dip in the core, even in the K_S band map, a factor of two to three lower surface brightness compared to the observations and the model clearly overestimates the surface brightness of the diffuse regions.

Fig. 10 Column density of H2 derived from Herschel observations and the comparison between column density profiles for selected models. The profiles have been computed as average values over 3 × 3 pixels along the arrow in the left panel.

Fig. 11 Observed map of the J band scattered light (left panel), the red numbers from one to four indicate the locations from where we have extracted NIR spectra, shown in the right panel. The solid lines show the observed spectra and the dashed lines show the spectra derived from the Default model. The dotted lines show the simulated spectra derived from the Default model, but the intensity of the scattered light has been multiplied by a factor of three before background subtraction.

4.3 NIR pectra

Figure 11 shows NIR spectra for observations and the Default model, for the four positions indicated in the figure. Point 1 corresponds to the brightest part of the J band map and the three other points tracedensity variations across the densest part of the cloud. In addition to the low intensities, the general shapes of the simulated spectra do not match the observations, the relative brightness of the K_S and J bands being too high. A factor of three increase in the scattered light would increase the net surface brightness by a larger factorbut would not improve the match with the SED shape. However, the observed H-band value has somewhat higher uncertainty because of the lack of a direct background sky brightness measurements at that wavelength.

Figure 12 shows the NIR spectra for all test cases and points 1, 3, and 4. Point 2 is similar to point 4 and its spectra are shown in Fig. E.11).

The error bars assume an uncertainty in the background sky brightness that is 30% in the H and 20% in the other bands. The Default-case simulated intensities are a factor two to three lower than the observed values. The models LG (Default model with grains up to 5 μm in size) and LGM (As LG but relative amount of large grains increased) produce higher intensities, but the intensity of the H band is ~ 30% lower than the observed value while the J and K_S bands are 30 to 40% brighter. The SIGMA model (two dust components, Default for diffuse and a component derived with SIGMA for dense LOS) is clearly overestimating the intensity of the scattered light in all three channels. In point 1, the THEMIS model (three dust component based on the THEMIS framework) is closest to the observations, but the intensity of the K_S band is 0.05 MJy sr⁻¹ too bright. However, in the point 3, the intensity and the shape of the spectra of the THEMIS model is within 10% of the observed intensity. In general the spectra from point 3 is more easily reproduced and the models LG, LGM, and DDust all produce approximately the correct SED shape but do not match the level of intensity. Like in the point 1, the SIGMA model is overestimating the intensity, but the shape of the spectra is correct.

For point 4, the COM models tend to underestimate the NIR intensities, with the exception of the LG and LGM models that overestimate the J and K_S bands, but underestimate the H band. However, models Scaled2 and Albedo give very similar results.

The increased FIR and submillimetre emissivity in the cases Scaled2 (emissivity for λ > 60 scaled by a factor of two) and Scaled4 (as Scaled 2 but scaled by a factor of four) decreases the column density, thus, the higher intensity of the H band compared tothe K_S band can beunderstood as a saturation effect. However, increasing the emissivity further does not improve the match as the column density becomes too low and the intensity of the scattered light is reduced. The increased emissivity of the grains decreases the J band optical depth, which for model Scaled4 is close to the NICER estimate, but the morphology of the surface brightness maps and the shape of the NIR spectra are further away from the observations. Thus, emissivity alone can not reconcile observations and simulations.

Fig. 12 J, H, and KS band intensities from observations (horizontal lines) and different models (symbols) for the map point 1 (lower panel), point 3 (middle panel), and point 4 (upper panel). The colours correspond to the J (purple), H (red), and KS (black) bands. All intensity values have been background subtracted. We assume a 20% uncertainty in background sky estimates for the J and KS bands and an uncertainty of 30% for the H band.

5 Discussion

Based on our modelling of the LDN1512 observations, its thermal dust emission can be fitted with many different assumptions on dust properties. However, we can not simultaneously fit both the dust emission and NIR scattered light with the COM models. The Default model was designed for high-latitude diffuse lines of sight, while the L1512 cloud has a dense central core. The disparity between the observations and our radiative transfer models based on the Default model is a clear indication for the need of evolution of the dust grains and we thus tested additional models, DDust, SIGMA, and THEMIS. We are able to reproduce the observed intensity of dust emission and scattered light at NIR wavelengths with the THEMIS model. However, the morphology of the surface brightness maps shows considerable deviation from the observations. In this section we discuss these results in more detail.

5.1 Uncertainties of the radiation field

The ISRF used in our simulations is based on the DIRBE observations and the shape of the Mathis et al. (1977) model. We did not test changes in the anisotropy of the external field. However, based on Gaia point sources, the contribution from the nearby stars is minimal.

Based on the intensity variation between the DIRBE pixels, the background uncertainty is of the order of ~ 20%, but for the H band, for which there are no direct DIRBE measurements, we have assumed a value of 30%. With the exception of the model SIGMA (model with two dust components, Default dust for diffuse regions and dust created with SIGMA for dense LOS), the spectra derived from both points 1 and 3 are systematically underestimating the intensity of the H band. For the THEMIS model (model with three dust components based on the THEMIS framework), the intensity of the J and K_S bands from point 3 are comparable with the observed values. However, the H band intensity from both points is still significantly below the observed values. For point 4, the J and H band intensities of the THEMIS model are comparable with the observed values, but the K_S band is over-estimated by a factor of two. Furthermore, the shape of the spectra of the model SIGMA agree with the observed spectra, although the intensity is clearly overestimated. Thus, the discrepancy between the observed and simulated values can be explained to a certain degree with uncertainties in the background sky estimate, but since the shape of the simulated spectra also depends on the optical depth of the cloud, uncertainties in the optical depth of themodel or variations in the relative abundances of the dust components will also affect the estimated intensity.

5.2 NIR extinction

In the models fitted to dust emission observations, the NIR extinction is typically three times higher than the direct NICER estimate of τ_J. In addition to NICER estimates using the Cardelli et al. (1989) extinction curve, we also calculated estimates for the extinction curves of the respective dust models. The values are computed as averages over a 5 × 5 pixel regionaround point 1 (the resolution of the optical depth maps was set to 40′′). For the Default model τ_Jext = 1.32, which is lower compared to the value derived using the Cardelli et al. (1989) extinction curve, τ_J,Card. = 1.50. For the models LG (model with grains up to 5 μm in size), LGM (as the LG model but the relative amount of large grains increased), and SIGMA, the optical depths are τ_J,ext = 2.19−3.0. For the model THEMIS the τ_J,ext = 1.42, which is within ~10% of the Cardelli et al. (1989) estimate.

The ratio for the Default model is τ_J,em∕τ_J,ext ~ 4.39, while it is lower, τ_J,em∕τ_J,ext = 1.55−1.66, for the models LG and LGM (Table 2). For the models with dust evolution, the ratios are smaller, 0.74 and 1.56 for THEMIS and SIGMA, respectively. Thus the models with larger grains, or some form of dust evolution, are more consistent between the sub-millimetre emission and NIR extinction.

5.3 Uncertainties of cloud models derived from FIR fits

We derived the cloud density distribution and strength of the external radiation field by fitting the observations of FIR dust emission. The estimated column density is sensitive to the dust temperature which in turn is sensitive to the strength of the radiation field. In this section, we quantify the dependences between the radiation field, the dust temperatures, and the optical depth of thecloud.

The effects of scaling and attenuation of the radiation field on the dust temperature can be solved from the equilibrium equation $${\displaystyle {\int }_0^{\infty }{Q}_{\text{abs}}}(\nu )\times I_{\text{ISRF}}\text{d}\nu ={\displaystyle {\int }_0^{\infty }{Q}_{\text{abs}}}(\nu )\times B_{\nu }(T)\text{d}\nu ,$$(6)

where Q(ν) are the dust absorption efficiencies, I_ISRF is the intensity of the radiation field, and B_ν(T) is a black-body function at temperature T. Figure 13 shows the dependence of the temperature of the Default and LGM dust grains on the energy density of the radiation field. Results are calculated for grains at an equilibrium temperature and using the Mathis et al. (1983) radiation field with a linear scaling factor ϵ and an attenuation by by e^−τ. We use values in the range of [0.1,10] for both ϵ and τ∕τ_J.

Compared to the Default model, the LGM model has more large grains and to reach similar temperature, a higher radiation field is required. Because of the higher emissivity, larger grains will have a lower temperature for a given radiation field and a lower column density is required to reach similar temperature as with the Default model.

Dust emission spectra are typically analysed as modified black-body radiation. However, the spectra will always deviate from this model, because of temperature variations in the sources and because the spectral index of dust opacity is not constant over the examined wavelength range. Figure 14 shows emission spectra obtained by multiplying the absorption cross sections of the Default and LGM models with a T = 15 K black-body function. The figure also shows the MBB fits using three points at 250, 350, and 500 μm and keeping both T and β as free parameters. The fitted temperatures T_c = 16.48 K and T_c = 15.85 K for the Default model and LGM model, respectively, are higher than the true temperature of 15 K. The fitted β values are β = 1.62 and β = 1.80, for the Default and LGM models respectively, lower than the β values of thedust models in the 250−350 and 350−500 μm intervals, β ~ 1.88 and β ~ 1.94 for the Default and LGM models, respectively. Thus, an observer relying on MBB fits would underestimate the column density (see also Fig. 10). Similar results of dust models containing only bare astrosilicates (Draine & Lee 1984; Draine & Li 2001, 2007) failing to reproduce the SED of the observed dust emission, have been discussed by Fanciullo et al. (2015), Planck Collaboration Int. XXII (2015), Planck Collaboration Int. XXIX (2016).

The errors in column density and radiation field will also affect the predictions of NIR scattering. However, the low intensity of our modelled NIR surface brightness is not caused by the uncertainties in the column density, but is related to the scattering efficiency of the dust models, or in the case of the evolved dust models (e.g. SIGMA and THEMIS) the chosen relative abundances of the different grain populations. Furthermore, MBB results are sensitive to temperature variations that can lead to a severe underestimation of dust column densities (Shetty et al. 2009; Juvela & Ysard 2012). For example, Pagani et al. (2015), who studied the cloud LDN 183, showed that Herschel observations can not be used to set strong constraints on the amount of very cold dust. Additional methods may be needed, such as observations of molecular lines and as well as NIR and MIR extinction.

The column densities based on the Default model vary within a factor of two, while the variation in the radiation field strengths is within a factor of ~1.5. The columndensity differences are not much larger between the SIGMA, THEMIS, and DDust models, but the differences in the radiation field strength reach a factor of two. However, the relationships between the true Q(ν)B(T_d) of the dust, the emission predicted by RT models, and the results of MBB fits are complex (Fig. 15). The emission from the dust grains might not follow any MBB curve, T_d and T_c differ and the fitted spectral index β differs from the opacity spectral index of the dust grains. These problems are acerbated when temperatures also vary within the beam.

We further quantify the discussion above by computing H band surface brightness maps using the Default model, with the assumption that we have underestimated the strength of the radiation field and overestimated the column density. Our aim is to see if by modifying these parameters, the Default model can be made to agree with the observations. We use a radiation field that is 1, 2, or 3 times and a column density that is 1, 0.6, or 0.3 times the value previously obtained for the Default model. Another component affecting the shape and strength of the NIR scattered light is the sky background behind the cloud, we assume that this is 1.0, 2.0, or 3.0 times the value derived from the DIRBE observations. The resulting H band surface brightness maps are shown in Fig. 16 and in Figs. E.6 and E.7. Increasing the strength of the radiation field will increase the surface brightness excess, but increasing the background sky brightness will decrease the surface brightness. The two values can be used to fix the level of the modelled surface brightness, but they are not enough to reproduce the observed surface brightness morphology (Fig. E.6). For that, also the cloud column density needs to be decreased.

Increasing the strength of the radiation field by a factor of two and simultaneously decreasing the column density by ~ 30%, (e.g. Fig. 16 panel G), the intensity of the simulated H band surface brightness is within 15% of the observed one. This can also be achieved by increasing the strength of the radiation field and the assumed background sky brightness by a factor of three and decreasing the column density by ~ 30% (e.g. Fig. 16 panel L). In both cases, the morphology of the surface brightness map in the central part is comparable with the observations, with a bright rim towards the south and decreasing surface brightness towards the core. However, the Veil (Fig. 9) is still prominent. Decreasing the density of the cloud further (Fig. E.7) the central region of the cloud becomes excessively compact, although the correct intensity can be reached by scaling the radiation field or the assumed intensity of the background sky brightness (e.g. Fig. 16 panels G, K, and L).

Although suitable parameters can be found to correct the H band intensity of the Default model, the required changes are substantial. Thus, it is evident that dust evolution needs to be taken into account, since the models SIGMA and THEMIS produce results closer to the observations without the need to fine tune other parameters.

Fig. 13 Comparison between the dust temperature and the energy density of the radiation field for the Default and LGM models. Thestrength of the radiation field has been normalised by the Mathis et al. (1983) model.

Fig. 14 Comparison of dust emission spectra and their MBB fits. The black and purple lines show the emission spectra for Default and LGM models with grains at 15 K temperature. The dashed lines show MBB spectra fitted to the 250, 350, and 500 μm points. The resulting fit parameters are Tc = 16.48, β = 1.62 and Tc = 15.85, β = 1.79, for the Default and LGM models, respectively. The red curve (scaled by a factor of 1.0 × 10−30 to fit the figure) is from a simulated Default model map pixel where the fitted temperature is similar to that of the black curve. The blue vertical lines show the 250, 350, and 500 μm frequencies.The relative differences between the black-bodies and the MBB fits near the 250, 350, and 500 μm points are shown in the plot below the emission spectra curves.

Fig. 15 Schematic overview between the assumptions that are made during a MBB fit and their comparison to the simulated spectra and a real cloud.

5.4 Effects of grain properties

Launhardt et al. (2013) showed that L1512 is a cold (T ≈ 12 K) core with high column density (N(H₂) ≈~2.0 × 10²²cm⁻²) and in the study by Lippok et al. (2013), the envelope of the core was better fitted by assuming a higher temperature compared to the core, which is consistent with the absence of internal heating. However, to explain the observed molecular abundances, Lippok et al. (2013) had to increase the hydrogen density derived by Launhardt et al. (2013) and the abundances of all modelled chemical species between a factor of two to three. In their chemical models Lippok et al. (2013) used a slightly coagulated grain model for the whole cloud with a coagulation time of 10⁵ yr at a gas density of 10⁵ cm⁻³ (Ossenkopf & Henning 1994). Lippok et al. (2013) discuss that by using a model of non-coagulated grains, they would have increased the hydrogen density by a factor of ~2.5.

Larger grains lead to stronger scattering as discussed by Steinacker et al. (2010) and Ysard et al. (2018), but the increase in grain size also increases the absorption coefficient. On the other hand, more complex or ‘fluffy’ grains can have a larger surface area, and have higher scattering efficiency with respect to their absorption efficiency (Lefèvre et al. 2016). This can be seen with the dense dust components of the models SIGMA and DDust (Fig. E.4), where the dense regions of the model have higher intensity and their morphology agrees better with the observations.

Our Default model consists of only PAH, silicate, and carbon grains, however, in dense regions of the ISM, gas phase freeze-out can create mantles on the grains. The initial studies trying to explain the detected coreshine required a high fraction of dust mass in large grains (Pagani et al. 2010; Steinacker et al. 2010; Andersen et al. 2013), with maximum grain size of the order of ~1 μm. However, to reach such large grain size through coagulation is difficult with respect to the time-scales, cloud densities, and turbulence (Steinacker et al. 2014). On the other hand, as discussed by Ysard et al. (2016), a dust model including mantle formation and low-level coagulation is able to reproduce the observed cloudshine levels without the need of very large grains. Our results using the THEMIS model show that the model can fit the observed emission, the intensity of the scattered light is within 15% of the observed values, and the morphology of the central region of the surface brightness maps is comparable to the observed surface brightness (Fig. 12. However, the intensity of the diffuse region surrounding the cloud is a factor of three to five higher than observed. Furthermore, as discussed by Bazell & Dwek (1990) and Ossenkopf (1991), porosity increases the absorption efficiencies at FIR and sub-millimetre wavelengths while the compositional inhomogeneities will have an effect upon the shape and strength of broadband features. Thus, these structural differences of the grains will affect both the true temperature and colour temperature derived from sub-millimetre flux ratios (Kruegel & Siebenmorgen 1994).

Thus it is evident that a larger maximum grain size or the relative amount of large grains in the dust mixture is not sufficient to reproduce the observed NIR surface brightness. Evolution of dust grains, in the form of aggregates and mantle formation, appears necessary. This is supported by our core temperature estimates that tend to overestimate the core temperature by ~ 1.5 K. However, the THEMIS model, with three different dust populations produces a colder core, with T_core = 7.6 K, similar to the estimated temperature of 8 K derived from line observations. However, even for the evolved grain models, the diffuse regions are still problematic as they are systematically brighter compared to the observations. The high surface brightness is likely related to the relative abundances of the dust populations, as for example, the THEMIS models has a relatively high abundance of the CMM population in the Veil and Striation regions, and similarly, the SIGMA model has a high abundance of the dense component in the Veil and Striations (Figs. A.1 and A.3). Decreasing the relative abundances of the dense dust components in these regions would decrease the surface brightness and improve the match with the observed surface-brightness morphology.

Fig. 16 Simulated H band surface brightness maps for the Default model with density scaled by a factor of 0.6 during the scattering computations and with different assumptions on the strength of the radiation field and on the sky brightness behind the cloud. Shown on the first row is the observed H band surface brightness map. Shown on the rows 2–4are the simulated H band maps, for which the strength of the radiation field has been scaled with a factor between 1 and 3, and the intensity of the background has been scaled between factors 1 and 3.0.

6 Conclusions

We have studied the cloud L1512 and modelled simultaneously the scattered NIR light at J, H, and K_S bands and the FIR emission at 250, 350, and 500 μm. We have used several dust models based on the Compiègne et al. (2011) dust and three separate dust models taking into account dust evolution. The radiation field used in the modelling is derived from the DIRBE observations. The NIR surface brightness is estimated using a density field and radiation field strength that are obtained from first fitting the FIR dust emission. The key result of our study are:

• The morphologies of the observed NIR surface brightness maps are in good agreement with the column density map derived from the Herschel observations. However, the low-column-density Veil seen above the cloud in the Herschel observations is not visible in scattered light.

• In the radiative transfer modelling, we can fit the observed emission with any of the tested dust models. The average fit residual in the 350 μm band are ± 1.5% and for the 250 and 500 μm bands at most − 15%. Depending on the model, the estimated H₂ column density at the cloud centre (point 1) ranges from 4.3 × 10²¹ cm⁻² to 1.6 × 10²² cm⁻², and the relative radiation field strength from 0.3 to 0.75.

• The core temperature of the radiative transfer models is on average ~1.5 K higher than the value ~8 ± 1 K, suggested by N₂H⁺ line observations. With T_core = 7.61 K, the THEMIS dust model gives the core temperature closest to the gas inferred value.

• The radiative transfer models matching the dust emission predict a J band optical depth that is on average three times higher than the value measured with the background stars. The exceptions are Scaled4 and THEMIS, with τ_J,em = 1.55 and τ_J,em = 1.05, respectively, in agreement with the value τ_J,Card = 1.5 derived using the Cardelli et al. (1989) extinction curve.

• The NICER estimates of τ_J,ext obtained with the NIR extinction curves of the tested dust models are mostly within 15% of the values obtained with the Cardelli et al. (1989) extinction curve. However, dust models containing larger grains (e.g. LG and LGM) increase τ_J,ext by a factor of two.

• For the models based on the Compiègne et al. (2011) model, the predicted surface brightness excess $I_{\nu }^{\Delta }$ in the central region of the cloud is a factor of two to four below the observed values and the morphology of the simulated scattered light maps does not match the observations. Increasing the maximum grain size of the dust grains or extending the width of the cloud along the line-of-sight will increase the intensity of the scattered light, but not enough. Thus, dust grain evolution (e.g. aggregates) is needed.

• Increasing the FIR emissivity by a factor of two (model Scaled2), increases the predicted NIR surface brightness by up to a factor of two. It also decreases the column density and produces more compact scattered NIR surface brightness maps. However, further increase of emissivity (model Scaled4) again decreases the NIR signal.

• The observed H-band surface brightness and its morphology could only be matched with the Default model by making substantial modifications to the values derived from the emission fitting and observations, indicating the need of changes in the dust properties.

• The observed thermal emission and scattered NIR surface brightness can be reasonably reproduced only by using dust models that take into account grain evolution. However, the J-band intensity and the H- and K_S-band intensities of the diffuse regions are far above the observed values.

It is easy to fit the dust emission alone but the simultaneous fitting of emission and scattering is challenging. The dust evolution must be taken into account, to produce sufficient amounts of scattered light and with the correct morphology. The uncertainty of the sky brightness behind the studied cloud can have a significant effect on both the intensity and morphology and should be constrained with a high precision before drawing conclusions on the dust properties.

Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Summary of dust models

A more detailed description of all dust models used in our radiative transfer modelling is provided in this section. Unless otherwise noted, the different dust models are based on the Compiègne et al. (2011) model. We use the size distributions and optical properties as included in the DustEM² (Compiègne et al. 2011).

Fig. A.1 Relative abundances of the diffuse and dense components for the DDust and SIGMA models across the innermost cells of the model cloud. Left panel: relative abundance of the diffuse component. Right panel: relative abundance of the dense component.

The changes that we have made to the Compiègne et al. (2011) model do not take into account any limits placed by the mineralogy or constraints placed by chemical abundances available in the ISM. The changes to for example the albedo or emissivity of the grains, are meant to be relatively small, but still large enough that differences between the models can be distinguished.

• Default: the Compiègne et al. (2011) model. Contains two populations of PAH grains with log-normal size distributions, a single component of small carbon grains with a log-normal distribution, and two components of large grains following power-law size distributions. The large grains consist of a population of carbon grains and a population of silicate grains. The average opacity spectral index in the range [250, 350] μm is β_250∕350 = 1.887.

• Albedo: the albedo of the grains at NIR wavelengths has been increased by 20% without changing the extinction cross sections during the scattering computations.

• Disttest: Compiègne et al. (2011) assumed a power-law size distribution for the large silicate and carbon grains. For this model the exponent factor of the power-law γ has been increased by 20% for both the silicate (original γ = −3.4) and carbon grains(original γ = −2.8) and β_250∕350 = 1.888.

• Gtest: the asymmetry parameter g of the grains has been increased by 15%.

• Scaled2: the emissivity of all wavelengths longer than 60 μm has been multiplied by a factor of two.

• Scaled4: as Scaled2, but the emissivity has been multiplied by a factor of four.

• Wide: the FWHM of the Gaussian distribution describing the line-of-sight density distribution has been increased from σ = 1.0 to σ = 1.73.

• LG: extended the maximum grain size of the Compiègne et al. (2011) model to 5 μm. The value of β_250∕350 is 1.903.

• LGM:as LG, but in addition increased the relative abundance of the large grains by a factor of two and decreased the relative abundance of the small carbon grains and PAH components by a factor of two. The value of β_250∕350 is 1.948.

• Ddust: combination of Default and LG dust models, the relative abundance of the latter increasing with density as $$S_{\text{LG}}=\frac{1}{2}+\frac{1}{2}\tan h\left(\frac{2.0\times \log (\rho )}{{\rho }_0}\right),$$(A.1)

  where ρ is the density of the modelled cloud and ρ₀ defines the limit where the relative abundances of the two dust components are equal. An example of the abundances is shown in Fig. A.1. In our model we have set ρ₀ = 1.45. For the average of the Default and LG models β_250∕350 = 1.955.

• SIGMA: a model with two dust components, for the diffuse component we use our Default model, and the dense component is built using SIGMA (Simple Icy Grain Model for Aggregates, Lefèvre et al. 2019). The dense component consists of aggregate silicate and carbon grains from the Min et al. (2016) model (note that we do not use the iron sulphide component) with added ices (Pollack et al. 1994). The final mixture consists of 56.8% silicates, 18.9% carbons, and 24.3% ices. The size distribution has a maximum grain size of 10.0 μm (see Fig. A.2), however for this work we have truncated the maximum grain size at 5 μm, and use a porosity factor of 0.7 for the grains. The porosity was chosen accordingly to the size distribution obtained from coagulation computation. Such a high porosity was necessary to obtain favourable conditions to stick dust grains together (Ormel et al. 2011; Hirashita & Nozawa 2013; Hirashita & Voshchinnikov 2014, Pagani et al., in prep.). In practice, there is a significant fraction of fluffy icy aggregates with a fractal degree close to 2. The dust distribution is representative of a coagulation at a constant density of 1 × 10⁶ cm³ with a constant porosity of 0.7 for 0.436 My. The simplification of dust distribution evolving at constant density and porosity will be discussed in a forthcoming paper (Pagani et al., in prep.) The relative abundances of the dense and diffuse dust components are set according to Eq. (A.1), with ρ₀ = 1.45. The threshold was so that the diffuse dust component would have a relative abundance of less than ~ 10% in the core of the model. For the average of the two dust components β_250∕350 = 2.024.

• THEMIS: a dust model as discussed by Köhler et al. (2015) and Ysard et al. (2016). The model consists of different dust population mixtures, that have a varying relative abundance based on the density of the model. In the following we adopt the naming convention as defined by Köhler et al. (2015). The diffuse regions of the model consist of mostly core-mantle (CM) grains which gradually evolve to core-mantle-mantle grains (CMM) as density increases. In the densest regions of the model, we assume that the grains have further evolved and are gradually replaced by aggregate CMM grains with additional ice mantles (AMMI). The relative abundances of the CM and CMM populations are set according to Eq. (A.1), with ρ₀ = 1.45. The relative abundance of the AMMI grains is then defined by taking all cells where ρ > 7.5, setting the relative abundance of these cells to 1.0 and smoothing the cells with a 3 × 3 Gaussian beam. We then reduce the relative abundance of the CMM population in these cells so that for each cell, the sum CM + CMM + AMMI = 1.0. An example of the relative abundances of the three components across the densest part of the cloud is shown in Fig. A.3. The size distributions of the different dust components are discussed by Köhler et al. (2015) and Ysard et al. (2016) (Figs. 1 and 2, respectively). The average β_250∕350 over the three dust components is 2.011.

Fig. A.2 Size distribution of the dust grains used for the dense component of the SIGMA model. For the emission fitting and scattering computations, the size distribution was truncated at 5 μm, as indicated by the purple line.

Fig. A.3 Relative abundances of the CM, CMM, and AMMI dust populations of the THEMIS model across the innermost cells of the model cloud. Left panel: relative abundance of the diffuse CM component. Centre: relative abundance of the CMM component. Right panel: relative abundance of the dense AMMI component.

Appendix B Stochastically heated grains

The effects of stochastic heating of small grains was not included in our modelling, because they are not directly connected to the submm emission and NIR scattering. However, in order to study how well our best-fit models agree with the observed emission in the MIR wavelengths, we computed two test cases using the models Default and THEMIS, including the stochastic heating. We did not fit the model to the observed emission, but rather used the best-fit parameters from the computations without SHGs and re-computed the emission including the SHGs.

The results of these computations are shown in Figs. B.1 and B.2. For both models, the emission in the SPIRE bands has decreased, by ~10−15%, because some of the energy is now emitted in the MIR wavelengths. The emission is at 100 μm 30–45% and at 160 μm 15% above the observed values. For the Default model, the morphology of the simulated 100 and 160 μm maps agreeswith the observed maps. For the THEMIS model, the bright rim in both maps is up to 40 MJy sr⁻¹ brighter than observed. The region where the CMM grains transition to the AMMI grains is clearly visible, especially in the 250 μm map.

An SED from the observations, computed as averages over a region corresponding to the red circle in Fig. 5 is shown in Fig. B.3 (the orange line). The observations between 3.6 and 8 μm are from the Spitzer IRAC instrument and the 24 μm data are from the Spitzer MIPS instrument. The MIPS data are corrected with a aperture correction of 2 MJy sr⁻¹. For the 60 μm observations, we use the improved reprocessing of the IRAS survey (IRIS) data (Miville-Deschênes & Lagache 2005) and the 100–500 μm observations are from the Herschel PACS (100 and 160 μm) and SPIRE (250, 350, and 500 μm) instruments.

An SED from the Default and THEMIS models with the SHGs included are shown in Fig. B.3. The background levels of the IRAC observations have not been calibrated to an absolute level. Thus, we assume an uncertainty of ± 30% for the simulated surface brightness at IRAC wavelengths and because of the uncertainty of the colour corrections, we assume an uncertainty of ± 30% for the 24 μm data. As discussed by Launhardt et al. (2013), the nominal calibration uncertainties for PACS and SPIRE are ~ 5 and ~ 7%, respectively,thus for the simulated MIR and FIR maps we assume a flat uncertainty of 10%. These uncertainties are highlighted with red and blue shading in Fig. B.3.

It is clear that the simulations do not precisely match the observations, although the shapes of the SEDs are in agreement. The THEMIS model is closer to the observations both in the NIR and FIR, although the Default model at 60 and 160 μm is closer to theobserved values.

Fig. B.1 Observed (first row) and simulated (second row) emission maps using the Default model at 100, 160, and 250 μm. Third row: difference between the observed and simulated maps. The simulated emission has been computed from the best fit parameters of the Default model and with taking into account stochastically heated grains.

Fig. B.2 Same as Fig. B.1, but for the THEMIS model.

Fig. B.3 Spectral energy distributions for the observations (orange line), the Default model (star symbols), and the THEMIS model (diamond symbols). The red and blue highlights show the assumed uncertainty in the modelled results for the Default and THEMIS models, respectively. The data points of the models are a sum of scattered surface brightness seen over the background level, estimated emission, and subtracting the background seen trough the cloud.

Appendix C Analysisof dust models based on the default model

In addition to the models discussed in the main part of the text, we have tested dust models that are modifications of the Default model, increased albedo (model Albedo), increased the value of the g parameter (model Gtest), and changes to the size distribution of the grains (Disttest). We have also tested a case where the LOS density of the cloud is wider compared to our assumed ’compact core’ model (model Wide). A description of these models is provided in Appendix A, and the results are summarised in Table C.1. The resulting dust emission and scattered surface brightness maps are shown in Appendix E. The scaling of the radiation field is similar to all models with K_ISRF ~ 0.45, expect for the model Wide which has a lower scaling factor of K_ISRF ~ 0.37.

Compared to the Default model the scattered surface brightness maps (see Figs. E.1–E.3) of models Albedo, Disttest, Gtest and Wide, show only minor differences. The most notable difference is for the model Albedo, which produces ~ 20% more surface brightness in all three bands compared to the Default model.

Increasing the line-of-sight extent of the cloud can produce more scattered light as the amount of illumination along the line-of-sight increases. We test a case, with the Default model, where the line-of-sight density distribution is extended,case Wide. The results indicate that the surface brightness is increased but only by ~ 10% (see Fig. E.3 second row).

Appendix D NIR and MIR observations

The MIR Spitzer observations were acquired from the Spitzer Heritage Archive. The cloud L1512 has been covered by three programmes, Programme id. 94 (PI: Charles Lawrence) and Programme id. 90109 (PI: Roberta Paladini), both using the InfraRed Array Camera (IRAC; see Fig. D.2) and in Programme id. 53 (PI: George Rieke) with the Multiband Imaging Photometer for Spitzer (MIPS). However, the programme 90109 was carried out during the warm mission, thus only the 3.6 and 4.5 μm channels were available. The IRAC observations are discussed in more detail by Stutz et al. (2009) and Steinacker et al. (2015), for the cold and warm missions, respectively. The MIPS observations are described by Rieke & Keene (2004) and Stutz et al. (2007).

The Spitzer 3.6 and 4.5 μm maps show extended surface brightness towards the central region of the cloud, but in both 5.8 and 8.0 μm maps the region is seen in absorption (see Fig. D.2). The surface brightness in the 3.6 and 4.5 μm maps can be caused by thermal emission by small grains, but the surface brightness is only seen from the dense central regions. If the surface brightness iscaused by thermal emission because of the high optical depth for the radiation heating the dust grains, one would expect it to be bright in the more extended region, not in the cloud centre. On the other hand, as discussed by Steinacker et al. (2010), in the current PAH emission models, the emission should increase towards the longer wavelengths and the Spitzer 4.5 and 5.8 μm bands, but in the Spitzer images the opposite is seen as the cloud in seen in absorption towards the longer wavelengths. Thus, we can conclude that the extended surface brightness seen in the 3.6 μm band, and at shorter wavelengths, is scattered light.

Shown in Figs. D.1 and D.2 are the observations from the WIRCam instrument and Spitzer space telescope. To study the diffuse signal, we have used the Source-Extractor (Bertin & Arnouts 1996, SExtractor) and Point Spread Function Extractor (Bertin 2011, PSFEx) to remove point sources. The extraction was carried out in three steps, in the first step, we use SExtractor to detect only the brightest point sources in each image. The detected bright sources are then analysed by PSFEx to construct a PSF for each detected source. The PSFs are then used in a second run with SExtractor to detect all point sources in the images. This second run produces an image of objects which we subtract from the original observationsresulting in an image in which stars appear as smooth holes. The NIR images are further smoothed over 6 × 6 map pixels to better show the morphology of the scattered light.

Fig. D.1 Colour maps show the NIR surface brightness at J band (left frame), H band (centre), and KS band (right). The surface brightness maps have been smoothed over 6 × 6 map pixels. The white contours show the N(H2) column density derived from the Herschel observations. The contour levels are 15, 45, and 75% of the peak column density of 1.12 × 1022 cm−2.

Fig. D.2 Spitzer observations covering the 3.6, 4.5, 5.8, and 8.0 μm bands. The surface brightness maps have been smoothed over 6 × 6 map pixels. The white contours show the N(H2) column density derived from the Herschel observations. The contour levels are 15, 45, and 75% of the peak column density of 1.12 × 1022 cm−2.

Appendix E Additional figures

Figures. E.1–E.4 show all of our surface brightness maps for the different dust models. Shown in Fig. E.5 are the individual components of the simulations.

Figures E.6–E.7 are the simulated surface brightness maps of the H band using the Default dust model. The rows and columns of the figures corresponding to different assumptions on the strength of the radiation field and the sky brightness behind the cloud. The model cloud column density is scaled by 1, 0.6, and 0.3 for the three figures, respectively.

Shown in Figs. E.8–E.10 are the residuals between our models and the observed emission (observed value minus the model prediction) for all dust models. Each row in the figures corresponds to a single model and the column show the reiduals in the 250, 350, and 500 μm bands. Additional NIR spectra extracted from point 2, see Fig. 11, is shown in Fig. E.11 and a schematic overview of our modelling process is shown in Fig. E.12.

Fig. E.1 Observed NIR surface brightness (first row) compared to model predictions with the Default (second row), Albedo (third row), and Disttest (fourth row) models.

Fig. E.2 Same as Fig. E.1, but for models Gtest (second row), Scaled2 (third row), and scaled4 (fourth row).

Fig. E.3 Same as Fig. E.1, but for models Wide (second row), LG (third row), and LGM (fourth row).

Fig. E.4 Same as Fig. E.1, but for models SIGMA (second row), THEMIS (third row), and DDust (fourth row).

Fig. E.5 Different components in the surface brightness maps for the Default model. Shown on the first row are the background subtracted observed surface brightnesses and the second row shows the simulated scattered light without the background subtraction. Shown on the third row is the component of the background that is seen trough the cloud (attenuated by a factor of e−τ) and the fourth row shows the background subtracted simulated surface brightness. The last row shows the difference between the observed and simulated background subtracted surface brightness.

Fig. E.6 Simulated H band surface brightness maps for the Default model with different assumptions on the strength of the radiation field and on the sky brightness behind the cloud. Shown on the first row are the observed surface brightness maps. Shown on the rows 2–4 are the simulated H band maps, for which the strength of the radiation field has been scaled with a factor between 1 and 3, and the intensity of the background has been scaled between factors 1 and 0.3. The density of the cloud is the same as in the Default model.

Fig. E.7 Same as Fig. E.6, but the cloud model has a 70% lower column density compared to the Default model.

Fig. E.8 Relative errors between the observed emission and our model predictions for 250, 350, and 500 μm bands. Each row corresponds to a different dust model. The rows from top to bottom correspond to models Default, Albedo, Disttest, and Gtest.

Fig. E.9 Same as Fig. E.8, but for models Scaled2, Scaled4, Wide, and LG.

Fig. E.10 Same as Fig. E.8, but for models LGM, SIGMA, THEMIS, and DDust.

Fig. E.11 J, H, and KS band intensities from observations (horizontal lines) and different models (symbols) for the map position 2. The colours correspond to the J (purple), H (red), and KS (black) bands. All intensity values have been background subtracted. We assume a 20% uncertainty in background sky estimates for the J and KS bands and an uncertainty of 30% for the H band.

Fig. E.12 Schematic overview of our modelling process to fit the dust emission and to compute an estimate for the scattered surface brightness.

References

All Tables

All Figures

	Fig. 1 Three colour image of the field L1512. The colours correspond to the NIR surface brightness at J (blue colour), H (green), and KS bands (red).All point sources have been masked.
In the text

	Fig. 2 Intensity of the radiation field, green line, as a function of wavelength. The purple line shows the Mathis et al. (1983) model. The black dots show the intensity of the DIRBE observations averaged over the HEALPIX map.
In the text

	Fig. 3 Model for the all-sky sky brightness distribution in the J band. The map has been rotated so that the direction towards L1512 is in the centre. The intensity of the model has been truncated at 20 MJy sr−1 for plotting and the image has been smoothed by a Gaussian with FWHM = 2°.
In the text

	Fig. 4 Comparison between the optical depth per H as a function of frequency for the Default, LGM, SIGMA, and THEMIS models. The yellow shaded regions show the frequency range from 250 to 500 μm and from 1.2 to 2.5 μm.
In the text

Fig. 5 Contours of J-band surface brightness on Herschel column density map. The plus signs show the locations of the 250, 350, and 500 μm emission maxima (red, green, and black, respectively). The white arrows indicate the projected direction of the three brightest nearby stars and the numbers next to the arrows indicate the 3D distance to the star and the Gaia G band average magnitude. The red circle shows the area that is used to compute the KISRF scaling factor.

In the text

	Fig. 6 J band extinction map and the optical depth at 250 μm. The white contours in the J band extinction map show the estimated uncertainties of 0.21, 0.3, 0.41 mag, from lowest to highest contour level.
In the text

	Fig. 7 Correlation between the optical depth at 250 μm and the J band extinction. The purple line is a least squares fit to the data and the Pearson correlation coefficient is 0.65 (p ≪ 0.01).
In the text

	Fig. 8 Observed (first row) and simulated (second row) emission maps using the Default model at 250, 350, and 500 μm. Third row: difference between the observations and simulations.
In the text

	Fig. 9 Observed (first row) and estimated surface brightness maps from the Default model at the J, H, and KS bands (second row). All maps have been background substracted.
In the text

	Fig. 10 Column density of H2 derived from Herschel observations and the comparison between column density profiles for selected models. The profiles have been computed as average values over 3 × 3 pixels along the arrow in the left panel.
In the text

Fig. 11 Observed map of the J band scattered light (left panel), the red numbers from one to four indicate the locations from where we have extracted NIR spectra, shown in the right panel. The solid lines show the observed spectra and the dashed lines show the spectra derived from the Default model. The dotted lines show the simulated spectra derived from the Default model, but the intensity of the scattered light has been multiplied by a factor of three before background subtraction.

In the text

Fig. 12 J, H, and KS band intensities from observations (horizontal lines) and different models (symbols) for the map point 1 (lower panel), point 3 (middle panel), and point 4 (upper panel). The colours correspond to the J (purple), H (red), and KS (black) bands. All intensity values have been background subtracted. We assume a 20% uncertainty in background sky estimates for the J and KS bands and an uncertainty of 30% for the H band.

In the text

	Fig. 13 Comparison between the dust temperature and the energy density of the radiation field for the Default and LGM models. Thestrength of the radiation field has been normalised by the Mathis et al. (1983) model.
In the text

Fig. 14 Comparison of dust emission spectra and their MBB fits. The black and purple lines show the emission spectra for Default and LGM models with grains at 15 K temperature. The dashed lines show MBB spectra fitted to the 250, 350, and 500 μm points. The resulting fit parameters are Tc = 16.48, β = 1.62 and Tc = 15.85, β = 1.79, for the Default and LGM models, respectively. The red curve (scaled by a factor of 1.0 × 10−30 to fit the figure) is from a simulated Default model map pixel where the fitted temperature is similar to that of the black curve. The blue vertical lines show the 250, 350, and 500 μm frequencies.The relative differences between the black-bodies and the MBB fits near the 250, 350, and 500 μm points are shown in the plot below the emission spectra curves.

In the text

	Fig. 15 Schematic overview between the assumptions that are made during a MBB fit and their comparison to the simulated spectra and a real cloud.
In the text

Fig. 16 Simulated H band surface brightness maps for the Default model with density scaled by a factor of 0.6 during the scattering computations and with different assumptions on the strength of the radiation field and on the sky brightness behind the cloud. Shown on the first row is the observed H band surface brightness map. Shown on the rows 2–4are the simulated H band maps, for which the strength of the radiation field has been scaled with a factor between 1 and 3, and the intensity of the background has been scaled between factors 1 and 3.0.

In the text

	Fig. A.1 Relative abundances of the diffuse and dense components for the DDust and SIGMA models across the innermost cells of the model cloud. Left panel: relative abundance of the diffuse component. Right panel: relative abundance of the dense component.
In the text

	Fig. A.2 Size distribution of the dust grains used for the dense component of the SIGMA model. For the emission fitting and scattering computations, the size distribution was truncated at 5 μm, as indicated by the purple line.
In the text

	Fig. A.3 Relative abundances of the CM, CMM, and AMMI dust populations of the THEMIS model across the innermost cells of the model cloud. Left panel: relative abundance of the diffuse CM component. Centre: relative abundance of the CMM component. Right panel: relative abundance of the dense AMMI component.
In the text

	Fig. B.1 Observed (first row) and simulated (second row) emission maps using the Default model at 100, 160, and 250 μm. Third row: difference between the observed and simulated maps. The simulated emission has been computed from the best fit parameters of the Default model and with taking into account stochastically heated grains.
In the text

	Fig. B.2 Same as Fig. B.1, but for the THEMIS model.
In the text

Fig. B.3 Spectral energy distributions for the observations (orange line), the Default model (star symbols), and the THEMIS model (diamond symbols). The red and blue highlights show the assumed uncertainty in the modelled results for the Default and THEMIS models, respectively. The data points of the models are a sum of scattered surface brightness seen over the background level, estimated emission, and subtracting the background seen trough the cloud.

In the text

	Fig. D.1 Colour maps show the NIR surface brightness at J band (left frame), H band (centre), and KS band (right). The surface brightness maps have been smoothed over 6 × 6 map pixels. The white contours show the N(H2) column density derived from the Herschel observations. The contour levels are 15, 45, and 75% of the peak column density of 1.12 × 1022 cm−2.
In the text

	Fig. D.2 Spitzer observations covering the 3.6, 4.5, 5.8, and 8.0 μm bands. The surface brightness maps have been smoothed over 6 × 6 map pixels. The white contours show the N(H2) column density derived from the Herschel observations. The contour levels are 15, 45, and 75% of the peak column density of 1.12 × 1022 cm−2.
In the text

	Fig. E.1 Observed NIR surface brightness (first row) compared to model predictions with the Default (second row), Albedo (third row), and Disttest (fourth row) models.
In the text

	Fig. E.2 Same as Fig. E.1, but for models Gtest (second row), Scaled2 (third row), and scaled4 (fourth row).
In the text

	Fig. E.3 Same as Fig. E.1, but for models Wide (second row), LG (third row), and LGM (fourth row).
In the text

	Fig. E.4 Same as Fig. E.1, but for models SIGMA (second row), THEMIS (third row), and DDust (fourth row).
In the text

Fig. E.5 Different components in the surface brightness maps for the Default model. Shown on the first row are the background subtracted observed surface brightnesses and the second row shows the simulated scattered light without the background subtraction. Shown on the third row is the component of the background that is seen trough the cloud (attenuated by a factor of e−τ) and the fourth row shows the background subtracted simulated surface brightness. The last row shows the difference between the observed and simulated background subtracted surface brightness.

In the text

Fig. E.6 Simulated H band surface brightness maps for the Default model with different assumptions on the strength of the radiation field and on the sky brightness behind the cloud. Shown on the first row are the observed surface brightness maps. Shown on the rows 2–4 are the simulated H band maps, for which the strength of the radiation field has been scaled with a factor between 1 and 3, and the intensity of the background has been scaled between factors 1 and 0.3. The density of the cloud is the same as in the Default model.

In the text

	Fig. E.7 Same as Fig. E.6, but the cloud model has a 70% lower column density compared to the Default model.
In the text

	Fig. E.8 Relative errors between the observed emission and our model predictions for 250, 350, and 500 μm bands. Each row corresponds to a different dust model. The rows from top to bottom correspond to models Default, Albedo, Disttest, and Gtest.
In the text

	Fig. E.9 Same as Fig. E.8, but for models Scaled2, Scaled4, Wide, and LG.
In the text

	Fig. E.10 Same as Fig. E.8, but for models LGM, SIGMA, THEMIS, and DDust.
In the text

	Fig. E.11 J, H, and KS band intensities from observations (horizontal lines) and different models (symbols) for the map position 2. The colours correspond to the J (purple), H (red), and KS (black) bands. All intensity values have been background subtracted. We assume a 20% uncertainty in background sky estimates for the J and KS bands and an uncertainty of 30% for the H band.
In the text

	Fig. E.12 Schematic overview of our modelling process to fit the dust emission and to compute an estimate for the scattered surface brightness.
In the text