© ESO, 2016

1. Introduction

Significant technological advances in the astronomical instrumentation during the last four decades enabled measurements of the far-infrared thermal dust emission (usually optically thin in that wavelength range) and hence estimates of the masses of dusty objects. Fitting the far-infrared and submillimeter flux or intensity distributions of optically thin sources can give their average temperatures and masses (Hildebrand 1983). This simple method has become standard in studies of Galactic star formation and a major source of our knowledge of the physical properties and evolution of self-gravitating cores and protostars. Although there are more sophisticated approaches (e.g., Kelly et al. 2012), a simple fitting of the observed spectral shapes remains the most widely used method in the observational studies of star formation (e.g., Könyves et al. 2015). Its inaccuracies, biases, and limitations need to be carefully investigated before reliable conclusions can be made on the physical properties and evolution of the observed objects.

Mass derivation from fitting total fluxes or pixel intensities involves a strong assumption of a constant temperature within an object. In addition to such poorly known parameters as the distance, the far-infrared opacity and its power-law slope, and the dust-to-gas mass ratio, the most problematic assumption is that a single color temperature obtained from the fitting is a good approximation of the mass-averaged physical dust temperatures. This may be true for only the simplest case of the lowest density starless cores, transparent in the visible wavelength range and thus practically isothermal, but it is clearly invalid for the protostellar envelopes that are centrally heated by accretion luminosity. Very sensitive dependence of the emission of dust grains on their temperature warrants careful investigation of the effects of nonuniform temperatures. There are papers that have investigated some aspects of the problem, notably the correlation between the estimated temperatures and power-law opacity slopes (e.g., Shetty et al. 2009a,b; Juvela & Ysard 2012, and references therein) and the inaccuracies of mass derivation and their effect on the resulting core mass function (Malinen et al. 2011).

The present purely model-based study simplifies the problem by removing the “observational layer” between the physical reality and observers. To investigate the intrinsic effects related to the physical objects, all observational intricacies (the complex filamentary backgrounds, instrumental noise, calibration errors, different angular resolutions across wavebands, etc.) were explicitly ignored. Measurement errors in intensities and fluxes are assumed to be nonexistent and the radiation emitted by the model objects is known to a high precision, limited only by their numerical accuracy. A grid of radiative transfer models of starless cores and protostellar envelopes was computed and their total fluxes and image intensities were fitted to derive the model masses. Known true values of the numerical models allow us to assess the qualities of the methods and fitting models, as well as the effects of nonuniform temperatures, far-infrared opacity slope, selected subsets of wavelengths, background subtraction, and angular resolutions. The main goal was to quantify how much the mass derivation methods are affected, what the realistic uncertainties of the temperatures and masses are, and what one could possibly do to improve the estimates. Although this study is completely independent of the instruments and wavebands used in actual observations, it employs six Herschel wavebands (70, 100, 160, 250, 350, and 500 μm; Pilbratt et al. 2010), for which a wealth of recent results in star formation has been obtained (e.g., Könyves et al. 2015, and references therein).

The radiative transfer models of starless cores and protostellar envelopes are presented in Sect. 2, the methods of mass derivation from fitting far-infrared and submillimeter observations are introduced in Sect. 3, the results of this work are presented in Sect. 4 and discussed in Sect. 5, the conclusions are outlined in Sect. 6, and further details are found in Appendices A–E.

2. Radiative transfer models

The models were computed with the 3D Monte Carlo radiative transfer code RADMC-3D by C. Dullemond¹. Spherical model geometry was chosen to simplify the problem by reducing the number of free parameters involved in the study: asymmetries in model density distribution would introduce dependence on viewing angle (e.g., Men’shchikov & Henning 1997; Men’shchikov et al. 1999; Stamatellos et al. 2004) and hence increase the uncertainties of derived parameters. Isotropic scattering by dust grains was considered.

Grids of models for starless cores and protostellar envelopes were constructed, covering the ranges of masses M (0.03−30M_⊙) and luminosities L (0.03−30 L_⊙) relevant for both low- and intermediate-mass star formation. The masses and luminosities were sampled at the values of 0.0316, 0.1, 0.316, 1, 3.16, 10, 31.6 (separated by a factor of ); for simplicity, they will be referred to as 0.03,0.1,0.3,1,3,10,30 (M_⊙, L_⊙). Although the luminosity of an accreting protostar depends on its mass, the goal is to separate the effects of masses and luminosities.

In addition to isolated models, their embedded variants were constructed by implanting the isolated models into the centers of larger spherical background shells of uniform densities, in order to simulate the fact that stars form within their dense parental clouds that shield the embedded objects from the interstellar radiation field. All models were put at a distance D = 140 pc of the nearest star-forming regions.

2.1. Dust properties

Properties of the real astrophysical dust grains are poorly known and they are unlikely to be universal in the different star-forming regions observed. The standard mass derivation methods ignore many complications related to the cosmic dust grains, assuming just a simple power-law opacity across all bands being fitted. For example, the presence of very small, stochastically heated grains is neglected (e.g., Desert et al. 1990); the contribution of these grains to the emission of starless cores and protostellar envelopes can become significant at λ ≲ 100 μm (e.g., Bernard et al. 1992; Siebenmorgen et al. 1992). For consistency with the mass derivation methods and previous studies of star formation, this model study adopts tabulated absorption opacities κ_abs for grains with thin ice mantles (Ossenkopf & Henning 1994), corresponding to coagulation time t = 10⁵ yr and number density n_H = 10⁶ cm^-3 (Fig. 2).

The opacity values at long wavelengths λ> 70 μm were replaced with a power law κ_λ ∝ λ^-2; the modification aimed at testing the widely used assumption on the power-law far-infrared opacities . At short wavelengths (0.1<λ< 1 μm), the opacities were extrapolated with a power law κ_λ ∝ λ^-0.87 based on the last tabulated values. Although dust scattering is unimportant in the far-infrared, scattering opacities were constructed to resemble the values and wavelength dependence κ_sca ∝ λ^-4 of typical dust grains. The resulting dust opacity at λ> 70 μm was parameterized by κ₀ = 9.31 cm² g^-1 (per gram of dust), λ₀ = 300 μm, and β = 2, with the maximum opacities limited by 10⁵ cm² g^-1 (Fig. 1).

2.2. Density distributions

The density structure of starless cores was approximated by an isothermal Bonnor-Ebert sphere (Bonnor 1956) with a temperature of 7 K and a central density of 5.2 × 10^-18 g cm^-3. This somewhat arbitrary choice of ρ(r) gives just a simple and convenient functional form (Fig. 2) resembling the observed flat-topped density profiles of starless cores (e.g., Alves et al. 2001; Evans et al. 2001). The issue of the gravitational instability (or stability) of the model cores is irrelevant for this study of the mass derivation methods. Protostellar envelopes were modeled as infalling spherical envelopes with the power-law densities ρ(r) ∝ r^-2 (e.g., Larson 1969; Shu 1977) around a central source of accretion energy (Fig. 2).

Model dust densities were scaled to obtain the desired grid of masses 0.03,0.1,0.3,1,3,10, and 30 M_⊙ using the standard dust-to-gas mass ratio η = 0.01. The outer boundary of all the models was placed at the same distance of R = 10⁴ AU, beyond which their density either changed to zero (isolated models) or remained constant until R_E = 3 × R (embedded models). The embedding cloud density was set equal to ρ(R) (Fig. 2), which corresponds to the denser models (i.e., more massive) being formed in a denser environment. Most of the mass of the model starless cores and protostellar envelopes (96% and 90%, respectively) is contained in their outer parts (0.1 R<r<R).

For the starless cores, the inner boundary was arbitrarily set to R₀ = 50 AU, as their densities are essentially constant and hence do not need to be resolved at smaller radii. The inner boundary of the dusty protostellar envelopes is defined by the dust sublimation temperature T_S ~ 10³ K. An exact value of T_S depends on the chemical composition and sizes of dust grains and so does the radius R₀ of the inner dust-free cavity. For the purpose of this study, it is adequate to adopt a single value T_S = 10³ K. With the model κ_ν and ρ(r) (Sect. 2.1), the resulting radiative-equilibrium temperatures (Fig. 3) lead to the inner boundaries of the dusty protostellar envelopes that are fairly accurately described by a simple formula, (1)The model space between the inner and outer boundaries was discretized by nonuniform grids with the relative zone sizes δlog r that smoothly varied from 0.6 to 0.02 (100−150 zones) for starless cores and from 0.002 to 0.06 (200−300 zones) for protostellar envelopes.

2.3. Radiation sources and optical depths

All models were illuminated from the outside by an isotropic interstellar radiation field (Black 1994) with the “strength” parameter G₀ = 1 (e.g., Parravano et al. 2003). The bolometric luminosity of the interstellar radiation field (ISRF) entering the isolated models at R amounted to L_ISRF = 1 L_⊙, whereas that crossing the boundary of embedding clouds at R_E was 9 L_⊙.

In addition to the external radiation field, the models of protostellar envelopes were assumed to be heated at their centers by a blackbody source of luminosity L_⋆ of 0.03,0.1,0.3,1,3,10, and 30 L_⊙ with an effective temperature of T_⋆ = 5770 K. Actual values of T_⋆ are unimportant, as the sources of accretion energy are surrounded by the completely opaque dusty envelopes reprocessing the hot radiation to T ≲ 10³ K very deep in their interiors.

Distribution of optical depths within dusty envelopes is one of the main parameters (along with the density structure) for the transfer of radiation and resulting radiative-equilibrium temperatures. All models are quite opaque at visible wavelengths, with radial optical depths τ_V ≈ 3−3 × 10³ for starless cores and τ_V ≈ 500−6 × 10⁵ for protostellar envelopes of different masses and luminosities (Fig. 1). At the far-infrared wavelength of 100 μm, starless cores with M< 3 M_⊙ are transparent, whereas the ones with M ≳ 3 M_⊙ are optically thick towards their centers. All protostellar envelopes are optically thick at 100 μm and some of them (with M ≳ 0.3 M_⊙) are even opaque at 500 μm towards their centers. Sizes of the dust-free cavities of protostellar envelopes increase with the luminosity of their central energy sources (cf. Eq. (1), Fig. 2), thus the optical depths of the envelopes decrease, approximately as .

High far-infrared optical depths of the model starless cores and protostellar envelopes are localized within relatively small spherical zones around their centers. Angular radii of the opaque dusty zones in the protostellar models can be described (at 70−500 μm) by a simple empirical expression (2)which can also be used (within a factor of 1.5−2) for the high-mass models of starless cores of 3,10, and 30 M_⊙, in which the opaque zone exists only at λ ≤ 70,160, and 250 μm, respectively.

The density profiles ρ(r) ∝ r^-2 of the protostellar envelopes are similar to those of the starless cores for r ≳ 0.1 R (Fig. 2). Therefore, whenever an inner opaque zone exists in the objects, its mass obeys m(ϑ) ∝ ϑ, hence the fractional mass is the fractional radius and, using Eq. (2), can be written as (3)where ϑ, D, and R are in units of arcsec, pc, and AU, respectively. At λ ≲ 100 μm, the opaque zone of high-mass objects extends over a large fraction of their mass. This means that the standard assumption of the far-infrared transparency is severely violated for massive objects. Protostellar envelopes with M< 3 M_⊙ have small opaque zones that contain little mass and thus they cannot substantially affect the standard methods of mass derivation.

2.4. Temperature distributions

The models of starless cores and protostellar envelopes acquire radiative-equilibrium dust temperatures T_d(r) shown in Fig. 3. In the adopted isotropic ISRF, the radiative-equilibrium temperature of dust grains with the model opacities from Fig. 1 is T_d = 17.4 K, the value that the isolated models and embedding clouds acquire at their outer boundaries in the limit τ_λ → 0. The lower mass models of starless cores are transparent and thus almost isothermal. Their higher mass counterparts develop steeper temperature gradients under the outer boundaries of the isolated models and embedding clouds and lower temperatures in their interiors (Fig. 3).

Displaying the same behavior under their outer boundaries, protostellar envelopes of all masses develop steep temperature gradients towards the inner boundary (Fig. 3). Higher accretion luminosities make the dust hotter and thus, with the adopted dust sublimation temperature T_S = 10³ K, the boundary of the inner dust-free cavity shifts towards larger radial distances (cf. Eq. (1)). An analytical approximation of the profiles T_d(r) for protostellar envelopes can be found in Appendix A.

Differences between the isolated and embedded models are highlighted by their different temperature distributions at the outer model boundary (Fig. 3). The temperatures of embedded models at r = R are significantly lower than those of the isolated models, owing to the absorption of ISRF in the embedding clouds (R<r ≤ R_E). The denser the embedding cloud is, the lower T_d(R) is and the greater the contrast to the isolated model (Fig. 3). As the bulk of the mass of the models is contained in the outer parts, the differences in the temperature profiles between the isolated and embedded models can greatly affect their observational properties, such as the images and total (integrated) fluxes.

2.5. Spectral energy distributions

After computing the self-consistent radiative-equilibrium dust temperature distributions T_d(r) from the radiative transfer models, observables – such as the intensity maps ℐ_ν and total fluxes F_ν – were obtained by a ray-tracing algorithm in separate runs of RADMC-3D. Effects of the Monte Carlo noise on F_ν, evaluated from the standard deviations about the azimuthally averaged intensity profiles I_ν(ϑ), are below 0.003% and 3% for the starless cores and protostellar envelopes, respectively, in all models and wavebands.

To emulate the standard observational procedure of flux measurements, F_ν were integrated from the background-subtracted model images ℐ_ν. The model background was evaluated as an average intensity within a δR′-wide annulus placed just outside the outer model boundary (R′>R). In practice, the annulus was one pixel in width (δR′ = 0.47″) and it was detached from the boundary by one additional pixel. For the isolated models, is the intensity of the isotropic ISRF, whereas for the embedded models, is determined by both and the transfer of radiation in the background cloud (R′ ≤ r ≤ R_E) along the rays passing through the annulus. Inaccuracies inherent in the standard algorithm of background subtraction are discussed in Sect. 5.2 and Appendix B.

Spectral energy distributions (SEDs) of the models of starless cores and protostars are shown in Fig. 4. The SED shapes depend on the density and temperature distributions (Figs. 2 and 3). Large differences between the SEDs for the isolated and embedded cores are mainly caused by differences in their temperature profiles near the model boundary. The SEDs of protostellar envelopes are affected by the same effects to a much lesser degree as their density profiles are centrally peaked and their temperature profiles are dominated by the internal radiation source. The SEDs of the models of different masses and luminosities show a large variety of shapes in the far-infrared domain (Fig. 4) due to varying optical depths and temperatures.

Providing a useful reference in our analysis, additional ray-tracing runs of RADMC-3D computed the fluxes of isothermal models. These are the same models described above (Fig. 2), in which self-consistent temperature profiles (Fig. 3) have been replaced with their mass-averaged values: (4)The resulting total fluxes of the isothermal models are denoted F_ν(T_M).

3. Fitting source fluxes and intensities

In observational studies, after obtaining multiwavelength images ℐ_ν and integrating background-subtracted (and deblended) fluxes F_ν of extracted sources, their spectral distributions need to be fitted to derive fundamental physical parameters, such as the source mass and luminosity.

The standard technique uses the well-known formal solution of the radiative transfer equation that can be written as (5)where I_ν is the observed specific intensity, T is the homogeneous temperature of an object, B_ν(T) is the blackbody intensity, and τ_ν is the optical depth of the object. After obtaining an image ℐ_ν ≡ I_νij, the total flux F_ν = ^∫ℐ_ν dΩ can be integrated over the solid angle Ω subtended by the object. For constant intensity, it reduces to F_ν = I_ν Ω. A critical assumption used in the derivation of Eq. (5) is that the object is homogeneous in temperature, whereas the temperatures of the astrophysical objects are actually nonuniform (cf. Fig. 3).

Two methods and two fitting models were explored in this work that have been used in observational studies of star formation to estimate source temperatures and masses.

3.1. Fitting total fluxes F_ν

In this method, the total fluxes F_ν are integrated from background-subtracted and deblended images ℐ_ν of source intensities and then are fitted to estimate source mass as one of the fitting parameters. With the adopted parameterization of the power-law opacity κ_ν = κ₀⁽ν/ν₀^)β, it is possible to write Eq. (5) in the form (6)where η is the dust-to-gas mass ratio and D is the source distance. The fitting model of Eq. (6) with five parameters (T, M, β, D, Ω) is referred to as modbody in this paper. After fitting F_ν and estimating the model parameters, the average column density N_H2 can be obtained from M = μm_HN_H2D²Ω, where μ = 2.8 is the mean molecular weight per H₂ molecule and m_H is the hydrogen mass.

With an additional assumption that measured fluxes F_ν represent optically thin emission², Eq. (6) can be written as (7)The fitting model of Eq. (7) with four parameters (T, M, β, D) is referred to as thinbody in this paper. By the definition (τ_ν ≪ 1), it produces only fits with the modified blackbody shapes κ_νB_ν(T) that are scaled up or down, depending on M. Obviously, the modbody fits with τ_ν ≪ 1 produce the same shapes as the thinbody model does, whereas the modbody fits with τ_ν ≫ 1 resemble a blackbody B_ν(T). In the intermediate (semi-opaque) cases, the short-wavelength parts of the fitted curves can be described by B_ν(T) while morphing into κ_νB_ν(T) at long wavelengths where the radiation becomes optically thin. With more flexible shapes, modbody can give better fits of the data, but it does not necessarily lead to good estimates of temperatures and masses.

After fitting fluxes with a modbody or thinbody model, an estimate of T_F and the corresponding mass M_F are obtained. For the realistic objects with strongly nonuniform temperatures T_d(r) (Fig. 3), emerging fluxes F_ν are heavily distorted from the simple shapes of the fitting models, hence these models are inadequate and an estimate of T_F does not guarantee that M_F is close to the true mass M. For the purpose of obtaining accurate M_F ≈ M, it is necessary (but not sufficient) to have T_F ≈ T_M, i.e., it is possible to interpret T in Eq. (7) as the mass-averaged T_M from Eq. (4). In fact, assuming τ_ν ≪ 1 in the far-infrared, the observed fluxes contain emission of all dust grains, which is proportional to the mass of dust at different T_d in the entire volume of an object: (8)Equations (7) and (8) can immediately be combined into a definition of the mass-averaged intensity B_νM. Since B_ν(T) ∝ T in the Rayleigh-Jeans domain, the equations are also readily converted into T_M from Eq. (4). In the model objects studied here, differences between B_νM and B_ν(T_M) quickly become negligible beyond the peak wavelength of the latter (λ≳2 λ_peak). Therefore, T_M is fully consistent with the fitting models at long wavelengths.

3.2. Fitting image intensities I_ν

In this method, it is possible to fit pixel intensity distributions I_νij of the background-subtracted and deblended images ℐ_ν of a source³, to derive a map of its column densities and then the source mass . It is convenient to express the modbody and thinbody models from Eqs. (6) and (7) as functions of the pixel column density N_H2: In this formulation, both models have only three fitting parameters (T, N_H2, β) in contrast to the case where total fluxes are fitted (five parameters for modbody in Eq. (6) and four parameters for thinbody in Eq. (7); see Sect. 3.1). Furthermore, limited angular resolutions of real images makes the results of fitting I_ν depend sensitively on the degree to which a source is resolved.

For fully resolved sources, such as the model objects used in this work, relatively small pixels sample completely independent intensities from different rays. For progressively lower angular resolutions, intensities within a beam become increasingly blended together. Radiation with different temperatures gets mixed not only along the line of sight, but also in the transverse directions, in the plane of the sky. For unresolved objects with intrinsic temperature gradients, radiation from the entire object becomes heavily blended, leading to strong distortions of their spectral intensity distributions.

An important assumption used in the derivation of N_H2ij is a constant temperature T_d(x,y,z) along the lines of sight within a certain radial distance from a pixel (i,j). The distance depends on the angular resolution of images: for less resolved sources, temperatures from a larger environment of the pixel contribute to its intensity. With low optical depths τ_ν ≪ 1 in the far-infrared, emission is observed from the entire column of dust grains at (i,j) with different temperatures T_d(z) along the line of sight. The reasoning associated with Eq. (8) can be applied to show that T in Eq. (10) is consistent with the column-averaged temperature (11)A mass-averaged temperature, equivalent to that from Eq. (4), can be obtained as .

3.3. Variable and fixed parameters

In most studies, the opacity slope β has been kept fixed in the fitting process to reduce the number of free parameters and improve the robustness of derived parameters. Following this practice, Sect. 4 presents and discusses only the results of fitting with a fixed opacity slope. When fitting intensities I_ν with β fixed, the number of free variable parameters becomes γ = 2 for both thinbody and modbody models (T, N_H2). When fitting fluxes F_ν, distance D is also assigned a fixed value to further reduce the degrees of freedom, although astronomical distances are poorly known. The number of free variable parameters is thus γ = 2 for thinbody (T, M) and γ = 3 for modbody (T, M, Ω).

In practice, after measuring F_ν = ^∫I_ν dΩ, the solid angle Ω over which I_ν were integrated is known⁴ and its value can be fixed, reducing γ for modbody to two free variables (T, M). In this model study, one could also keep Ω = π (R/D)² constant, as the true values of R and D are known; however, modbody would then become completely equivalent to thinbody. Indeed, fixing Ω of transparent objects at accurate (or even overestimated) values means that the optical depths in Eq. (6) are very small (τ_ν = κ_νηMD^-2Ω^-1 ≪ 1), which effectively converts modbody into thinbody. Only when fixing strongly underestimated values Ω ≪ π (R/D)², the far-infrared τ_ν become large enough to produce any noticeable differences between modbody and thinbody. This work investigates qualities of two different models, hence Ω was allowed to vary in all modbody fits of F_ν.

When fitting pixel intensities I_ν instead of F_ν, the far-infrared τ_νij within most of the image pixels (i,j) are small, even for perfectly resolved sources. For poorly resolved or unresolved sources, radiation within the beams gets diluted and maximum values of τ_νij in the images become smaller. All models of starless cores and protostellar envelopes contain 96% and 90% of their masses, respectively, in their outer parts (r ≳ 0.1 R, Fig. 2). Intensities in the outer parts of the source images come mostly from the pixel columns of dust with τ_νij ≪ 1 in the far-infrared. Only in the models with M ≳ 3 M_⊙ do they become substantially affected by the radiation from the central opaque zone (Sect. 2.3). As a result, the masses derived from fitting I_ν are almost the same (within ~20%) for both fitting models and hence only the thinbody results are presented for this method.

3.4. Data points and their subsets

Fitting was executed for a set of the total model fluxes {F_λi} or pixel intensities {I_λi}(i = 1,2,...,6) at the Herschel wavelengths λ_i of 70, 100, 160, 250, 350, and 500 μm. In this model-based study, the intensities and fluxes of numerical models have essentially no measurement errors. It makes sense, however, to make their uncertainties resemble typical observational values, hence to get an idea of realistic inaccuracies of the estimated parameters (masses, temperatures). Before the fitting, the model intensities and fluxes were assigned an additional (optimistic) uncertainty of 15%, a value similar to the levels of calibration errors in real observations (e.g., with Herschel). The above uncertainties were associated with the exact data points to see how typical data uncertainties translate into the resulting error bars of the derived parameters. Extra uncertainties come from the fact that the dust-to-gas ratio η, reference opacity κ₀, and distance D, which are used in the fitting models (Sects. 3.1, 3.2) but held constant, are actually poorly known. Conservatively assuming that the quantities have random and independent uncertainties of 20%, the latter were added in quadrature to those of the derived masses, for the same purpose of obtaining the total resulting mass uncertainties.

To isolate the effects of temperature gradients in starless cores and protostellar envelopes (Fig. 3), the fitting was done for several subsets of data, removing some (or none) of the shortest-wavelength points from the fitting process. The data subsets are denoted Φ_n = {Y_λi}, where Y_λi is either F_λi or I_λi and n is the number of the longest wavelengths used in the fitting⁵. Fits of total fluxes were considered successful (acceptable) and their results are shown below, if χ² ≤ n−γ + 1, with the last term added to allow testing χ² for zero degrees of freedom (n = γ). Fits of image intensities were considered successful, if the same goodness condition was fulfilled in all pixels within an object. These results, as well as the somewhat less reliable results with n−γ + 1 <χ²< 10 in some pixels are presented below. Details of the fitting algorithm can be found in Appendix C.

4. Results

This section describes derived parameters for both starless cores and protostellar envelopes, obtained from acceptable fits for all subsets Φ_n (Sect. 3.4) for both modbody and thinbody (Sect. 3.1). Results for the isothermal models are presented in Appendix D. To evaluate the effects of the uncertain far-infrared opacity slope, results are shown for β = 2 used in the radiative transfer modeling and for two other β values (1.67, 2.4), differing from the true value by a factor of 1.2. Results obtained with variable fitting parameter β are described in Appendix E.

Masses derived from fitting images ℐ_ν of objects with temperature gradients must depend on their angular resolutions (Sect. 3.2). To investigate this effect, the model images with pixels of 0.47″ were convolved with Gaussian beams of 1, 36, and 144″ (FWHM) and then resampled to 1, 12, and 48″ pixels, respectively. For the objects with diameters of 142″ (2 × 10⁴ AU, Fig. 2), the three variants represent resolved, partially resolved, and unresolved cases.

In this paper, the term uncertainties refers to the error bars of measured or derived quantities, the term inaccuracies (sometimes simply errors) refers to the deviations of the derived quantities from their model values, and the term biases denotes variable systematic dependences of inaccuracies across the ranges of model parameters (M, L, T_M).

4.1. Selected examples

Examples of the fits of F_ν for the isolated starless cores and protostellar envelopes with masses of 0.03, 0.3, and 3 M_⊙ are shown in Fig. 5. Although the qualitatively similar plots for embedded models are not presented, their derived parameters and uncertainties are described in Sect. 4.2.

Flux distributions of the isolated starless cores (Fig. 4) are similar to those of the modified blackbodies κ_νB_ν(T_M). The fits for a low-mass core with M = 0.03 M_⊙ shown in Fig. 5 are identical for all subsets Φ_n since the core is nearly isothermal, with T_d(r) very similar to its T_M = 16.3 K. Fluxes of the higher-mass cores of 0.3 and 3 M_⊙ display larger deviations from the fluxes F_ν(T_M) of isothermal models for larger subsets Φ_n (n = 3 → 6). The shapes of F_ν become “hotter” because of the steeper temperature profiles (Sect. 2.4) at the outer boundary (Fig. 3). For a massive core of 3 M_⊙, discrepancies between F_ν and F_ν(T_M) at λ ≲ 160 μm reach factors ≳ 5.

Flux distributions of an isolated protostellar envelope with M = 0.03 M_⊙ (Fig. 4) display various shapes that are quite different from those of κ_νB_ν(T_M), whereas for a more opaque envelope of 3 M_⊙ they become similar to the modified blackbody shapes. The protostellar fits (Fig. 5) show greater deviations for larger subsets Φ_n (n = 3 → 6), much larger than those of starless cores. Differences between F_ν and F_ν(T_M) reach orders of magnitude at λ ≲ 100 μm. The shapes appear much “hotter” owing to T_d ~ 100−10³ K (Sect. 2.4) deep inside the envelopes (Fig. 3). The lower-mass protostellar envelopes are more transparent and the hot emission greatly distorts F_ν at λ ≲ 250 μm.

4.2. Properties derived from fitting fluxes F_ν

Isolated starless cores, β = 2 (Fig. 6). For the low-mass, transparent cores (M → 0.03 M_⊙), quite accurate values T_F ≈ T_M and M_F ≈ M are derived for all subsets Φ_n. For the denser, more opaque cores (M → 30 M_⊙), derived T_F and M_F become more over- and underestimated, respectively, as the spectral shapes of F_ν become much wider and distorted towards shorter wavelengths (Fig. 4). The biases and inaccuracy of the estimates depend on the subset Φ_n, with the least inaccurate T_F and M_F obtained for the thinbody fits of Φ₂. However, the biases of the parameters across the entire mass range remains fairly strong. Derived masses of the starless cores are underestimated within a factor of 2 for 1 <M ≤ 3 M_⊙ and factor of 5 for 3 <M ≤ 30 M_⊙.

Embedded starless cores, β = 2 (Fig. 6). For the low-mass, transparent cores (M → 0.03 M_⊙), M_F are underestimated by a factor of 1.35 for all subsets Φ_n, although T_F are quite accurate because the standard observational procedure of background subtraction ignores the fact that embedding backgrounds tend to be rim-brightened at their outer boundary R (Appendix B, Sect. 5.2). The embedded cores have T_d(r) that are quite flat across their boundary for all masses (Fig. 3). Having no flux distortions caused by nonuniform temperatures (Fig. 4), the F_ν peaks of the most massive cores (M → 30 M_⊙) move towards the longest wavelength (λ₆ = 500 μm), which leads to T_F and M_F that are under- and overestimated, respectively.

Isolated protostellar envelopes, β = 2 (Fig. 7). Emission of the hot dust with T_d ~ 100−10³ K greatly skews their F_ν towards shorter wavelengths (Fig. 4). This becomes especially significant for the lower mass, more transparent envelopes (M → 0.03 M_⊙, L_⋆ → 30 L_⊙) that produce hotter dust over a much larger volume (Fig. 3). The thinbody fits of larger subsets Φ_n (n = 3 → 6), lead to errors in T_F and M_F that reach factors of 1.6 and 3, respectively. The smallest subset Φ₂ is unaffected by the hot emission and it produces fairly accurate thinbody estimates of T_F and M_F (for all M and L_⋆) within factors of 1.1 and 1.3, respectively.

Embedded protostellar envelopes, β = 2 (Fig. 7). Results are qualitatively similar to those of the isolated envelopes, although with larger inaccuracies. Derived M_F are underestimated by at least a factor of 1.5, mostly due to over-subtraction of the rim-brightened embedding background (Appendix B, Sect. 5.2). Although the envelopes have T_d(r) that are quite flat across their boundaries (Fig. 3), their derived parameters are greatly affected by the skewed F_ν owing to the hot dust deep in their interiors. The thinbody fits of large subsets Φ_n (n = 3 → 6) lead to inaccuracies in T_F and M_F as large as factors of 1.4−1.8 and 3−5, respectively. The most accurate T_F and M_F, obtained for the smallest subset Φ₂, are underestimated within factors of 1.2 and 2.

Effects of the adopted opacity slope β on the estimated parameters are similar for both starless cores (Fig. 6) and protostellar envelopes (Fig. 7). Although detailed behavior of the differences with respect to the above results for true β = 2 depends on the subset Φ_n, clear general trends can be seen. Shallower slopes (β = 1.67) lead to an increase in T_F and thus M_F becomes smaller, whereas steeper slopes (β = 2.4) lead to a decrease in T_F and hence M_F becomes larger, in both cases by a factor of approximately 2.

The thinbody fitting model produces much better overall results than modbody does. Parameters estimated with modbody become so incorrect that they may be considered completely unusable. The importance of estimating accurate mass-averaged temperatures T_M for deriving correct masses M_F is illustrated by the isothermal models presented in Appendix D.

4.3. Properties derived from fitting images ℐ_ν

This section presents results for both starless cores and protostellar envelopes, obtained from successful fits of the background-subtracted ℐ_ν for all subsets Φ_n, for only the thinbody fitting model. Derived modbody masses are practically the same as the thinbody masses, because the bulk of the model mass is in optically thin regions (Sect. 3.3). Effects of the adopted far-infrared opacity slopes are the same as when fitting F_ν (Sect. 4.2): under- or overestimating β by a factor of 1.2 gives masses M_ℐ that are systematically under- or overestimated by a factor of 2. The method of fitting images ℐ_ν, thereby deriving , and afterwards integrating source mass M_ℐ brings clear benefits for well-resolved starless cores with nonuniform temperatures, compared to the other method (Sect. 4.2) of first integrating total fluxes F_ν from ℐ_ν (losing all spatial information) and then estimating M_F from the fitting model.

Isolated starless cores, β = 2 (Fig. 8). For the fully resolved models, derived T_ℐ and M_ℐ have fairly good accuracy and little bias for acceptable fits, although the range of the latter for larger Φ_n (n = 3 → 6) shrinks to the lowest masses. As the transparent low-mass cores (M → 0.03 M_⊙) are almost isothermal, derived T_ℐ and M_ℐ perfectly agree with T_M and M for any subset Φ_n. Massive cores with more variable T_d(r) (Fig. 3) also have significant variations of T_d(z) along the line of sight at pixel (i,j). Emission of hot dust skews the spectral shapes of I_νij towards shorter wavelengths, even more so at the high-mass end. The most accurate masses are obtained for Φ₂, whereas larger Φ_n (n = 3 → 6) give increasingly incorrect T_ℐ and M_ℐ. With degrading angular resolutions, the inaccuracies and biases increase, especially for M → 30 M_⊙ and larger Φ_n (n = 3 → 6). As expected, in the limiting case of unresolved objects the results approach those obtained with the method of fitting fluxes F_ν (Fig. 6).

Embedded starless cores, β = 2 (Fig. 8). For the fully resolved models, derived T_ℐ have fairly good accuracy and little bias for the acceptable fits, although the range of the latter for larger Φ_n (n = 3 → 6) shrinks to even lower masses than for the isolated models. Showing no particularly large bias over almost the entire range of model masses, M_ℐ are underestimated by a factor of 1.3 owing to the standard observational procedure of background subtraction (Appendix B, Sect. 5.2). Derived parameters of the models do not depend on angular resolutions, as they have relatively flat T_d(r) across their boundaries (Fig. 3), hence the spectral distortions of I_νij are negligible.

Isolated protostellar envelopes, β = 2 (Fig. 9). For the fully resolved models, derived T_ℐ and M_ℐ are very accurate across all masses and luminosities. With degrading angular resolutions and for larger Φ_n (n = 3 → 6) the inaccuracies and biases increase quite considerably. The accretion energy released in the envelopes centers heats the dust to T_S ~ 10³ K, making T_d(r) strongly nonuniform. For the lines of sight passing through the inner radial zones, the hot emission skews the I_νij shapes towards shorter wavelengths. For the unresolved envelopes, the results become similar to those obtained with the method of fitting fluxes F_ν (Fig. 7).

Embedded protostellar envelopes, β = 2 (Fig. 9). For the fully-resolved models, derived T_ℐ are slightly less accurate for the acceptable fits than in the case of the isolated envelopes. The range of the latter in more massive models for larger Φ_n (n = 3 → 6) shrinks towards higher L_⋆. The most accurate M_ℐ, obtained for Φ₂, are underestimated by a factor of 1.45, mostly because of the over-subtraction of the rim-brightened background (Appendix B, Sect. 5.2). For the partially resolved and unresolved envelopes, the most accurate M_ℐ (for Φ₂) are underestimated by factors of 1.5−2, whereas fitting larger Φ_n (n = 3 → 6) leads to errors by factors of 4−5.

5. Discussion

Spectral flux and intensity distributions of the radiative transfer models of the starless cores and protostellar envelopes (0.03−30 M_⊙,L_⊙) were fitted using the modbody and thinbody models. Derived values of the fitting parameters were then compared to their true values to quantify the qualities of the mass derivation methods, fitting models, and various sources of errors.

As shown in Sect. 4, large intrinsic inaccuracies and biases need to be taken into account when applying the methods of mass derivation to the observed sources. In addition to being affected by nonuniform temperatures, estimated masses are also affected by the adopted value of β and subset of data points Φ_n, as well as by the removal algorithm of the background emission of an embedding cloud. In the method of fitting fluxes F_ν, the masses depend on the fitting model, whereas in the method of fitting images ℐ_ν, they depend on the angular resolution.

The results of this purely model-based work discussed below may be directly applicable only to sources with very accurate measurements (with negligible errors). Real observations deal with images of relatively faint, crowded sources on strong and variable backgrounds, obtained with quite different angular resolutions, and thus they carry much larger measurement errors. Observations are substantially affected by various statistical and systematic errors, depending on the adopted source extraction method (e.g., Men’shchikov et al. 2012; Men’shchikov 2013) and especially the treatment of background subtraction and deblending. Implications for the real-life studies are considered below, whenever possible.

5.1. Mass derivation methods

In the first method, source fluxes F_ν are integrated from the images ℐ_ν, their spectral distribution is fitted, and source mass M_F is estimated from the fitting model. In the second method, the pixel spectral shapes I_ν of the images ℐ_ν are fitted and the source mass M_ℐ is integrated from the resulting image of column densities. For unresolved sources and the thinbody fitting model, the methods give very similar levels of inaccuracy, whereas for resolved images, the methods differ quite substantially.

When fitting F_ν, the observed source emission from its entire volume is blended in the spatially integrated fluxes that retain no spatial information. For the models with strongly nonuniform T_d(r) (Fig. 3), resulting heavy distortions of the spectral shapes of F_ν (Fig. 4) from those of the fitting models lead to large systematic errors in estimated parameters (Figs. 6 and 7).

When fitting ℐ_ν, it is very beneficial to have a higher angular resolution. For fully resolved objects, pixels (i,j) sample independent I_νij from different columns of dust. For the transparent lower mass models (M ≲ 0.3 M_⊙), derived M_ℐ are quite accurate (Figs. 8, 9). For lower resolutions, the intensity of each pixel (i,j) blends with that of its larger surroundings within the beam, not only along the line of sight. The contamination of I_νij by the more distant areas, leads to a substantial degradation of T_ℐ and M_ℐ, especially when fitting large Φ_n (n = 3 → 6). Thus, the benefits of this method are vanishing with decreasing angular resolutions.

Multiwavelength Herschel images have been used to reconstruct radial temperature and density profiles of well-resolved sources (Roy et al. 2014). Whenever such reconstructed densities are accurate enough, they can be used to obtain masses of the nearby sources. Results of this study demonstrate, however, that the simple method of fitting images ℐ_ν is able to deliver accurate masses for spatially resolved sources (Sect. 4.3).

5.2. Background subtraction

Stars form in the densest parts of interstellar clouds, hence the embedded models of starless cores and protostellar envelopes must be more realistic than the isolated models. Although the spherical uniform-density embedding clouds are idealized, in a first approximation they account for the absorption and re-emission of ISRF, leading to realistic temperature profiles within the model objects. However, the presence of surrounding material makes it necessary to subtract its contribution to study the properties of the starless cores and protostellar envelopes alone. In observational practice, backgrounds are estimated by an average intensity in a narrow annulus placed just outside a source (cf. Sect. 2.5). Subtraction of such a flat background is not quite accurate as a transparent embedding cloud around any object always tends to be rim-brightened and resembles a crater, in contrast to a distant, physically unrelated back- or foreground. This effect is discussed in detail in Appendix B.

The actual observable depths of the background craters may be shallower, when the local (filamentary) background itself is embedded in a less dense but more extended cloud or is seen in projection onto a distant, physically unrelated back- or foreground. The rim-brightening effect gets diluted, if the column densities of the source-embedding background and of the other unrelated clouds are similar. Poorer angular resolutions also tend to smear out the effect for less resolved sources. Realistic temperature gradients within the embedding backgrounds (Fig. 3) can either reduce or increase the crater depths by ~ 10% for starless cores and protostellar envelopes, respectively (Appendix B).

For unresolved sources, the observational algorithm of background subtraction is likely to overestimate fluxes as stars are born within the gravitationally unstable densest peaks of the parent clouds. Large beams blend the object’s emission with that of its mountain-like environment, spreading the mix downhill, towards the valleys of lower cloud densities. The real background under an unresolved source must be hill-like, whereas the background values from an annulus tend to come from more distant valleys. The problem is aggravated in crowded regions, where no local source-free annuli around overlapping sources can be found and where one needs to deblend sources. Angular resolution degrades with wavelengths, hence the degree of flux overestimation becomes strongly biased towards longer wavelengths.

5.3. Nonuniform temperatures

Both fitting models make a sensitive assumption that the objects have a uniform temperature T, which seems to make them inadequate for the applications to starless cores and protostellar envelopes with nonuniform T_d(r). For the purpose of the derivation of accurate masses, however, the uniform T can be interpreted as an appropriate average quantity. In the methods of fitting F_ν and ℐ_ν, the temperature is consistent with T_M and T_Nij from Eqs. (4) and (11), respectively (cf. Sects. 3.1 and 3.2). In other words, to estimate masses M_F or M_ℐ that are accurate (≈M), it is necessary that the fitting returns T_F or T_ℐij as close as possible to the average values T_M or T_Nij, respectively. This is clearly demonstrated in Appendix D by the accurate masses obtained for the isothermal models with T_d(r) = T_M.

The inhomogeneous temperatures tend to distort the spectral shapes of F_ν and I_νij of the objects towards shorter wavelengths (Figs. 4 and 5). With a strong dependence of the dust emission peak on temperature (κ_νPB_νP(T) ∝ T⁵), the radial zones with higher T_d(r) make a much greater contribution to the observed spectral shapes. Therefore, the shapes are skewed mainly owing to the emission of those parts of the objects that have T_d(r) >T_M or T_d(r) >T_Nij. In other words, distortions of the spectral shapes are caused by the dust with excess temperatures above the average values.

This is further demonstrated by additional ray-tracing observations of the models, in which the excess temperatures were removed: T_d(r) → min{T_d(r),T_M}. Derived masses of these mostly isothermal models (not shown) are almost as accurate as those of the fully isothermal models (Appendix D), only within a few percent lower. As is expected, there is almost no dependence on the subsets Φ_n (n = 2 → 6), which indicates that the spectral shapes are indeed not distorted.

5.4. Fitting models

When fitting images ℐ_ν, both fitting models are equivalent and estimated parameters are indistinguishable (Sect. 3.3). When fitting fluxes F_ν, the results of this work show that thinbody generally returns far more accurate masses than modbody does for both isolated and embedded variants of starless cores and protostellar envelopes (Figs. 6 and 7).

Although the modbody fits often look better (i.e., they have smaller χ² values), they generally bring parameters that are much more inaccurate. Indeed, the spectral shapes of F_ν are skewed towards short wavelengths by emission from their hotter parts. With more free parameters, modbody describes more flexible shapes, between B_ν(T) and κ_νB_ν(T). It is able to produce better fits of the distorted spectral shapes of objects with nonuniform T_d(r) and hence it always tends to produce significantly over- and underestimated T_F and M_F, respectively (Sect. 4.2). Furthermore, most of the modbody fits have τ_ν ~ 1 even in the far-infrared, which is fundamentally inconsistent with the radiative transfer models whose fluxes F_ν represent optically thin emission (τ_ν ≪ 1).

The thinbody model produces the best overall results and smallest biases and inaccuracies in derived T_F and M_F for the isolated and embedded starless cores and protostellar envelopes (Figs. 6−9). The thinbody fits are, by definition, optically thin in the far-infrared and thus consistent with the radiative transfer models. Only two variable fitting parameters of thinbody contribute to better robustness of T_F and M_F, compared to modbody with one extra free parameter.

Contrary to what is usually assumed in observational studies, the results show that it must be counterproductive to aim at precise fitting of the peaks and shorter wavelength parts of I_ν and F_ν. When the distorted shapes are reproduced more accurately, the estimates of the temperatures and masses are less accurate.

5.5. Opacity slopes

The standard methods of mass derivation ignore the presence of very small stochastically heated dust particles, assuming just a simple power-law opacity across all bands, and so do the radiative transfer models in this study. Emission of such very small grains within the real objects could enhance fluxes at 70 and 100 μm and, in effect, skew their spectral shapes farther towards short wavelengths, leading to more heavily overestimated temperatures and underestimated masses.

Various compositional and structural properties of real cosmic dust grains in different environments may lead to far-infrared opacity slopes that are different from β ≈ 2 (expected for small compact spherical grains) and even to wavelength-dependent β_λ. This study explored three constant values (1.67, 2.0, 2.4) to probe their influence on the accuracy of derived masses. Fixing β in the fitting process reduces the number of free parameters and improves the consistency (reduces biases) of derived parameters for objects with different physical properties (M, L_⋆).

With the correct β value, masses derived with the thinbody fits are generally off the true mass M, the magnitude of discrepancy depending on how much and in what direction derived temperature deviates from T_M. When β is over- or underestimated by a factor of 1.2, derived masses become over- or underestimated within a factor of 2 with respect to the masses obtained using the true value β = 2. This is a direct consequence of the temperatures being under- or overestimated, correspondingly, a behavior that is easy to understand. In contrast to the thinbody fits, no clear trends with respect to the inaccuracies in the adopted β value can be found for modbody, except that it generally returns greatly over- and underestimated T_F and M_F.

To quantify the effects of freedom in this fitting parameter, additional fits with variable β were performed (Appendix E). As is expected, they showed much greater biases and inaccuracies in derived parameters (Figs. E.1 and E.3), as the extra degree of freedom also makes the resulting β values incorrect (Fig. E.2), the magnitude of error depending on the true values of M and L_⋆. It is possible to compare these results with those obtained in previous studies focused on the relationships between the derived β_F and T_F (Shetty et al. 2009a,b; Juvela & Ysard 2012, and references therein). The present models of starless cores and protostellar envelopes show that the correlations of the two quantities may be both positive and negative (Fig. E.4), with almost no correlation in the case of isothermal models. They must be induced by deviations of the spectral shapes of F_ν (Fig. 4) from F_ν(T_M) (Fig. D.1), caused by the nonuniform T_d(r) (Fig. 3). For the protostellar envelopes, the correlations are non-monotonic and they may either be strongly negative or positive, depending on the luminosity.

5.6. Data subsets

For nonuniform profiles T_d(r) of starless cores and protostellar envelopes (Fig. 3), better parameters are estimated with thinbody when using smaller subsets of data (n = 6 → 2) as the latter are less affected by the skewed spectral shapes. The most accurate masses are obtained by fitting just two of the longest wavelength data points; in most cases, however, a subset Φ₃ produces very similar results. Larger subsets Φ_n (n = 2 → 6) may give slightly better M_F only when fixing an incorrect β value for the lower-mass starless cores (M ≲ 1 M_⊙, Fig. 6). Using the inadequate fitting model with an incorrect β, larger Φ_n can constrain T_F to better resemble T_M. For overestimated β, derived M_F always shift to higher values (Fig. 7), which offsets the general opposite trend to underestimate M_F and thus may give more accurate results.

Inaccuracies of the data points in real observations are usually more substantial than those assumed in this work, aggravated by the systematic uncertainties that may lead to both over- and underestimated F_ν (Sect. 5.2). Background-subtracted and deblended I_ν at each wavelength with different angular resolutions have independent and different systematic errors. The latter must be large and uncertain on the bright and structured backgrounds in star-forming regions. Moreover, unresolved sources are likely to include emission from clusters of objects.

Results of this model study are directly relevant to real observations only in the simplest case (which is rare) of accurate measurements with negligible errors. A blind application of the findings to real complex images may lead to incorrect results if the above caution is ignored and small subsets Φ₂ of data points with large and independent measurement errors are fitted. For such data, it would be safer and more appropriate to fit a larger subset of the longest wavelength data (Φ₃ or Φ₄, depending on the quality of measurements). Distortions of the observed spectral shapes towards shorter wavelengths is an intrinsic property of both starless cores and protostellar envelopes, affecting all sources in star-forming regions, independently of the level of measurement errors.

It beyond the scope of this model-based work to give general recipes to observers on how to select data points to fit. This study highlights the intrinsic behavior of the mass derivation methods by eliminating the “observational layer” (with all its complications and uncertainties) between the objects and the observer. It is important to realize that the peak and shorter wavelength shapes of I_ν and F_ν are most skewed by the temperature excesses within objects (Sect. 5.3) and their influence has to be minimized to obtain accurate results. In view of the strong dependence of the results on Φ_n, it is advisable to examine fits of all subsets of data points for each observed source to estimate the robustness of the results and to possibly choose the fits giving the best mass estimate.

5.7. Mass uncertainties

To make a bridge between this purely model-based study with no measurement errors and actual observational studies and see how typical statistical errors in the input data would translate into those of the derived masses, this work assigned (fairly optimistically) ±15% errors to the model intensities and fluxes, and adopted ±20% errors in η, κ₀, and D.

The uncertainties in derived masses returned by the fitting algorithm are 40−70%, depending on the subset Φ_n of fluxes (Figs. 6−9, D.2, D.3). For the acceptable fits of larger subsets (n = 2 → 6), the derived mass uncertainty is dominated by the 20% errors of the parameters η, κ₀, and D, because the effect of the 15% measurement errors becomes smaller for the fits constrained by a larger number of independent data points. For smaller subsets (n = 6 → 2), the fits are less constrained, hence the contribution of the 15% error bars to the derived mass uncertainty becomes larger. Different subsets Φ_n give very similar results only for fully resolved sources with the method of fitting images ℐ_ν (Figs. 8 and 9).

In real observations, statistical measurement uncertainties in I_ν and F_ν are larger than the ±15% errors assumed in this study. Furthermore, it would be more realistic to adopt uncertainties of η, κ₀, and D of at least ±50%, which would raise the derived mass uncertainties well beyond 100%. By including the mass inaccuracies (of a factor of 2) induced by a 20% uncertainty in β and systematic errors (of factors of at least 2) caused by the nonuniform temperatures within the observed sources, it is clear that the absolute values of masses derived from fitting are inaccurate and uncertain (within a factor of at least 2−3). It is possible to neglect the uncertainties in η, κ₀, and D, if the focus is on studying relative properties of a population of objects all at roughly the same distance within a certain star-forming cloud with homogeneous dust properties. Apart from this, however, one has to derive accurate absolute values of the most fundamental parameters to make correct and physically meaningful conclusions.

It is quite important to carefully estimate mass uncertainties: without realistic error bars, derived masses are meaningless and correct conclusions are unlikely. To go one step further and obtain an idea of the actual errors of derived masses, it is possible to construct radiative transfer models of the observed population of sources, distribute the model sources over the observed images and extract them, and finally derive their masses. Comparing derived masses with the fully known model properties, reasonable estimates of the actual errors in derived masses are obtained.

6. Conclusions

This paper presented a model-based study of the uncertainties and biases of the standard methods of mass derivation (fitting fluxes F_ν and images ℐ_ν), widely applied in observational studies of the low- and intermediate star formation. To focus on the intrinsic effects related to the physical objects, all observational complications leading to additional flux or intensity errors (filamentary and fluctuating backgrounds, instrumental noise, calibration errors, different resolutions, blending with nearby sources, etc.) were assumed to be nonexistent. As a consequence, results of this work are directly relevant only for the simplest case of bright isolated sources on faint backgrounds with negligible measurement errors. The real mass uncertainties for starless cores and protostellar envelopes are likely to be larger than those found in this work.

Embedding backgrounds of physical objects tend to be rim-brightened (i.e., resembling craters), their depths depend on the sizes of the object and embedding cloud. The standard observational procedure of flat background subtraction may give systematically underestimated I_ν and F_ν, and hence masses for resolved sources. Poorer angular resolutions at longer wavelengths tend to systematically overestimate I_ν and F_ν, and hence masses for unresolved objects, as their emission gets blended with that of the mountain-like background and possibly with other objects within the same beam.

Temperature excesses above average values T_M and T_Nij are the primary reason for the skewness of the spectral shapes of F_ν and I_νij towards shorter wavelengths. Depending on M, L_⋆, β, Φ_n, fitting model, and angular resolution, they lead to overestimated temperatures and various biases. With the method of fitting F_ν, masses become underestimated by factors 2−5. When fitting ℐ_ν, similarly large inaccuracies are found only for unresolved objects, whereas with better angular resolutions they decrease and become very small for well-resolved objects.

When fitting ℐ_ν, both fitting models are equivalent and estimated M_ℐ are indistinguishable. When fitting F_ν, thinbody gives far more accurate M_F than modbody does. The latter causes such great biases and inaccuracies in T_F and M_F that modbody must be considered unusable.

Fixing opacity slope β reduces biases in derived parameters. When β is too high or low by a factor of 1.2, derived masses become over- or underestimated by a factor of 2 with respect to those obtained using the true β = 2. Qualitatively, this behavior is caused by the natural tendency of steeper β to produce lower temperatures, hence higher masses. Quantitatively, the factors are approximate and they may depend on some of the assumptions used in this study. Mass derivation with a free variable β should be avoided, as it tends to lead to very strong biases and erroneous masses.

Subsets Φ_n of data points strongly affect derived masses, except when fitting images ℐ_ν of fully resolved sources. Given the nonuniform T_d(r) of the model objects, the most accurate masses are estimated with thinbody using subsets that are as small as possible (n = 6 → 2). In real observations with substantial independent errors in different wavebands, it should be much safer and more accurate to fit slightly larger subsets (Φ₃ or even Φ₄). Those data points that are on the peak of their spectral distribution or on the short-wavelength side should be ignored, whenever possible, to improve the accuracy of derived masses. In practice, it is advisable to investigate fits of all subsets of data for each observed source, to verify robustness of the results and to possibly choose the best mass estimate.

Inaccuracies of the derived masses depend on the mass range considered. Dividing the range of 0.03−30 M_⊙ at 1 M_⊙ into the low- and high-mass objects and considering unresolved or poorly resolved sources with β = 2, the following conclusions drawn. Masses of the isolated low- and high-mass starless cores are underestimated by factors 1−1.3 and 1.3−4, respectively. The mass inaccuracies increase towards the high-mass end and for larger subsets Φ_n (n = 2 → 6). Masses of the embedded low-mass cores are underestimated by a factor of 1.4. They are more biased towards the high-mass end, changing from under- to overestimated within a similar factor. Masses of the protostellar envelopes are considerably biased over the range of 0.03−30 L_⊙ and their inaccuracies strongly increase for larger subsets of Φ_n (n = 2 → 6). Masses of the isolated and embedded envelopes become underestimated by factors 2−3 and 3−5, respectively. Masses of the low-mass starless cores are likely to be determined much more accurately than those of protostellar envelopes.

Typical mass uncertainties returned by the fitting algorithm are 40−70%, adopting statistical errors of 15% for model intensities (fluxes) and optimistically assuming that η, κ₀, and D were known to within 20%. If more realistic statistical errors in the measurements and parameters of at least 50% are adopted, the mass uncertainties increase well beyond 100%. Larger subsets Φ_n (n = 2 → 6) of independent data points are beneficial in somewhat reducing the resulting mass uncertainties. On the other hand, the larger subsets are also highly undesirable, because they escalate the systematic mass inaccuracies by at least a factor of 2 as a result of nonuniform temperatures. Smaller subsets Φ_n (n = 6 → 2) are able to minimize the systematic errors caused by the temperature variations, but they increase the chances of getting incorrect masses in the case of inaccurate data measurements in real observations.

Without extremely accurate flux measurements and knowledge of the free parameters (η,κ₀,β,D), and without radiative transfer simulations to have an idea of the actual mass errors, it would be reasonable to assume that the absolute values of masses of the unresolved or poorly resolved objects are inaccurate to within at least a factor of 2−3. This may be less problematic, if the relative properties are studied of a population of objects within a star-forming cloud, hopefully with the same distance and dust opacities. Ultimately, however, accurate absolute masses are necessary to make correct, physically meaningful conclusions.

There are several ways to improve mass estimates: (1) using a multiwavelength source extraction method measuring the most accurate, least biased background-subtracted and deblended fluxes across all wavebands; (2) selecting the best sources from the extraction catalogs, with the most accurately and consistently measured fluxes over at least three longest wavelengths; (3) using the thinbody fitting model for the purposes of temperature or mass derivation; (4) estimating the model parameters η, κ₀, β, and D as accurately as possible and always performing fitting with β fixed; (5) for resolved sources, fitting their background-subtracted (and deblended) images and integrating source masses from column densities; (6) fitting all subsets of data points for each source and choosing the smallest possible subset (Φ₃ or Φ₄) that gives the most accurate temperatures and masses; (7) using radiative transfer models to simulate observed images, extracting the model sources, deriving their masses, and comparing them to the true model values to have an idea of the actual errors for the derived masses of observed sources.

---

Acknowledgments

This study employed SAOImage DS9 (by William Joye) developed at the Smithsonian Astrophysical Observatory (USA), the CFITSIO library (by William D. Pence) developed at HEASARC NASA (USA), and SWarp (by Emmanuel Bertin) developed at Institut d’Astrophysique de Paris (France). Radiative transfer code MC3D-sph version 3.12 (by Sebastian Wolf, Wolf 2003) was used to compute the first generation of the models in this work. The plot utility and ps12d library used in this work to draw figures directly in the PostScript language were written by the author using the PSPLOT library (by Kevin E. Kohler) developed at Nova Southeastern University Oceanographic Center (USA) and the plotting subroutines from the AZEuS MHD code (by David A. Clarke and the author) developed at Saint Mary’s University (Canada). Collaborative work within the Herschel Gould Belt and HOBYS key projects was very beneficial. Useful comments on a draft made by Pierre Didelon, Arabindo Roy, Alana Rivera-Ingraham, Sarah Sadavoy, Philippe André, and by the anonymous referee helped improve this paper.

References

Appendix A: Protostellar envelopes temperatures

With the adopted κ_ν and ρ(r) (Sects. 2.1 and 2.2), the radial temperature profiles of protostellar envelopes (Fig. 3) can be approximated by a combination of two power laws: (A.1)where r is in AU and parameters A and B depend on mass and accretion luminosity: Equation (A.1) describes the temperatures induced by the central accretion energy source (ignoring ISRF), valid for T_d ≲ 300 K. The first term in Eq. (A.1) approximates the steepest profiles in the inner semi-opaque region, the second term represents temperatures in the transparent outer part of the envelopes. An approximate borderline between the two regimes can be estimated directly from Fig. 3 as (A.2)

Appendix B: Rim-brightened backgrounds

In contrast to the emission of the distant and physically unrelated backgrounds or foregrounds, embedding backgrounds resemble craters. The central spherical region occupied by an object (r ≤ R) does not belong to the embedding cloud (R<r ≤ R_E) and thus must be considered empty when determining the object’s background. Rim brightening for uniform-density transparent isothermal clouds with such a cavity depends only on their relative radial thickness (R_E−R) /R. It can be quantified by the ratio f_S of intensities (or column densities) along the lines of sight passing through the rim and the center of the cavity: (B.1)According to Eq. (B.1), the background under embedded objects can be overestimated from a few percent to a factor of several, hence background-subtracted values and masses may become substantially underestimated. For the present models with R = 10⁴AU and R_E/R = 3, the background and masses are over- and underestimated by f_S = 1.41, respectively. The value is the discrepancy of derived masses seen for the embedded starless cores and protostellar envelopes (Figs. 6−9, D.2, D.3). The effect becomes much stronger for very thin shell-like embedding clouds, whereas it vanishes for extended background clouds. For the size ratios R_E/R of 1.3 and 30, the factor f_S takes the values of 2.77 and 1.03, respectively. Numerical examples of the rim-brightened backgrounds for both isothermal and non-isothermal spherical embedding clouds are shown in Fig. B.1.

Realistic temperature profiles in the embedding clouds bring only minor quantitative changes, not altering qualitatively the rim brightening effect (Fig. B.1). Steep positive or negative temperature gradients in the dense shells around embedded starless cores and protostellar envelopes (Fig. 3) tend to slightly reduce or increase the brightening effect (by ~ 10%), respectively.

Actual geometry of the real background clouds is of minor importance, the only relevant assumption is that the embedded object (hence, its background cavity) has a convex shape. For instance, assuming a plane-parallel embedding cloud with thickness 2 R_E along the line of sight, it is possible to obtain a slightly different expression than Eq. (B.1) for the rim brightening factor: (B.2)Plane-parallel geometry makes the brightening factors f_P somewhat larger than f_S from Eq. (B.1), with the difference being stronger for thinner shell-like clouds. For instance, the size ratios R_E/R of 30,3, and 1.3, correspond to the brightening factors f_P of 1.03,1.46, and 3.55, respectively.

Depths of the background craters may be quite dissimilar for different objects in real observations. The observations show that the interstellar medium is strongly filamentary and that stars tend to form in narrow, very dense filaments (e.g., Men’shchikov et al. 2010; André et al. 2014). For an object embedded in a cylindrical filament of radius R_E in the plane of the sky, the brightening factor f_C is intermediate between f_S and f_P, depending on the position angle of the radius-vector from the center of the object to its outer boundary R. It is easy to see that f_C = f_P along the filament’s axis, whereas f_C = f_S in the orthogonal direction, across the filament.

The widths of the embedding filaments appear to be similar to the sizes of embedded objects (Men’shchikov et al. 2010). Assuming their cylindrical geometry, the embedding filaments of starless cores and protostellar envelopes are likely to have R_E only a factor of about 2−3 larger than R. For such narrow filamentary backgrounds of resolved objects, the rim brightening effect is quantified by factors f_C ≈ 1.8−1.4.

Background rim brightening may be observable only when imaging the nearby resolved sources unaffected by other distant back- or foregrounds. With this effect in action, the standard observational algorithm of background subtraction underestimates F_ν by factors similar to f_S or f_P. With poorer angular resolutions, the rim of the background crater gets smeared out and thus the brightening effect eventually vanishes for unresolved sources which also have an opposite trend to produce underestimated background and hence overestimated I_ν and F_ν (cf. Sect. 5.2).

Appendix C: Fitting procedure

The nonlinear least-squares fitting algorithm used in this work employs the Levenberg-Marquardt method (Press et al. 1992) that minimizes χ² residuals between the model and data. The method requires a user to provide initial guesses for model parameters. Tests have shown that an arbitrary choice of the initial values of the fitting models (Sects. 3.1, 3.2) does not guarantee convergence to the global χ² minimum. A fully automated fitting procedure has been designed to overcome this problem and avoid any need to make arbitrary initial guesses.

The algorithm explores the multidimensional parameter space of the model with a large number of trial fittings of the spectral distributions of data points. The parameter space is discretized in logarithmically equidistant steps δlog p, covering all relevant initial values of temperature T (4−10³ K), mass M (3 × 10^-3−3 × 10²M_⊙), column density N_H2 (10¹⁷−10²⁷ cm^-2), and solid angle Ω (1.85 × 10^-13−1.85 × 10^-5 sr). Large initial discretization steps δlog p ~ 10 and the above ranges of parameters are adaptively refined in an iterative binary search down to δlog p ≈ 0.1, accelerating the algorithm in finding the global minimum. Data points are fitted using all combinations of the model parameters⁶ in the adaptively refined parameter space and their initial values are found that converge to the globally smallest χ² in a fully automated procedure.

Algorithms described in this paper were written as a versatile and robust FORTRAN utility fitfluxes that efficiently estimates modbody or thinbody parameters from fitting either total fluxes of cataloged sources or intensities of multiwavelength images. The code is easy to install and use and it is freely available from the author upon request.

Appendix D: Results for isothermal models

The importance of estimating T_F that approaches the mass-averaged temperature T_M from Eq. (4) for deriving accurate masses is shown by isothermal models, those described in Sect. 2 and used throughout this paper where self-consistent (radiative-equilibrium) profiles T_d(r) were replaced with T_M. The isothermal models were then observed and imaged in a ray-tracing run of the radiative transfer code. The resulting SEDs (Fig. D.1) are essentially the modified blackbody shapes κ_νB_ν(T_M) ν and the same is true for the spectral shapes of image pixels. The temperature excesses above T_M (Sect. 5.3) that greatly distorted the model SEDs (Fig. 4) towards shorter wavelengths do not exist in the isothermal models. Consequently, the isothermal shapes of F_ν and I_νij bring much more expected and accurate results.

With the correct value β = 2, derived M_F of the isolated starless cores and protostellar envelopes derived with thinbody agree with the true masses M (Fig. D.2). The same results are obtained for modbody, with the exception of the minimal subset Φ₃ for some models. An inspection of the problematic fits show that T_F is somewhat overestimated because of an additional degree of freedom in modbody and a formal search for a globally best fit in its parameter space. The fits in question have the globally lowest χ² value, but they correspond to small values of Ω and, therefore, to high optical depth τ_ν ~ 1. However, there are also other very good fits with somewhat larger χ² that do produce accurate T_F = T_M with τ_ν ≪ 1, consistent with the models. The problem seems to be just a simple consequence of the finite accuracy of the numerical models and their fluxes.

With the same β = 2, the derived masses of the embedded models (Fig. D.2) are almost uniformly underestimated by a factor of approximately 1.3. The reason for the difference with respect to the isolated models is the conventional approach to background subtraction. Emission of a transparent cloud embedding a physical object tends to be rim brightened (Appendix B, Sect. 5.2). An average intensity in an annulus may overestimate background, from a few percent to a factor of several, hence may underestimate the background-subtracted intensities I_ν, fluxes F_ν, and derived masses (see Sects. 4.2 and 4.3).

For an inadequate fitting model, skewed by β values fixed below or above its correct value, the fits are obviously biased to over- or underestimate T_F and hence to under- or overestimate M_F, correspondingly (Fig. D.2). The inaccuracy is within a factor of 2, independently of whether β is a factor of 1.2 lower or higher. Fitting larger subsets Φ_n (n = 2 → 6) for the isothermal models may give somewhat better derived parameters, as they better constrain T_F with an incorrect slope β.

Mass derivation from images ℐ_ν of the isothermal models delivers results that are similar to those described above for both isolated and embedded variants (Fig. D.3). Dependence on the adopted β is the same as described above, hence only the results with correct β = 2 are presented. With no temperature deviations from T_M in the isothermal models, derived M_ℐ are very accurate for all angular resolutions, in contrast to the results with the self-consistent T_d(r) (Figs. 8 and 9). Derived M_ℐ for the embedded models are practically identical to M_F, being underestimated by a factor of 1.3 owing to the background rim-brightening effect (Appendix B, Sect. 5.2).

Appendix E: Results for free variable β

In some applications of the mass derivation methods, the opacity slope β has been allowed to vary along with the other fitting parameters (T, M, Ω). To quantify effects of the extra degree of freedom, additional fits with variable β were performed in this study. Although the parameters for both fitting models were derived and analyzed, only the much less incorrect thinbody results are presented and discussed here.

Figure E.1 compares derived T_F and M_F of the isolated, embedded, and isothermal variants of starless cores and protostellar envelopes with the true model values (T_M, M). Although the isolated cores and envelopes display behavior that is qualitatively similar to the β = 2 case (Sect. 4.2), the biases in T_F and M_F towards denser (more massive) cores and envelopes, as well as over L_⋆ for the latter, become much stronger. For example, for the isolated starless cores with M ≳ 10 M_⊙, derived T_F are overestimated by a factor of 2, whereas M_F are underestimated by a factor of 15. For the embedded protostellar envelopes of M = 3 M_⊙ with L_⋆ ≲ 1 L_⊙, temperatures are overestimated by a similarly large factor and M_F underestimated by a factor of 8.

Such errors are caused by the derived β_F whose values for starless cores are systematically lowered towards higher mass models, and are underestimated by a factor of 4 (Fig. E.2). In the case of the embedded protostellar envelopes, the values of β_F are progressively underestimated towards lower luminosities, up to a factor of 2 (Fig. E.2). The very large errors in β_F are, in turn, caused by the temperature excesses over T_M (Sect. 5.3), which is highlighted in Fig. E.1 by the accurate results for the isothermal models. As in the fixed β case, errors in derived parameters for the embedded starless cores are smaller than those for the isolated cores, whereas the behavior is opposite for the embedded protostellar envelopes (cf. Figs. 6, 7). However, the biases over the masses and luminosities in Fig. E.1 become stronger than in the fixed β case, which again is attributed to the additional biases in the derived β_F values (Fig. E.2).

The method of fitting images ℐ_ν delivers results that are similar to those described above, for both isolated and embedded variants (Fig. E.3) of partially resolved and unresolved objects. Derived masses for the fully resolved objects are much more accurate and they do not depend on the subset of data points Φ_n, for the reasons discussed in Sects. 4.3 and 5.1. With degrading angular resolutions, accuracy of the estimated parameters deteriorates to the levels obtained from the method of fitting total fluxes F_ν (Fig. E.1). Larger beams heavily blend emission with nonuniform temperatures from different pixels, distorting their spectral distribution towards shorter wavelengths. In both methods, derived masses become systematically much less accurate (greatly underestimated) when fitting larger subsets Φ_n (n = 3 → 6). However, even the smallest subsets Φ₃ show very significant inaccuracies and different biases that depend on the mass and luminosity of an object.

The above results obtained with free fitting parameter β can be compared with those from the previous studies focused on the relationship between derived temperature and opacity slope (cf. Shetty et al. 2009a,b; Juvela & Ysard 2012, and references therein). Figure E.4 shows the intrinsic dependencies between β_F and T_F for the isolated and embedded starless cores and protostellar envelopes, and for the isothermal versions of the isolated models. The correlations between β_F and T_F take variety of shapes, from a strongly positive to a strongly negative correlation, with practically no correlation for the isothermal models. This suggests that they must be caused by different kinds of deviations of the spectral shapes of F_ν (Fig. 4) produced by the nonuniform T_d(r) (Fig. 3) from F_ν(T_M) (Fig. D.1). Smaller subsets Φ_n (n = 6 → 3) bring less correlated β_F and T_F than the large subsets do. For the protostellar envelopes, the correlations are non-monotonic and they may either be strongly negative or positive, depending on the luminosity of the central energy source.

All Figures