EDP Sciences
Free Access
Issue
A&A
Volume 592, August 2016
Article Number A61
Number of page(s) 16
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201525792
Published online 25 July 2016

© ESO 2016

1. Introduction

Stars ultimately form through the gravitational collapse of cold and dense molecular cloud cores, irrespective of how these cores were formed in the first place. The initial temperature and density structure of such gravitationally bound cloud cores are important properties that determine the onset and course of the collapse. Yet, until very recently we had little direct observational evidence of the internal temperature structure of such cores. Consequently, the derivation of density profiles from dust emission data was often done with the simplifying assumption of isothermality (e.g., Ward-Thompson et al. 1999; Kirk et al. 2005; Launhardt 2005; Launhardt et al. 2010) or relied on weakly constrained model temperature profiles (e.g., Evans et al. 2001; Zucconi et al. 2001; Shirley et al. 2005). Only more recently have some studies tried to constrain the temperature distribution of starless molecular cloud cores directly from observational data of the thermal dust continuum emission (e.g., Ward-Thompson et al. 2002; Schnee & Goodman 2005; Stamatellos et al. 2007; Stutz et al. 2010; Nielbock et al. 2012; Beuther et al. 2012; Lippok et al. 2013; Launhardt et al. 2013; Pitann et al. 2013; Schmalzl et al. 2014).

The thermal emission from dust grains is indeed the most robust tracer of the temperature and density structure of such cold and dense cloud cores (see discussion in Launhardt et al. 2013). At the typical temperature of starless cores, 620 K, the dust grains emit thermal radiation mainly at far-infrared (FIR) wavelengths. The spectral shape of this emission depends on the temperature (and density) distribution along the line-of-sight (LOS), the optical depth of the emitting layer, and the properties of the dust grains. Longward of λ ≈ 200 μm, the thermal emission is usually optically thin even at column densities of 1025 cm-2 and thus traces well the structure in the interior of the dense cores, provided the dust temperature and the optical properties of the grains are known.

Thus, to derive the density structure from the thermal dust emission, cores must be observed in the FIR, in particular, toward the peak of their thermal spectral energy distribution (SED). Since the FIR is difficult to observe from the ground, most previous studies generally lacked these key observations or had to rely on relatively low-resolution FIR flux measurements from space observatories like IRAS (12 at 60100 μm), ISO (13 at 60200 μm), Spitzer (18′′40′′ at 70160 μm), or AKARI (60′′90′′ at 65160 μm). The Herschel1 space observatory (Pilbratt et al. 2010) was the first facility to cover most of the FIR wavelength range with high sensitivity and at an angular resolution (6′′36′′ at 70500 μm; Aniano et al. 2011) that is comparable to the largest ground-based single-dish millimeter telescopes. Thus, with Herschel observations and complementary long-wavelengths ground-based data, we are able to derive the dust temperature structure of molecular cloud cores and to put more robust constraints on their density structure than what was possible in the pre-Herschel era.

With the goal of constraining the physical conditions in molecular cloud cores during the prestellar and the earliest protostellar stages of star formation, we initiated the Herschel guaranteed time key project “earliest phases of star formation” (EPoS; Ragan et al. 2012; Launhardt et al. 2013). As part of this project, we observed a sample of Bok globules at five continuum wavelengths between 100 μm and 500 μm. Bok globules are small, nearby, and relatively isolated molecular clouds with a single core or at most very few dense cores, which makes them the most simply-structured and easily observable star-forming units in the Galaxy and ideal laboratories to study in detail the initial conditions of isolated star formation. An important selection criterion for the globules targeted in the EPoS survey was that they are located in regions with exceptionally low FIR background confusion noise to allow for deep observations at 100 μm. This wavelength is crucial for deriving the temperature of the cold dust since it constrains the short-wavelength side of the SED peak. This also implies that our targets are all located outside the Galactic plane at latitudes of | b | ≈ 3.4...12.4°.

Using a robust temperature mapping algorithm that applies modified black-body fits to the dust emission spectra, we have already shown that these isolated cores are thermally dominated by heating from the interstellar radiation field (ISRF); they have warm (1420 K) exteriors and cooler interiors (<1114 K, Stutz et al. 2010; Launhardt et al. 2013). In addition, because of their isolation and simple structure, we were able to use a ray-tracing inversion technique to reconstruct the 3D temperature and density profiles of the globules (Nielbock et al. 2012; Lippok et al. 2013; Schmalzl et al. 2014). We found that the temperature profiles of the starless cores drop to 713 K toward their centers, while their outer detectable rims are typically warmer by 5.5 ± 2.5 K.

In this paper, we compare the dust temperature profiles derived using this ray-tracing inversion technique with self-consistent radiative transfer models and explore the physical conditions that lead to the observed temperature distributions. The paper is structured as follows. In Sect. 2, we describe the sources and data used in this study. In Sect. 3, we describe the overall strategy as well as the dust model and modeling tools. In Sect. 4, we compare the results for the dust temperature profiles derived by the ray-tracing inversion method with those predicted by the self-consistent radiative transfer models. The results and uncertainties are discussed in Sect. 5. Finally, Sect. 6 summarizes and concludes the paper.

2. Sources and data

Table 1

Source list.

Based on the results of earlier studies, we selected for the EPoS survey 12 nearby, isolated, and previously well-characterized Bok globules, all located in regions of exceptionally low cirrus confusion noise (see Launhardt et al. 2013, and references therein). Each of these globules contains one to two embedded cores which, in some cases, are of different nature or evolutionary status. For this study of the temperature structure of starless cores, we selected those six globules from the EPoS sample that contain at least one starless core. We note that three of the six selected globules also contain a protostar in addition to the starless core (see Table 1 and Fig. 5). However, for this study we mask out the protostars (see Sect. 4.1), analyze the starless cores only, and discuss the uncertainties introduced by the presence of the nearby protostars in Sect. 5.1. Throughout this paper, we use the same core name convention as in Launhardt et al. (2013), except where explicitly referring to the entire globule or where distinction between core and globule in single-core globules would be meaningless. Source names, coordinates, and distances of the selected cores are summarized in Table 1. Total gas masses and outer radii of the cores are in the ranges 2.614 M and 0.20.4 pc (see Table 2). We use maps of the dust continuum emission in the PACS 100 μm, 160 μm (Poglitsch et al. 2010), and SPIRE 250 μm, 350 μm, and 500 μm bands (Griffin et al. 2010) that were obtained as part of the Herschel guaranteed time key project EPoS and processed as described in Launhardt et al. (2013). These data are complemented by ground-based (sub-)mm observations at 450 μm, 850 μm, and 1.2 mm that are presented in Launhardt et al. (2010, 2013) and references therein.

3. Modeling methods

3.1. General strategy

The ray-tracing inversion technique was developed to infer the dust temperature and density structure of starless cores directly from the observed dust emission maps without the need to make assumptions about the specific physical conditions (Nielbock et al. 2012). Apart from scattering (which is neglected by the ray-tracing), both the ray-tracing inversion and self-consistent radiative transfer use the same basic equations to relate the dust properties (opacity, temperature, and optical depth) to the effectively emitted flux density. However, the two methods have different approaches: The ray-tracing inversion assumes pre-defined analytical parameterized prescriptions for the temperature and density profiles before iteratively optimizing their parameter values by independently fitting the flux densities (SEDs) at each map pixel. The profile equations and iteration method are described in Sect. 3.3 (see also Nielbock et al. 2012; Lippok et al. 2013; Schmalzl et al. 2014). The self-consistent radiative transfer, on the other hand, also uses pre-defined density profiles (e.g., from a physical model or here derived by the ray-tracing inversion technique), but involves a model of the ISRF to calculate self-consistently the corresponding equilibrium temperature distribution.

The ray-tracing technique has the advantage that it does not need to make assumptions about the physical conditions, but only assumes a dust opacity model and free parameterized profile shapes for the temperature and density. It can also capture deviations from spherical symmetry and reproduce more complicated source structures. With this approach, a large and densely sampled parameter space can be probed with relatively little computational effort. On the other hand, the ray-tracing technique does not directly provide constraints on the physical conditions (e.g., the ISRF) that lead to the observed temperature and density distributions. This short-coming can, however, also be advantageous when compared to self-consistent radiative transfer modeling. If the range of parameters probed by a given setup of forward modeling is too small, the resulting “best-fit” model would not reveal the actual physical conditions that characterize the observed target. In contrast, the ray-tracing inversion technique yields temperature and density profiles that are likely closer to the real distributions and may thus lead to better hints toward the actual physical conditions when afterwards verified by self-consistent radiative transfer models. Thus, the ray-tracing inversion technique is well-suited to exploring the dust emission data and derive the source structures, but should then be complemented by a radiative transfer model to better characterize the physical conditions and environments of the cores.

For exploring our Herschel and complementary ground-based continuum data of isolated molecular cloud cores, we therefore chose to first apply the ray-tracing inversion technique, and afterwards verify the results by self-consistent radiative transfer models. As we show in Sect. 3.3, the analytical parameterization of the radial density profiles used for the ray-tracing inversion indeed reflects the actual physical configurations of the cores and the solutions found by the ray-tracing inversion for both the density and the temperature is practically always unique. With the results in hand, we now determine the physical conditions that best match the inferred temperature profiles. For this purpose, we adopt the azimuthally averaged density profiles derived with the ray-tracing inversion as input for self-consistent 1D radiative transfer modeling and calculate the equilibrium temperature distributions for these density profiles. We use the same dust opacity model as in the ray-tracing inversion (Sect. 3.2), but vary both the total strength of the ISRF and the extinction by a surrounding envelope, as described in Sect. 3.4. We then determine how well and for what combination of ISRF strength and envelope extinction the predicted equilibrium temperature profiles agree with the ray-tracing results.

3.2. The dust model

For the purpose of this paper, we adopted a dust opacity model from Ossenkopf & Henning (1994). Their models assume dust grains consisting of a mixture of silicates and amorphous carbon with different levels of coagulation and ice layer coverage around the agglomerates. Starting with the Mathis-Rumpl-Nordsieck (MRN) size distribution (Mathis et al. 1977), the dust grains are processed (coagulated) within a certain time and at a certain density. We chose the model for moderately processed grains with a coagulation time of 105 yr at a density of 105 cm-3 and with thin ice mantles (hereafter called OH5a2). These properties are closest to those typically observed in the starless cores in our sample (see Table 2 and Lippok et al. 2013). The mass absorption coefficient of this dust model at λ = 1.2 mm is κ1.2 mm = 0.79 cm2 g-1 of dust. It is thus somewhere in the middle of the range of values covered by several other published dust models (≈0.2−2 cm2 g-1; e.g., Ossenkopf & Henning 1994; Weingartner & Draine 2001; Ormel et al. 2011). For converting the dust mass into hydrogen mass, we adopt a mean hydrogen-to-dust mass ratio in the solar neighborhood of MH/Md = 110 (e.g., Sodroski et al. 1997). To obtain the total gas mass, Mg, the hydrogen mass must be multiplied by the mean atomic weight per H nucleus of the ISM, for which we adopt a value of μ = 1.4 (Cox 2000), i.e., Mg/Md ≈ 150.

The radiative transfer calculations require dust opacity values over the wavelength range between 90 nm and 10 mm to sample the full spectrum of the ISRF (Sect. 3.4 and Fig. 4). The dust models of Ossenkopf & Henning (1994), however, cover only the wavelength range between 1 μm and 1.3 mm. Therefore, we extrapolate the dust opacities into the ultra-violet using the prescription of Cardelli et al. (1989,Eq. (1), Table 3, Col. 5), and from 1.3 to 10 mm using a simple power-law. Finally, since Ossenkopf & Henning (1994) do not give scattering efficiencies, we augmented the OH5a model with scattering efficiencies following the approach of Young & Evans (2005) and using the albedos from the WD3.1 model (Weingartner & Draine 2001). The corresponding opacity spectrum is shown in Fig. 1. The albedos from the WD5.5B model (Weingartner & Draine 2001) would probably be the more appropriate equivalent to the coagulated OH5a dust model, but given the impossibility of calculating scattering efficiencies for such complex dust models as Ossenkopf & Henning (1994), for this paper, we use the WD3.1 albedos and restrict ourselves to a discussion of the related uncertainties in Sect. 5.1.3.

thumbnail Fig. 1

Dust opacity models (extinction/absorption coefficient per g of dust) used in this paper. Black lines show the absorption (solid) and extinction (dashed) coefficients of the modified OH5a model (Ossenkopf & Henning 1994). The gray line shows the OH1 model (see Sect. 5.1.3; for clarity, only the absorption coefficient spectrum is shown). Vertical dashed lines mark the wavelength range of the original OH models, i.e., values outside this range are extrapolated (see text). Arrows indicate the wavelength bands of the dust emission data used in this paper.

Open with DEXTER

3.3. Ray-tracing inversion method

The ray-tracing inversion technique derives a best-fit (χ2 minimization) dust temperature and volume density distribution to the uncertainty-weighted dust continuum flux density maps of molecular cloud cores. To work properly, flux density maps at three or more wavelengths, optimally distributed around the peak of the thermal SED, are required. The method needs as input parameterized descriptions of the functional form of the density and temperature profiles, as well as a dust opacity model (Sect. 3.2). Since the profile parameter values are fit independently for each LoS, moderate deviations from spherical or elliptical symmetry can be easily reproduced. The technique was first described in Nielbock et al. (2012) and initial results for individual starless cores from the EPoS sample were presented in Nielbock et al. (2012), Lippok et al. (2013), and Schmalzl et al. (2014). Here we only list the analytical profile formulas used in this method to illustrate the meaning of the profile parameters listed in Table 2.

thumbnail Fig. 2

Modified normalized Plummer profiles (Eq. (1)) with parameters chosen to mimic: (1) a simple power law; (2) a flat-density core with smooth outer edge; (3) an isothermal BES; and (4) a profile like observed in B 68.

Open with DEXTER

For the density profile, we use a universal modified Plummer-like function, which can, for a given choice of parameters, mimic various kinds of profiles (Fig. 2), including a simple (truncated) power-law, a nearly constant density sphere (with smooth edge), or a Bonnor-Ebert (BE) density profile (Plummer 1911; Ebert 1955; Bonnor 1956; Whitworth & Ward-Thompson 2001; Nielbock et al. 2012). This modified Plummer profile is described by (1)where nH = 2 n(H2) + n(H) is the total number density of hydrogen nuclei. This profile (i) accounts for an inner flat-density core inside r0 with a peak density n0 = Δn + nout (usually nout ≪ Δn); (ii) approaches a power-law with exponent η at rr0; (iii) turns over at r2r0n/nout)1 /η into a flat-density halo approaching nout; and (iv) is cut off at rout. The tenuous envelopes, which we have shown in Launhardt et al. (2013) to surround most globules and which become often evident as cloudshine3 halo, are actually neither azimuthally symmetric nor always fully spatially recovered by our observations of the dust emission. Therefore, rout is only estimated from the circularized 500 μm emission profiles, which show a similar spatial extent as the cloudshine (see Figs. A.1 through A.12 in Launhardt et al. 2013). The values of the free profile parameters Δn, r0, η, and nout are iterated in the ray-tracing inversion until convergence is achieved between LoS (input) and PoS (plane of sky, output), as described in Nielbock et al. (2012).

For the temperature structure of the starless cores, we use an empirical profile function that resembles the radiation transfer equation for an externally heated cloud. It couples the local temperature to the effective optical depth toward the outer “rim” at which the ISRF impacts: (2)with ΔT = ToutTmin and the frequency-averaged effective optical depth (3)τ0 is an empirical (i.e. free) scaling parameter that accounts for the a priori unknown mean dust opacity and the SED of the UV radiation of the ISRF. The temperature at the core center, T0T(r = 0), converges to Tmin if τISRF ≫ 1. As for the density profile, the values of the free profile parameters, Tout, ΔT, and τ0, are iterated until convergence is reached between LoS input and PoS output. Figure 3 compares the equilibrium temperature profile for B 68, self-consistently modeled with radiative transfer, with an analytical profile fit using Eq. (2). The good agreement between the two profiles suggests that Eq. (2) provides a reasonable parameterization of the temperature profile of an externally heated starless core.

Table 2

Properties of the studied cores.

thumbnail Fig. 3

Comparison of the self-consistently calculated equilibrium temperature profile for B 68 (solid line) and an analytical profile fit following Eq. (2). The little kinks in the radiative transfer temperature curve at r ≈ 2000 and 7000 au are computational artefacts due to the high optical depth.

Open with DEXTER

We also tested the uniqueness of the ray-tracing solution in χ2 space. While the solution was found to be genuinely unique for most pixels and most density/temperature profiles, the confidence regions can become quite elongated (banana-shaped) and even break up into two or three minima for certain configurations and pixels, in particular in regions with steep density and temperature gradients (see also Juvela & Ysard 2012, for a discussion of this problem in a slightly different parameter space). However, the occasional jumping of the solution to one of these secondary (unphysical) minima can be easily circumvented by implementing a very low-weighted near-infrared (NIR) extinction map (or any other smooth column density proxy), which effectively suppresses the secondary minima in the χ2 plane without affecting the value of the primary solution. Thus, the solution is always unique in practice. Moreover, the relative smoothness of the resulting mid-plane dust temperature and integrated column density maps (Fig. 5), which are composed of independently fit map pixels with different wavelength coverage (due to different map sizes at different wavelengths) also shows that the solution of the ray-tracing inversion is stable against both noise in the individual flux maps as well as wavelength coverage (but see wavelength coverage requirement mentioned above).

3.4. Self-consistent radiative transfer modeling

We use the radiative transfer code HYPERION (Robitaille 2011) to self-consistently calculate the equilibrium temperature distributions of the starless cores. As inputs we use the azimuthally averaged radial density profiles derived with the ray-tracing inversion technique and the same OH5a dust opacity model. The only free parameter in the modeling is sISRF, the relative total strengths of the ISRF. In addition, we allow NH(rsym), the column density of the surrounding envelope that shields the core from the ISRF, to vary within its observational constraints (see Table 2). The best-fit values of sISRF and NH(rsym) are then determined by comparing the calculated temperature profiles with the ray-tracing results as described in Sect. 3.5.

In our physical model for the radiative transfer, the dust is heated by the ISRF alone and cools by emitting thermal radiation. Additional heating via, for example, collisions with molecules, which are themselves heated by cosmic rays, is neglected since previous studies showed that this effect is very weak even at the highest densities found in the starless cores considered here (< 1.5 K, Goldsmith 2001; Evans et al. 2001). Heating by a potential compression of the cores is also neglected since this becomes relevant only at lower temperatures and higher densities than present in our cores (Keto & Caselli 2010).

Since we only compare radial profiles and carry out 1D radiative transfer calculations, the model cores are spherically symmetric and the ISRF is considered isotropic. However, the (approximate) spherical symmetry of the actual cores breaks down at large radii in most cases (see Figs. A.1 through A.12 in Launhardt et al. 2013), such that the 1D approximation may no longer hold even with azimuthal averaging. We therefore determine (visually) and list in Table 2 for each core a radius rsym inside which the core can safely be considered spherically symmetric. Temperature profiles from the two methods are only compared inside rsym. To still account for the attenuation of the ISRF by the material outside of rsym, we estimate the mean attenuation of the ISRF at rsym from the observed LoS column density, NH(rsym), assuming spherical symmetry. The uncertainty range of NH(rsym), which is mostly related to the actual deviations from spherical symmetry of the envelope and which we derive from the azimuthal scatter of NH values at rsym, is accounted for by varying rout in Eq. (1). The observed values and uncertainty ranges of NH(rsym) are also listed in Table 2.

thumbnail Fig. 4

ISRF model developed by the GALPROP consortium for the solar neighbourhood (Ackermann et al. 2012) with the CMB already included (solid line). For comparison we also show the model by Black (1994) as parameterized in Zucconi et al. (2001).

Open with DEXTER

For the ISRF, we use a model that was developed by the GALPROP consortium (Ackermann et al. 2012). From their four-component model for the Milky Way (stellar, scattered, transient, thermal), we obtained a table with the SED at the grid position corresponding to the Sun’s location in the Milky Way (RGC = 8.5 kpc). Added to that model is the cosmic microwave background (CMB) as a 2.72 K black-body. Figure 4 shows the SED of our adopted ISRF model, along with that of the “classical” model by Black (1994). Compared to this older model, the GALPROP model has a 1.5 times higher flux at wavelengths shortward of 8 μm (mostly stellar contribution), which is the most important spectral range for dust heating. Both the total and relative contributions of the different components in the GALPROP ISRF model have uncertainties that are difficult to quantify. Furthermore, individual components of the ISRF may vary differently locally on spatial scales smaller than sampled by the model, for example, due to the proximity to luminous stars, to star-forming regions, or to molecular clouds that can also shield a region from the general ISRF (not only for locations inside the cloud, but also if a dark cloud is located between a globule and the galactic plane). Since we have no means of constraining such local variations of the individual components better than done in the GALPROP model, we took the simplifying approach to account for the uncertainty of the local ISRF by introducing the free scaling parameter sISRF, which multiplies the total strength of the GALPROP ISRF (i.e., excluding the CMB). In Sect. 5.1.4, we also test a different spectral shape of the ISRF and discuss its effect on the results.

3.5. Comparing ray-tracing and self-consistent radiative transfer results

With the setup described above, we calculated radiative transfer models for grids of the two free parameters sISRF and NH(rsym). The latter is only free within the observational constraints. The value of sISRF was varied in 19 steps between 0.1 and 5.0 for all sources. This total range is somewhat larger than the range of ISRF scaling factors listed by Ackermann et al. (2012) for different cosmic ray models and different locations in the Galaxy. Since all our sources are located in the Solar neighbourhood, are relatively isolated, and are not located close to star-forming regions or luminous individual stars, it should thus cover the uncertainty in our knowledge of the local strength of the ISRF.

The total number of models per source depended on the uncertainty range of NH(rsym) (Sect. 4.1 and Table 2) and ranged from 57 to 209. The best-fit values of sISRF and NH(rsym) (which we label and ) are then determined by means of χ2 minimization of the differences between the temperature profiles at rrsym predicted by the radiative transfer models and those derived with the ray-tracing inversion. In addition to deriving the best-fit parameter values of sISRF and NH(rsym), we evaluate the goodness of the fits by verifying if the best-fit equilibrium temperature distributions agree with the temperatures of the ray-tracing inversion technique within its uncertainty of K (see Sect. 5.1) everywhere within rsym.

4. Modeling results

4.1. Ray-tracing results

thumbnail Fig. 5

Mid-plane dust temperature maps (color, all maps on the same scale) of all six globules studied in this paper, derived with the ray-tracing inversion technique and overlaid with contours of the hydrogen column density (white: 1021 and 1022 cm-2, black: 0.5 and 5 × 1021 cm-2). Yellow squares indicate the center positions of the starless cores as listed in Table 1 and in Launhardt et al. (2013). Asterisks indicate the positions of embedded protostars (see Launhardt et al. 2013).

Open with DEXTER

thumbnail Fig. 6

Radial profiles of all six starless cores, derived with the ray-tracing inversion technique and plotted up to rsym. a) Mid-plane hydrogen volume number density vs. radius; b) hydrogen column density vs. impact radius; c) mid-plane dust temperature vs. radius.

Open with DEXTER

We derived the dust temperature and density structure of all six starless cores using the ray-tracing inversion technique and the OH5a dust opacity model to fit the dust continuum flux density maps in up to eight bands between 100 μm and 1.2 mm. The initially resulting maps of mid-plane dust temperature and density were already shown in Lippok et al. (2013). Like in Launhardt et al. (2013) and in Lippok et al. (2013), we derive a strong negative correlation between central dust temperature and peak column density, which suggests that the thermal structure of the cores is dominated by external heating from the ISRF and shielding by dusty envelopes.

For the present paper, we recalculated the temperature and density maps using an improved convergence scheme for the ray-tracing algorithm. This resulted in slightly different profile parameters from those of Lippok et al. (2013, Table 3), albeit qualitatively and quantitatively very similar dust temperature and density maps. The new values of the profile parameters are, however, within the formal uncertainties of the Lippok et al. (2013) results and the differences are thus not significant. The respective new mid-plane dust temperature maps and contours of the (integrated) column density are shown in Fig. 5. The azimuthally averaged radial profiles of the mid-plane volume density and dust temperature as well as of the integrated column density are shown in Fig. 6. The corresponding radial profile fit parameters according to Eqs. (1) and (2) are listed in Table 2. Note that η (Eq. (1)) is only listed for completeness here; its value has a large uncertainty since it is sensitive to how the relative weights of individual data points (pixels) in the profile fitting are scaled with radius and signal-to-noise ratio (see also uncertainty discussion in Sect. 5.1.3). Only these radial profiles are used as input for and comparison with the self-consistent 1D radiative transfer modeling.

thumbnail Fig. 7

Cumulative radial mass distribution of the six starless cores derived with the ray-tracing inversion technique and after masking out the nearby protostars in CB 17, CB26, and CB 244.

Open with DEXTER

Figure 7 shows the cumulative radial gas mass distributions of the six cores, derived from the column density maps after masking the overlapping protostars and their respective envelopes in CB 17, CB 26, and CB 244. The size of masking regions was chosen to be a compromise between avoiding the local surroundings of the warm protostars and using as much as possible the extended emission from the cold cores since the protostellar envelopes neither have sharp boundaries nor can they clearly be separated in the emission maps.

The physical outer radii of the globules are all within the narrow range (6.5 ± 2.5) × 104 au (Table 2), and, with the exception of CB 244, the cumulative mass distributions are basically flat beyond r ≈ 4.5 × 104 au (Fig. 7). This means that the more extended envelopes, which are partially visible as cloudshine (see Launhardt et al. 2013, Figs. A.1 through A.12), do not contribute much to the total mass. To ensure good comparability of the sources, we therefore chose to derive the total source masses from the cumulative mass distributions (Fig. 7) within a fixed radius of 5 × 104 au around the column density peaks of the starless cores. The resulting total gas masses are listed in Table 2 and are in the range 2.614 M.

thumbnail Fig. 8

Central density vs. total gas mass for the six starless cores (data from Table 2). The dashed line indicates the maximum stable density of a pressure-supported, self-gravitating modified (nonisothermal) BES as calculated by Keto & Caselli (2008, their model with photoelectric heating at the core boundary taken into account). Error bars represent the mean relative uncertainties of 30% on the cloud mass and 50% on the central density (see Sect. 5.1).

Open with DEXTER

To better characterize the cores and to help interpreting the results, we revisited the stability assessment of Launhardt et al. (2013) and Lippok et al. (2013). Figure 8 shows the central H2 volume density vs. the total gas mass for the six starless cores (see Table 2) along with the stability criterion by Keto & Caselli (2008) for pressure-confined, self-gravitating modified (non-isothermal) Bonnor-Ebert spheres (BES, Ebert 1955; Bonnor 1956). CB 17, CB 26, CB 27, and B 68 are located close to the boundary between stable and unstable, making their evolutionary path unpredictable. CB 4 is clearly a subcritical core and must be purely pressure-confined (or transient). CB 244 (SMM2) is a clearly supercritical core and thus a good place to look for infall signatures. These conclusions agree well with the analysis of Launhardt et al. (2013).

4.2. Radiative transfer results

We have already shown in Sect. 3.3 and Fig. 3 that using Eq. (2) for the temperature distribution in the ray-tracing inversion leads, in principle, to qualitatively and quantitatively very similar temperature profiles derived with this technique and predicted by the radiative transfer models. Before comparing both model results in more detail for the individual cores in Sect. 4.3 and discussing the uncertainties and implications in Sect. 5, we demonstrate here the systematic effects of varying the two free parameters used in the radiative transfer modeling, sISRF and NH(rsym).

thumbnail Fig. 9

Effect on the radial equilibrium temperature distribution predicted for B 68 from scaling the total flux of the GALPROP ISRF by factors sISRF = 0.5,1,2,3 (with OH5a dust and NH(rsym) = const. = 1.5 × 1021 H cm-3; black lines). Blue dots show the mid-plane dust temperature distribution inferred with the ray-tracing inversion.

Open with DEXTER

thumbnail Fig. 10

Effect on the radial equilibrium temperature distribution predicted for B 68 from varying the envelope extinction of B 68 via NH(rsym) for values of 0.4,0.8,1.7,2.5,4.2 × 1021 H cm-3 (with OH5a dust and sISRF = const. = 2.2; black lines). Blue dots show the mid-plane dust temperature distribution inferred with the ray-tracing inversion.

Open with DEXTER

Figure 9 compares the equilibrium temperature distributions in B 68 for different ISRF scaling factors sISRF (with OH5a dust and NH(rsym) = const. = 1.5 × 1021 H cm-3). We find that the total strength of the ISRF scales the temperature approximately equally everywhere throughout the entire core. Figure 10 compares temperature profiles of B 68 for different envelope extinctions (now leaving sISRF = const. = 2.2). In contrast to varying the total power of the ISRF, we find that changing the envelope extinction mainly affects the temperature at large radii only, while the effect on the core temperature is marginal. The reason for this behavior is that the outer layers mainly absorb the UV part of the ISRF, while the radiation at longer wavelengths, that heats the interiors, is nearly unaffected. Therefore, the temperature profiles for different envelope extinctions, but constant sISRF, converge towards the core centers. Hence, for a given dust model and spherical geometry, envelope extinction (or variations in the UV-to-IR flux ratio of the ISRF) controls the core-envelope temperature contrast, whereas scaling the entire ISRF mainly affects the mean temperature of the entire globule.

4.3. Comparison of ray-tracing and radiative transfer results

All six cores show a positive radial temperature gradient with cool interiors and warmer envelopes, as expected and predicted by radiative transfer models for clouds that are externally heated by the ISRF and shielded by dust. This general behavior is readily reproduced for all sources by the ray-tracing inversion of the dust emission maps (Fig. 6c). On a more quantitative side, both the central core temperatures of 7.512 K and the outer envelope temperatures of 13.519 K determined with the ray-tracing method (Table 2) are well within the ranges that are predicted by self-consistent radiative transfer models (Fig. 11).

Table 3

Best-fit values of sISRF and NH(rsym).

thumbnail Fig. 11

Comparison of the temperature distributions derived with the ray-tracing inversion technique (blue dots and shaded regions) to the best-matching thermal equilibrium temperature distributions (black lines). Blue dots represent pixel values of the mid-plane temperature map derived with the ray-tracing inversion method. Shaded regions indicate the uncertainty of the mean temperature profile which we estimate as K (Sect. 5.1). All models assume OH5a dust. Best-matching scaling factors sIRSF for the total flux of the ISRF are also indicated in each panel.

Open with DEXTER

Table 3 lists for all six cores the individual best-fit parameter values and , along with the minimal values from comparing the ray-tracing and radiative transfer results. From the χ2 maps, we also derive approximate mean relative 1σ uncertainties of % and % for all sources. Despite the qualitatively different effects on the equilibrium temperature distribution from varying sISRF and NH(rsym) (Figs. 9 and 10), the two parameters remain partially correlated and the uncertainties are therefore not completely independent. However, we find that solutions with sISRF deviating by more than a factor of two from our best-fit value do no longer fit the data at all. This should be kept in mind when interpreting the absolute uncertainty values for the individual sources listed in Table 3, which we derived from the mean relative uncertainties mentioned above.

The corresponding comparison between the temperature profiles from the two methods is shown in Fig. 11. The quality of the fits differs between the sources, which is reflected in both the values (Table 3) and in the plots (Fig. 11). The profiles for CB 4 and B 68 show excellent agreement. For CB 27 and CB 244, the respective profile shapes differ slightly, though the predicted equilibrium temperatures are still within the uncertainty range of the ray-tracing results of K everywhere within rsym (Sect. 5.1). For CB 17 and CB 26, this latter quality criterion could not be strictly met. In both cores we “measure” (with the ray-tracing inversion) a smaller temperature gradient than what is predicted by the radiative transfer models for even the largest allowable values of NH(rsym). In Sect. 5.1.2 we discuss possible reasons for these discrepancies.

With the exception of CB 4 (), the most distant source in our sample (Table 1) and the only core with clearly subcritical central density (Fig. 8), we find that the GALPROP ISRF generally results in temperatures that are significantly lower than what we derive from the data with the ray-tracing inversion. Indeed, the other five sources were best fit with . This trend was also found by, for example, Shirley et al. (2005) and Launhardt et al. (2013). However, we must caution that this finding only holds for the applied dust opacity model and dust models with other NIR–to–FIR opacity ratios may lead to other values of (see discussion of uncertainties in Sect. 5.1.3). Therefore, we refrain here from drawing general conclusions on the total strength of the local ISRF. The lower value for CB 4 may thus not necessarily suggest a weaker local ISRF, but may well be related to a different dust opacity law in this low-density cloud.

For most sources, we also had to increase the envelope extinction to the maximum value allowed by the observational constraints (Table 2). However, since the estimated uncertainty of the best-fit NH(rsym) values is of the same order as the uncertainty of the observational constraints, we do not consider this trend significant.

5. Discussion

5.1. Uncertainties

We estimate the total uncertainty of the dust temperature in the starless cores to be K at T ≈ 10 K along the entire radial profile of the sources. The most significant contribution to the uncertainty of the central temperature was found to come from the symmetry assumptions made in the ray-tracing inversion method. The most significant contribution to the uncertainty of the outer (envelope) temperature comes from the irregular geometry of the envelopes. The relative uncertainty of the central density is estimated to be σn,rel ≈ 35%, with the most significant contributions coming from the inversion method and the dust opacities. However, this latter statement assumes that we consider only reasonably applicable dust opacity models and not the entire range of the still observationally poorly constrained dust opacity models that are in the literature. The relative uncertainty of the values for the total core mass is estimated to be σM,rel ≈ 30%. The formal uncertainty on the derived relative strength of the ISRF is estimated to be σsISRF ≈ 20% when considering the modified OH5a dust opacity model only. However, as mentioned above (Sect. 4.3) and discussed further in Sect. 5.1.3, the actual uncertainty of the ISRF strength may be significantly larger when dust opacity models with different NIR–to–FIR opacity ratios are allowed. In the following, we discuss separately the contributions to these uncertainties coming from the data and the ray-tracing inversion method, the source properties, the dust opacity law, and the assumed ISRF.

5.1.1. Data and inversion method

The uncertainties introduced by the observing methods and data calibration and reduction procedures are analyzed and discussed in detail in Launhardt et al. (2013). Although the flux calibration of the Herschel data has been improved since then, the uncertainties are dominated by data reduction issues such as, for example, spatial filtering and background determination, such that the assessment in Launhardt et al. (2013) is still considered valid. Hence, we estimate the data-related uncertainty of the dust temperatures inferred by the ray-tracing inversion to be σT< ± 0.5 K at 10 K.

The uncertainties of the temperature mapping and ray-tracing inversion method have already been assessed in three preceding papers (Nielbock et al. 2012; Launhardt et al. 2013; Lippok et al. 2013). The uncertainties of the dust temperature estimation in the envelopes, where LoS temperature gradients are small, are dominated by the data-related uncertainties (see above). Towards the core centers, where flux densities are much higher, but also LoS temperature gradients, the local uncertainty of the derived core temperature is dominated by the fine-tuning of the profile parameters in the ray-tracing inversion, which is related to the symmetry assumptions and the imperfect convergence between PoS and LoS profile parameters (Sect. 3.3). Based on the analysis in Nielbock et al. (2012), we estimate the inversion method-related uncertainty of the central core temperatures to be K.

Since the uncertainty of the LoS mean (optical depth-weighted) temperature is much smaller than that of the local central temperature (see Launhardt et al. 2013), the method-related uncertainty of the column density is also small. From Monte Carlo tests, we estimate it to be σN0 ≲ ± 10%. Due to this balancing along the LoS, the method-related effect on the value of the central density is also smaller than one would infer from the uncertainty of the central temperature only. From Monte Carlo tests, we estimate it to be σn0 ≲ ± 20%. However, the latter two statements only hold for the most nearby and fully spatially resolved cores as we explain below.

We also tested the effects of beam smoothing by omitting the 500 μm data and working at the 350 μm resolution of 25′′, since the 500 μm Herschel beam of 36.̋4 (Aniano et al. 2011) might not fully resolve the coldest parts of the core centers (see also Nielbock et al. 2012; Schmalzl et al. 2014). For example, for the two most distant sources in our sample (CB 4 and CB 17), the flat-density core diameters (2 × r0, Table 2) are only 2.7 times larger than the beam. We did indeed find an increase of the central density and peak column density by up to 50% for these two sources when working with the smaller beam. However, the resulting values of the central temperature were always within this uncertainty range and we did not find any systematic trend towards lower temperatures. Since the temperature rises faster towards larger radii than the density decreases, that is, the radial temperature profiles are typically 3040% narrower than the density profiles, this effect does not arise from unresolving the cold, flat-density core. It can rather be attributed to smoothing the steep drop-off of the radial temperature profiles and thus effectively re-distributing the mass along the LoS in the radius range ≈5000−20 000 au (see Fig. 6). The amplitude of these effects, which are related to the aforementioned fine-tuning and convergence of PoS and LoS profile parameters, thus depends on the angular size and morphology of the individual cores and cannot be uniquely quantified. This analysis also illustrates one short-coming of our method and the advantage of forward-modeling which would take full advantage of the higher angular resolution of the shorter-wavelength data.

5.1.2. Source-specific properties

Even in the most simple-structured star-forming regions like Bok globules, the envelopes of the starless cores are not strictly spherical (e.g., Myers et al. 1991; Launhardt et al. 2013; Lippok et al. 2013). Furthermore, three out of six globules in our sample contain a nearby secondary core with an embedded protostar. In several cases, the starless cores are also not located in the center of the globules. We partially eluded these problems by masking those regions that clearly deviate from the symmetric structure. We also set an outer radius, rsym, for the comparison of ray-tracing and radiative transfer results where the assumption of symmetry breaks down or where the signal in the continuum maps gets too weak to allow for reliable derivations of the structure. While this procedure is to some extent arbitrary, we verified that the results on and are nearly unaffected when varying rsym within reasonable ranges. Nevertheless, this radially increasing deviation from spherical symmetry leads to an uncertainty also for the outer parts of the azimuthally averaged temperature profiles. From the azimuthal scatter of the best-fit mid-plane temperature values at rsym, we estimate this uncertainty to be also K.

The most important specific properties of the individual globules are described in Launhardt et al. (2010) and Launhardt et al. (2013) and summarized in Table 1. Here we only discuss those characteristics that are relevant for the interpretation of results from this paper. CB 4 and B 68 are the two cores where the temperature profiles from the ray-tracing technique and the radiative transfer models had excellent agreement. These two cores are also the most round ones. Additionally, these sources are both single cores4 without a nearby protostar, suggesting that they have the simplest structures in our sample.

Conversely, CB 27 and CB 244-SMM2, both of which showed moderate agreement between the two modeled temperature profiles, are more complex. CB 27 is the most elliptical globule in our sample (see Fig. 5), which may explain the deviations between the 1D radiative transfer models and the ray-tracing results. CB 244-SMM2 has a nearby bright protostellar core (SMM1, 90′′ or 18 000 au to the west), and our separation of the two cores may not have been perfect, such that the recovered dust emission toward the starless core may be enhanced by the presence of this protostar. Therefore, the ray-tracing inversion technique may not accurately portray the temperature profile of this source. In addition, but to a smaller extent, anisotropic heating of the starless core by this protostar may also play a role. However, we do not consider the local enhancement of the radiation field due to the protostar in SMM1 significant since it is located at a projected distance of 18 000 au, has an estimated bolometric luminosity of only 1.8 L (Launhardt et al. 2013), and its outflow cones (through which most of the NIR/MIR luminosity escapes) are oriented perpendicular to the connecting line towards SMM2 (see Yun & Clemens 1994).

Finally, for CB 17-SMM and CB 26-SMM2, the ray-tracing technique and the radiative transfer models produced inconsistent temperature profiles. For CB 17-SMM, a low-luminosity Class I young stellar object (YSO) is located ≈10′′ (2500 au) northwest from the center of the starless core (Schmalzl et al. 2014). This YSO impairs the flux profile reconstruction for the starless core and leads to larger systematic deviations than for the other globules as this blend probably leads to an overestimation of the inner temperature of the starless core in the ray-tracing inversion. This YSO may also locally enhance the radiation field to which the core SMM is exposed. However, as in the case of CB 244, we consider this effect less significant than the blending mentioned above since the YSO has a bolometric luminosity of only 0.12 L (Launhardt et al. 2013) and is not embedded in the core SMM (Chen et al. 2012). In addition, the core of CB 17-SMM is less resolved than most other globule cores, because of the combination of a relatively small physical size and a large distance (250 pc), such that beam-smearing may actually not be negligible for this core.

CB 26 is a globule with multiple cores and a relatively complex visual appearance (Launhardt et al. 2013, Fig. C.4). Moreover, C18O(2 – 1) spectra of the SMM2 core show multiple velocity components (Lippok et al. 2013, Fig. 3) which suggest that this core may actually be a pole-on view of a longer filament or the projection of two or more diffuse cores. If CB 26-SMM2 is indeed a by-chance projection of two or more diffuse cores in a filament seen pole-on, both the ray-tracing inversion and the 1D radiative transfer would fail to reproduce its physical structure and internal temperature distribution. For this reason, we do not further consider CB 26-SMM2 for the following discussion of the effects of the dust opacity law and the ISRF.

5.1.3. Dust opacities

The derivation of a temperature from the SED of the thermal dust emission is sensitive to the spectral slope β of the dust opacity law at FIR through mm wavelengths. However, most previous observations of cold and dense molecular cloud cores lack robust detections over this wide wavelength range, resulting in poorly constrained estimates of β. Thus, many studies instead adopt appropriate modeled dust opacity laws, most of which have slopes in the range β ≈ 1.9 ± 0.1 (e.g. Ossenkopf & Henning 1994; Weingartner & Draine 2001; Ormel et al. 2011). Within this range, the dust temperatures inferred from the SEDs would vary by less than ± 0.5 K (see Launhardt et al. 2013). However, the combination of these slight temperature differences with the much larger differences in the absolute values of the FIR – mm dust opacities between the various dust models results in an uncertainty of the density, column density, and mass of about a factor of three. This would in turn affect the shielding and energy balance of the cores and thus the predicted equilibrium temperature distributions and sISRF values (see also discussion below). Although the absolute scaling value for the FIR – mm opacity is unknown, one can generally exclude models with unprocessed ISM-type dust grains or models with extremely coagulated grains that also have very thick ice mantles. Such extreme cases are unlikely to represent the conditions in molecular cloud cores.

Ideally, we would use different dust opacity laws for the core centers and envelopes, as both regions differ in temperature and density and hence in the expected degree of grain processing. While the use of two different dust opacity laws would have little effect on the observed dust temperature profiles, there would be a more significant effect on the density profiles (and hence on the equilibrium temperature distribution and ISRF strength inferred from radiative transfer models). We tested such a two-layered dust model on B 68, using OH5a opacities at nH> 104 cm-3 (r> 1.5 × 104 au) and OH1 opacities outside. Compared to our single-opacity results with OH5a opacities alone (see Table 2), we found that the two-layered model gave consistent results for T0, a 0.8 K increase for Tout, a 25% increase of n0, a 50% increase of N0, and a lower value of η (4 instead of 5). This may partially explain our relatively large values of η in Table 2. Regarding the best-fit values of η, it should also be mentioned that there is a relatively strong degeneracy between r0 and η, such that the η values listed in Table 2 should not be over-interpreted. For the two-layered opacity model of B 68, we could, for example, still obtain reasonable fits to the data by forcing r0 to 20% below the formal best-fit value, which would decrease the value of η further from four to three. Nevertheless, we lack sufficient constraints to use such a multi-layered dust model, and we caution that our comparison for B 68 is not robust quantitatively. It only qualitatively reveals the trend of the systematic uncertainties introduced by using a single-layer dust opacity model for such cloud cores.

This well-known problem has been discussed recently by, for example, Pagani et al. (2015). Their analysis, like ours above, shows that the lack of good constraints on the dust opacity law and its possible change along the LoS can lead to significant uncertainties in the column densities and masses derived from dust emission data alone, in particular when very cold (Tin< 10 K) cores are involved. However, our more sophisticated treatment of the temperature structure of the starless cores as compared to Pagani et al. (2015) leads to smaller uncertainties than derived by these authors. Furthermore, we can show that the larger uncertainty on the values of the dust opacity and temperature in the central core region with high density has very little effect of the estimate of the total core mass. The radial profile and mass distribution diagrams for B 68, CB 27, and CB 244-SMM2 (Fig. 12) clearly demonstrate that the cold and dense innermost regions of this core actually contribute very little to the total mass. Most of the core mass comes from a shell with ≈5000−15 000 au radius, where the dust temperature is ≈2−5 K higher than at the core centers and the densities are a few 104 cm-3. If expressed in temperature intervals (e.g., 1 K), most of the mass actually sits in the outer low-density envelope where the dust temperature is ≥ 13−15 K. This shows that, despite of the decreased sensitivity and increased uncertainty for the temperature and column density of the coldest dust, there is in practice no “hidden mass” problem in such cores. Nevertheless, the degeneracy between dust opacity and temperature cannot be solved with dust emission data alone and we can only derive an estimate of the related uncertainties in the derived column densities and masses.

thumbnail Fig. 12

Radial profiles of volume density, dust temperature, and relative mass per ΔR = 1000 au shell (left from top to bottom), as well as relative mass per logarithmic density interval and per temperature interval ΔT = 1 K (right) for B 68 (filled squares), CB 27 (open squares) and CB 244-SMM2 (filled triangles).

Open with DEXTER

Furthermore, the NIR–to–FIR opacity ratio is another key parameter of the dust opacity law. Although this ratio does not affect the derivation of a temperature from the thermal SEDs, it greatly affects the shielding of the cloud from the ISRF and the subsequent energy balance, which in turn affects the equilibrium temperature distribution predicted by self-consistent radiative transfer models and thus the ISRF strength deduced. This NIR–to–FIR opacity ratio is very uncertain, and the aforementioned opacity models have large scatter; they vary within a factor of four. Since different models have different trends when going from ISM-type to more processed dust, we cannot easily predict in which direction the radiative transfer temperature profiles would change and how the comparison with the ray-tracing results would be affected. Consequently, we cannot exclude currently that our conclusions on the value of (Sect. 4.3), which we derive from comparing the radiative transfer temperature profiles with the ray-tracing results, might be significantly affected by this uncertainty in the NIR–to–FIR opacity ratio of the applicable dust model. Thus, we need better observational constraints on the NIR–to–FIR opacity ratio, a subject which we will analyze and discuss in more detail in a forthcoming paper.

Related to this and continuing in the same theme is the uncertainty in the applicable albedos (or scattering efficiencies) mentioned in Sect. 3.2. The difference between the applied WD3.1 albedos and the possible alternative WD5.5B (Weingartner & Draine 2001) would be an about 60% increased extinction efficiency at the peak of the albedo spectrum at λ ≈ 1 μm. Although this is much smaller than the aforementioned general uncertainty in the NIR–to–FIR opacity ratio, the application of WD5.5B albedos instead of those from WD3.1 would result in systematically slightly increased shielding of the cold core centers in the radiative transfer modeling and thus imply the derivation of even higher values of sISRF (see also discussion in Sect. 5.1.4.)

Finally, laboratory experiments indicate that the dust opacities in the FIR and mm-range have a more complicated behavior than assumed by current dust models. Boudet et al. (2005) and Coupeaud et al. (2011) showed that the spectral index β in the FIR and mm-range depends on temperature and on wavelength. Both studies found an increase of β above two at low temperatures and at wavelengths above a few hundred μm, depending on the composition of the grains. These findings can be explained with a quantum mechanical model of amorphous grains (Paradis et al. 2011). Since robust observational constraints on such a (predicted) behavior are still lacking, we do not take into account this effect, but only mention here that it could contribute additional uncertainties. Although the changes would probably be small compared to other uncertainties, this effect would slightly lower the central core temperature inferred from the SED (by at most a few tenths of K), slightly increase the local core density, and also slightly affect the energy balance in the self-consistent radiative transfer modeling.

Table 4

Comparison of NH peak column densities with previous dust emission studies.

5.1.4. Interstellar radiation field

Not only the cloud structure is asymmetric, the IRSF can also be. It could be enhanced towards the direction of luminous stars, star-forming regions, or the galactic plane. On the other hand, large molecular clouds located between the galactic plane and the globules can shield the globules from the general ISRF. In fact, it was found in Launhardt et al. (2013) that most globules have azimuthally varying outer temperatures of their envelopes with one side often being warmer than the other. However, Launhardt et al. (2013) did in general not find a connection between directions of increased temperatures and galactic heating sources. In our modeling approach, we should not introduce significant systematic errors into the analysis by assuming an isotropic ISRF, because local variations of the ISRF should be coupled to local variations of the dust temperatures and thus an azimuthally averaged temperature profile should to first order correspond to the azimuthally averaged ISRF.

We also estimated the dependence of the results on the spectral shape of the ISRF by making the same analysis as described in Sect. 4.3 for radiative transfer models that assume the ISRF model by Black (1994, see also Fig. 4. We find that sISRF is by a factor 1.3−1.5 higher for these models than for models assuming the GALPROP ISRF, which roughly corresponds to the difference in the flux at wavelengths shortward of 8 μm (Sect. 3.4 and Fig. 4). The relative shape of the temperature profile (in–out contrast) is nearly unaffected. Although the comparison with the Black (1994) ISRF comprises a very simple and limited test of the effects from varying the spectral shape of the ISRF, we conclude that the uncertainty in our knowledge of the exact spectral shape of the ISRF is negligible compared to the uncertainty of its total flux at UV to IR wavelengths. Furthermore, as discussed in the previous section, we do not consider our results on the absolute strength of the ISRF robust since the degeneracy with the observationally only poorly constrained NIR–to–FIR dust opacity ratio cannot be resolved based on the available data.

5.2. Comparison with previous studies

In Table 4, we compiled a comparison of NH peak column densities (i.e., NH2 converted where applicable) derived from this work with earlier studies of the dust emission from the respective sources. We note that, in addition to the differences in wavelength coverage, data quality, and treatment of the dust temperature, the studies also differ in angular resolution (12′′36′′) and dust opacity law used (see last column in Table 4). Nevertheless, the comparison shows that, with the few exceptions discussed below, the results of the different studies are consistent with each other if one considers the methods involved, that is, simple single-temperature SED fits result in somewhat lower column densities than the ray-tracing inversion and assumptions of a too low fixed temperature results in higher column densities. The ≈40% higher column density derived by Schmalzl et al. (2014) for CB 17 can be attributed entirely to the different angular resolutions that were used in these two studies (25′′ vs. 36′′; see also discussion in Sect. 5.1.2). The higher column densities derived for B 68 by Nielbock et al. (2012) and Lippok et al. (2013) are not related to a smaller beam size, but can only be attributed to the problems of fine-tuning the profile parameters when the steep radial temperature gradient is not fully resolved (see discussion in Sect. 5.1.2) and the improved PoS – LoS convergence scheme for the ray-tracing algorithm mentioned in Sect. 4.1 and used for this paper.

To our knowledge, B 68 is the only core in our sample with a previous characterization of its dust temperature structure by other authors. Bianchi et al. (2003) compared SCUBA 850 μm and SIMBA 1.2 mm maps of this globule to the NIR extinction map of Alves et al. (2001b) and calculated the dust temperature distribution with a model by Gonçalves et al. (2004). The correlation of the dust emission to the NIR extinction showed a flattening toward large AK and they showed that models assuming a positive temperature gradient can describe the data better than models assuming an isothermal cloud, although their data did not constrain the temperature of B 68. Bergin et al. (2006) also calculated the dust temperature distribution of B 68 using the radiative transfer code of Zucconi et al. (2001) and the density structure that was derived by Alves et al. (2001b) from NIR extinction measurements. They deduced a central dust temperature of about 8 K, which agrees well with our analysis. Recently, Roy et al. (2014) used an inverse-Abel transformation based technique to infer the dust temperature and density structure of B 68 from Herschel data. Apart from the treatment of the outer halo (Sect. 3.3 and Eq. (1)), and the use nH2 instead of nH, their results agree well with ours from this paper. In Nielbock et al. (2012), we also used radiative transfer modeling to test the hypothesis that the increased FIR emission on the side facing the galactic plane (see Fig. C.8 in Launhardt et al. 2013) could be the result of a stronger ISRF coming from this direction. While the models supported this assumption qualitatively, they could not constrain the degree of anisotropy in the ISRF.

Another study that investigated the density and temperature structures of 20 dense cores in the L 1495 cloud in Taurus based on Herschel data was recently presented by Marsh et al. (2014). Although there is no overlap between their sources and ours, their COREFIT algorithm differs from our ray-tracing inversion in several aspects (e.g., forward-modeling vs. inversion starting from data, assumption of spherical geometry vs. allowing for moderate deviations from even elliptical geometry, etc.), and they use a different model for the ISRF from us, these authors derive conclusions that are very similar to ours. They estimate central core dust temperatures in the range 6−12 K (our range spans 7.5−11.9 K) that are also negatively correlated with peak column density, suggesting external heating by the ISRF and dust shielding of the central cores. These authors also encounter the problem that radiative transfer models require un up-scaling of the ISRF with respect to the COREFIT results to explain their central core temperatures. However, due to the different modeling approaches and ISRF models involved, a quantitative comparison with our results is difficult. As in our paper, they don’t find a conclusive explanation for this behavior, but can only suggest that further study is necessary.

Last but not least, we want to mention that other authors have also attempted to derive gas kinetic temperature distributions of similar starless cores. For example, Pagani et al. (2007) employed 1D non-LTE molecular line radiative transfer modeling to infer the radial profile of the kinetic gas temperature in the L 183 starless core from observations of N2H+ and N2D+. They find that that the gas in the core center is very cold (7 ± 1 K) and thermalized with the dust and increases outward up to about 12 K at 1.3 × 104 au. While both the trend of outward increasing temperature as well as the absolute values of the temperature compare well to our results on the dust temperature in our cores, a direct comparison is difficult for the following reasons. Due to the abundance and excitation profiles of the respective molecules (see also Lippok et al. 2013), the sensitivity at larger radii, where we still can measure dust temperatures, is lost. Furthermore, at larger radii and lower densities, the dust and gas are no longer thermally coupled and the gas temperature is expected to rise significantly above the dust temperature.

6. Summary and conclusions

With the Herschel space observatory, it has become possible to spatially resolve nearby starless cores in the FIR and to spectrally sample the peak of the SED of the thermal dust radiation with high sensitivity for the first time. Such observations allow us to break the density-temperature degeneracy in the interpretation of their dust continuum emission and to derive tight constraints on the dust temperature and density distribution of the cores. For this purpose, we have developed a ray-tracing inversion technique with which we can reconstruct the 3D temperature and density structure of starless cores directly from the data with a minimum of assumptions and without restriction to a specific physical model (Nielbock et al. 2012; Lippok et al. 2013; Schmalzl et al. 2014).

In this paper, we use Herschel observations of six starless cores from the EPoS sample to compare the observed temperature profiles from our ray-tracing inversion technique to the equilibrium temperature profiles from self-consistent 1D radiative transfer models. The only free parameters of the radiative transfer models are the relative strength of the ISRF (which we scale freely) and the selective extinction of the ISRF by a tenuous envelope (which we vary within the observational constraints). We restrict our analysis to one specific dust opacity model (OH5a), which we consider the most physically meaningful dust model for the core interiors, and only discuss how the application of other dust models would affect the results.

We derive physical outer radii of the six globules within the range (6.5 ± 2.5) × 104 au, total core masses within a fixed radius of 5 × 104 au in the range 2.6−14 M, central volume H number densities in the range (2.5−32) × 104 cm-3, peak column densities in the range (1.0−9.2) × 1022 cm-2, and density profiles that are characterized by a flat-density core with radii in the range (6−17) × 103 au, power-law fall-offs that are typically somewhat steeper than BES, and extended tenuous outer halos with typical dsnities of a few hundred H cm-3. We show that most of the core mass is neither located in the cold core centers, nor in the extended tenuous envelopes, but within a shell with ≈5000−15 000 au radius, where the dust temperature is ≈2−5 K higher than at the core centers and the densities are a few 104 cm-3. All starless cores are found to be significantly colder inside than outside. Central core temperatures are in the range 7.5−11.9 K and show a strong negative correlation with peak column density. Outer envelope temperatures are in the range 13.5−19 K and core – envelope temperature differences are in the range 2.4−9.6 K. Altogether this suggests that the thermal structure of the cores is dominated by external heating from the ISRF and shielding by dusty envelopes.

We find that for four of the six cores, the dust temperature distributions inferred directly from the dust emission data with the ray-tracing inversion method can be reproduced well with self-consistent radiative transfer models. The best agreement is achieved for relatively round sources without nearby (blending) protostars, for which both the ray-tracing inversion and 1D radiative transfer work best (CB 4 and B 68). Slight discrepancies in the detailed temperature profile shapes, albeit within the uncertainty ranges, are encountered for sources that are very elliptical (CB 27) or have nearby protostellar cores that partially overlap (in projection) with the starless core (CB 244). In the first case, the 1D approximation in the radiative transfer model was not optimal and it remains to be shown if the temperature structure inferred with the ray-tracing inversion can be reproduced better with 3D radiative transfer. In the second case, the blending by a nearby protostellar core and imperfect masking leads to larger uncertainties in the derived flux density profiles.

For two cores, the discrepancies between the temperature profiles inferred with the ray-tracing inversion and 1D radiative transfer exceed the formal uncertainties. Of these, CB 17 was likely affected by strong blending by a nearby protostar, resulting in large uncertainties in the reconstructed flux density profiles of the starless core. For the other core, CB 26, we found that this object may be a super-projection of three or more filamentary cores, such that neither 1D radiative transfer is applicable, nor does the ray-tracing inversion result in a reliable temperature profile.

We also confirm preliminary earlier results from various other studies which found that the usually adopted canonical value of the total strength of the ISRF in the solar neighbourhood is too low when an OH5-type dust model is invoked. This becomes evident by radiative transfer models providing core temperatures that are lower than observed if the ISRF is not increased. However, since the dust opacity model parameter that most strongly affects the conclusion on , the NIR–to–FIR dust opacity ratio, is poorly constrained observationally, we cannot resolve this degeneracy and are unable to draw robust conclusions on the actual strength of the local ISRF.

In summary, we conclude that for cores with not too complex morphology and no blending by additional nearby sources, the ray-tracing inversion technique infers temperature and density profiles that can be well-reproduced with self-consistent radiative transfer models. Moreover, the ray-tracing inversion technique can naturally account for moderate deviations from spherical symmetry. Hence the method can also be applied to other cores, provided their geometry is not too complex and the data are of similar quality and cover a similar wavelength range to those in this study.


1

Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

2

We note that this opacity model is slightly different from the often-used OH5 model which assumes coagulation at a density of 106 cm-3 (which is higher than what we actually observe in the cores). The “OH5a” model is not tabulated in the Ossenkopf & Henning (1994) paper, but is available at ftp://cdsarc.u-strasbg.fr/pub/cats/J/A+A/291/943

3

Stellar photons that are scattered at small dust grains in the optically thin halos (e.g., Foster & Goodman 2006).

4

Ignoring the very low-mass secondary core at the tip of the south-western trunk in B 68 (see Alves et al. 2001a; Nielbock et al. 2012).

Acknowledgments

We thank A. Strong and T. Porter for advice regarding the details of the GALPROP ISRF model used in this paper. We also thank A. M. Stutz, Y. L. Shirley, J. Steinacker, and E. Keto for helpful discussions and comments on the paper draft. The work of J.K. and S.E.R. was supported by the Deutsche Forschungsgemeinschaft priority program 1573 (“Physics of the Interstellar Medium”). H.L. and M.N. were funded by the Deutsches Zentrum für Luft- und Raumfahrt (DLR). PACS has been developed by a consortium of institutes led by MPE (Germany), including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).

References

All Tables

Table 1

Source list.

Table 2

Properties of the studied cores.

Table 3

Best-fit values of sISRF and NH(rsym).

Table 4

Comparison of NH peak column densities with previous dust emission studies.

All Figures

thumbnail Fig. 1

Dust opacity models (extinction/absorption coefficient per g of dust) used in this paper. Black lines show the absorption (solid) and extinction (dashed) coefficients of the modified OH5a model (Ossenkopf & Henning 1994). The gray line shows the OH1 model (see Sect. 5.1.3; for clarity, only the absorption coefficient spectrum is shown). Vertical dashed lines mark the wavelength range of the original OH models, i.e., values outside this range are extrapolated (see text). Arrows indicate the wavelength bands of the dust emission data used in this paper.

Open with DEXTER
In the text
thumbnail Fig. 2

Modified normalized Plummer profiles (Eq. (1)) with parameters chosen to mimic: (1) a simple power law; (2) a flat-density core with smooth outer edge; (3) an isothermal BES; and (4) a profile like observed in B 68.

Open with DEXTER
In the text
thumbnail Fig. 3

Comparison of the self-consistently calculated equilibrium temperature profile for B 68 (solid line) and an analytical profile fit following Eq. (2). The little kinks in the radiative transfer temperature curve at r ≈ 2000 and 7000 au are computational artefacts due to the high optical depth.

Open with DEXTER
In the text
thumbnail Fig. 4

ISRF model developed by the GALPROP consortium for the solar neighbourhood (Ackermann et al. 2012) with the CMB already included (solid line). For comparison we also show the model by Black (1994) as parameterized in Zucconi et al. (2001).

Open with DEXTER
In the text
thumbnail Fig. 5

Mid-plane dust temperature maps (color, all maps on the same scale) of all six globules studied in this paper, derived with the ray-tracing inversion technique and overlaid with contours of the hydrogen column density (white: 1021 and 1022 cm-2, black: 0.5 and 5 × 1021 cm-2). Yellow squares indicate the center positions of the starless cores as listed in Table 1 and in Launhardt et al. (2013). Asterisks indicate the positions of embedded protostars (see Launhardt et al. 2013).

Open with DEXTER
In the text
thumbnail Fig. 6

Radial profiles of all six starless cores, derived with the ray-tracing inversion technique and plotted up to rsym. a) Mid-plane hydrogen volume number density vs. radius; b) hydrogen column density vs. impact radius; c) mid-plane dust temperature vs. radius.

Open with DEXTER
In the text
thumbnail Fig. 7

Cumulative radial mass distribution of the six starless cores derived with the ray-tracing inversion technique and after masking out the nearby protostars in CB 17, CB26, and CB 244.

Open with DEXTER
In the text
thumbnail Fig. 8

Central density vs. total gas mass for the six starless cores (data from Table 2). The dashed line indicates the maximum stable density of a pressure-supported, self-gravitating modified (nonisothermal) BES as calculated by Keto & Caselli (2008, their model with photoelectric heating at the core boundary taken into account). Error bars represent the mean relative uncertainties of 30% on the cloud mass and 50% on the central density (see Sect. 5.1).

Open with DEXTER
In the text
thumbnail Fig. 9

Effect on the radial equilibrium temperature distribution predicted for B 68 from scaling the total flux of the GALPROP ISRF by factors sISRF = 0.5,1,2,3 (with OH5a dust and NH(rsym) = const. = 1.5 × 1021 H cm-3; black lines). Blue dots show the mid-plane dust temperature distribution inferred with the ray-tracing inversion.

Open with DEXTER
In the text
thumbnail Fig. 10

Effect on the radial equilibrium temperature distribution predicted for B 68 from varying the envelope extinction of B 68 via NH(rsym) for values of 0.4,0.8,1.7,2.5,4.2 × 1021 H cm-3 (with OH5a dust and sISRF = const. = 2.2; black lines). Blue dots show the mid-plane dust temperature distribution inferred with the ray-tracing inversion.

Open with DEXTER
In the text
thumbnail Fig. 11

Comparison of the temperature distributions derived with the ray-tracing inversion technique (blue dots and shaded regions) to the best-matching thermal equilibrium temperature distributions (black lines). Blue dots represent pixel values of the mid-plane temperature map derived with the ray-tracing inversion method. Shaded regions indicate the uncertainty of the mean temperature profile which we estimate as K (Sect. 5.1). All models assume OH5a dust. Best-matching scaling factors sIRSF for the total flux of the ISRF are also indicated in each panel.

Open with DEXTER
In the text
thumbnail Fig. 12

Radial profiles of volume density, dust temperature, and relative mass per ΔR = 1000 au shell (left from top to bottom), as well as relative mass per logarithmic density interval and per temperature interval ΔT = 1 K (right) for B 68 (filled squares), CB 27 (open squares) and CB 244-SMM2 (filled triangles).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.