© The Authors 2023

1 Introduction

Thermal dust emission at far-infrared and millimetre wavelengths is used to trace the structure and physical conditions of interstellar matter (ISM), especially in connection with star formation (SF). The targets range from a Galactic diffuse medium and local molecular clouds to entire galaxies. The measured intensities represent emission from a range of conditions that depends on both the nature of the sources and the angular resolution of the observations. The continuum emission also carries information about the dust itself, which is interesting due to the significant and systematic evolution that takes place in the dust properties during the SF process (Draine 2003).

The long-wavelength dust emission is often modelled as modifed black-body (MBB) emission, a black-body spectrum multiplied by a power law representing the frequency dependence of the dust absorption cross section κ_ν. The analysis may assume both a single temperature T and optically thin emission, which is only an approximation of the emission of any astronomical source. If the assumptions of a single temperature and optically thin emission are both violated, an analytical estimation of source properties is hardly possible, and more detailed modelling of the observations is needed. Even for optically thin sources, all observations consist of emission from dust at different temperatures within the telescope beam, both along the line of sight (LOS) and over the sky. The frequency dependence of κ_ν is also more complex than a single power law (Paradis et al. 2010; Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014; Gordon et al. 2014; Köhler et al. 2015; Juvela et al. 2015a). However, given the limitations as to the number of observed frequencies and the signal-to-noise ratio (S/N) of the data, the spectral energy distributions (SEDs) have to be modelled with just a few free parameters, which include at least the intensity and one temperature. The simultaneous fitting of temperature and β is complicated by the inherent anti-correlation of these parameters, which makes the individual parameters very sensitive to noise (Shetty et al. 2009a; Kelly et al. 2012; Juvela & Ysard 2012b; Juvela et al. 2013b). The problem is further complicated by temperature variations, which tend to systematically decrease the observed β values relative to the actual β of the dust (Shetty et al. 2009b; Juvela & Ysard 2012a). If dust emission is only used to measure column densities, the results are likely to be more robust if β is kept constant, even if the real variations of β may then result in up to tens of percent errors in the column density estimates. Spectral indices can only be constrained with high-S/N data that cover a wide frequency range, and even then the estimates can still be biased due to the temperature variations.

The single-component MBB model can be extended by allowing multiple temperature or β components. In studies of dust properties, one may further allow changes in β as a function of frequency (Paradis et al. 2012; Gordon et al. 2014; Galliano 2018). However, many sources are known to have wide temperature distributions, with smaller effects expected from the β variations. Therefore, it is more common to use multiple temperature components while the dust optical properties are assumed to be constant or modelled with one free parameter at most. This is the case especially in studies of external galaxies, where most temperature structures remain unresolved (Galametz et al. 2012; Chang et al. 2021). However, temperature variations can be significant even in local SF clouds, from ~6 K inside cold cores (Crapsi et al. 2007; Harju et al. 2008; Pagani et al. 2015) to more than 100 K in hot cores (Motte et al. 2018; Jørgensen et al. 2020). In Galactic studies, the single-component MBB model is still the main tool, although it suffers from the well-known tendency to underestimate column densities. The actual three-dimensional temperature fields can only be estimated with radiative transfer modelling or, in the case of geometrically simple objects, with tools such as the inverse Abel transform (Roy et al. 2014).

The calculation of multi-component MBB fits poses some technical challenges. The increased number of free parameters may require regularisation, for example by using parameter priors in the Bayesian framework. Especially the problem of Τ – β correlations has led to the use of hierarchical models, where the hyper-parameters of the parameter distributions are fitted together with the parameters (i.e. Τ and β) of individual measurements (Kelly et al. 2012; Galliano 2018; Lamperti et al. 2019). The efficiency of priors and hierarchical models still varies from case to case. In a complex SF region, the priors must accommodate a wide range of temperatures and provide correspondingly weaker constraints for the individual measurements. There is also the risk of biasing the estimates for rare or unexpected sources, such as an individual hot (or cold) core seen against a field of colder (or warmer) background (Juvela et al. 2013b).

Marsh et al. (2015) have presented the point process mapping (PPMAP) method, which fits dust observations with a combination of fixed SED templates, such as MBB functions (Marsh et al. 2017; Howard et al. 2019). The emission is represented by points in the state phase, and a requirement of using the minimum set of points to fit the data (to an accuracy of χ² ~ 1) also acts as regularisation. One important point is that PPMAP uses all observations at their native resolution. The more typical approach is to convolve all observations to a common lower resolution, which then allows the SED of each map pixel to be analysed independently. The use of the native data resolution brings additional challenges to the SED fitting. First, it requires higher S/N observations since the data are used at a higher resolution (more noise per resolution element) and the short-wavelength information of small-scale intensity variations is combined with weaker constraints on the SED shape, which partly depend on observations at longer wavelengths and lower angular resolutions. Second, the fitting can no longer be done pixel by pixel; it has become a global optimisation problem with a much higher computational cost.

In this paper, we use simulations of interstellar clouds to compare the precision of column density estimates calculated with the single-component MBB model and with some multi-component models. For the latter, we examine both the use of SED templates (here MBB functions with fixed temperatures) and the use of MBB functions with temperatures directly as free parameters. The observations are used at their original resolution by including the convolution with the telescope beam as part of the model. In view of the associated increased computational cost, we also pay attention to the organisation of calculations and the potential benefits of using graphics processing units (GPUs). The fits were mainly carried out with direct numerical optimisation, but some tests were also carried out using Markov chain Monte Carlo (MCMC) methods. The simulated observations include simple toy models with two temperature components and more complex cloud models where the synthetic observations were created through radiative transfer modelling.

The contents of the paper are the following. In Sect. 2, we present the multi-component MBB models and discuss some details of the SED fits. Section 3 presents the main results from the tests, which include toy models with ad hoc temperature distributions (Sect. 3.1), Bonnor-Ebert spheres with internal and external heating (Sect. 3.2), a cloud filament with stochastically heated grains (Sect. 3.3), and a large-scale cloud model based on magnetohydrodynamic (MHD) simulations (Sect. 3.4). We discuss the results in Sect. 4 before presenting the main conclusions in Sect. 5.

2 Methods

Dust emission can be approximated as MBB emission, where the observed intensity I_ν is $$I_v=B_v\left(T,\beta \right)\left(1-e^{-{\tau }_v}\right)\approx B_v\left(T\right){\tau }_v=B_v\left(T\right)\text{Σ}{\kappa }_v\left(v_0\right){\left(\frac{v}{v_0}\right)}^{\beta },$$(1)

where Τ is in principle the dust temperature but more precisely the colour temperature that depends on all temperatures along the LOS. In Eq. (1) the optical depth τ is written as the product of mass surface density Σ and the dust mass absorption coefficient κ_ν, which is further assumed to follow a power law with the opacity spectral index β. The constant ν₀ is an arbitrary reference frequency. In the following, Eq. (1), the modified black-body function, is written as MBB(T, β).

When the emission originates from a source with a distribution of temperatures, the intensities may still be approximated with a sum of MBB function, $$I_v\approx {\displaystyle \sum _{i}W\left(T_i\right)\text{MBB}\left(T_i,{\beta }_i\right)},$$(2)

where the weights W(Τ_i) are the fitted parameters. We use two versions of Eq. (2). In the N-MBB models, the fit is based on directly on the above equation at each map position r̄, $$I\left(\overline{r},v\right)={\displaystyle \sum _i^N{W}_i\left(\overline{r}\right)}{\text{MBB}}_v\left(T_i\left(\overline{r}\right),{\beta }_i\left(\overline{r}\right)\right).$$(3)

The model includes N MBB functions, with W_i(r̄) and Τ_i(r̄) as the free parameters while the spectral index β_i(r̄) is kept constant. With N=1, Eq. (3) describes the standard single-component MBB fit (1-MBB).

The alternative SED models N-TMPL correspond to $$I_v={\displaystyle \sum _i^N{W}_i\left(\overline{r}\right)C_i\left(v\right)}.$$(4)

Here C_i(ν) are fixed spectral templates. These can be MBB functions for some fixed values of Τ and β, the corresponding functions for some distributions of Τ and β, or any other templates that could be expected to describe the observed emission. In the tests, we set C_i(ν) equal to MBB function at a given discrete temperature. Because the temperatures remain fixed during the fitting, the model is different from Eq. (3) and the scaling factors W are the only free parameters of the N-TMPL models.

In Eqs. (3)–(4), the left side is the observed intensity I_ν, which is the true sky surface brightness convolved with the telescope beam. In 1-MBB models, we follow the standard analysis, where all observations are convolved to a common lower angular resolution, the fit itself does not include convolutions, and each map pixel can be fitted completely independently. This analysis thus results in estimated intensity and temperature maps at the same common resolution. In all other cases, the model and the observations are compared at the original resolution of the observations. This means that the model predictions have to be convolved with the telescope beam on each step of fit. The resolution of the resulting column density maps is now not as well defined. The band with the highest resolution provides information of the sky brightness at that frequency and that resolution, while at least two bands are needed to trace temperature variations. The sensitivity to emission from dust at different temperatures varies between the bands, making the effective resolution of the resulting τ map less precise.

When convolution is part of the optimisation, the fit minimises the weighted least squares errors $${\chi }^2={\displaystyle \sum _k{\left(\frac{I_k-M_k\otimes B_k}{\delta I_k}\right)}^2,}$$(5)

where k refers to the frequency band, M_k is the unconvolved model prediction (right-hand side of Eq. (3) or Eq. (4)), B_k is the telescope beam, and ‘⊗’ stands for convolution. The problem is computationally harder, since the model-predicted surface brightness maps are convolved with the telescope beams on every step of the optimisation and the fit parameters for the whole map are connected via the convolutions. The convolution can be done in real space, with direct summation over the beam footprint, or using Fourier transformation (fast Fourier transforms, FFTs). We use mainly the latter option, which is faster, but more limited in its ability to weight individual measurements or to handle missing values and map boundaries.

When beam convolutions are part of the fitting procedure, all parameters for the entire map are combined into a single optimisation problem. The large number of data points precludes the use second order information (i.e. Hessian matrices), which would lead to excessive memory usage. We use mainly conjugate gradient algorithms, where only the gradient of the χ¹ function is needed (relative to each of the fitted parameters). The problem remains feasible, because the telescope beams couple data directly only over small distances. Therefore, although all fitted parameters are in principle connected, this connection is strong only over short distances.

If we denote the observed surface brightness values with S, the χ² function of N-TMPL models becomes $${\chi }^2={\displaystyle \sum _{i,k}{\left(\frac{S_{i,k}-{\left[{\text{Σ}}_{\text{c}}{W}_c{C}_{c,k}\otimes B_k\right]}_{i,k}}{\delta S_{i,k}}\right)}^2,}$$(6)

where the indices i, k, and c refer to the position, the frequency band, and the SED template, respectively. Since the T_C values are pre-selected constants, the number of free parameters (W_c) is the same as the number of the fitted SED components. For the N-MBB models the minimised function is $${\chi }^2={\displaystyle \sum _{i,k}{\left(\frac{S_{i,k}-{\left[({\text{Σ}}_{\text{c}}M\left(I_{\text{c}},T_{\text{c}},{\beta }_{\text{c}}\right)\otimes B_k\right]}_{i,k}}{\delta S_{i,k}}.\right)}^2.}$$(7)

Here M stands for the MBB function, and the term in square brackets denotes the convolved model-predicted map at one frequency. The index i is omitted within the square brackets, but the result of the convolution is a function of position (index i) and frequency (index k). When β is kept constant, the number of free parameters is two times the number of SED components.

For better performance and to guarantee sufficient accuracy of the gradient estimates, it is important that the gradients are calculated analytically (Appendix A).

By default, the optimisation of N-TMPL models is performed in real space. However, the templates C_c,f appear in Eq. (4) only as constant factors, and this makes it possible to move the optimisation to Fourier space. If the Fourier transformation is denoted by ℱ, one minimises $${\chi }^2={\displaystyle \sum _{i^{\prime },k}\frac{{\left[ℱ\left(S_{i^{\prime },k}\right)-{\text{Σ}}_{\text{c}}{C}_{\text{c},k}ℱ\left(W_{i^{\prime },\text{c}}\right)\right]}_{i^{\prime },k}^2}{\text{d}{S}_k^2}},$$(8)

where the index i′ now denotes an element in the Fourier space. The free parameters are the real and imaginary parts of the Fourier amplitudes that correspond to the real-space weight maps W_i,c Because the elements ℱ(W_i,c) result from the Fourier transformation of a real-valued function, only n + (n − 1)N nonredundant elements need to be optimised, where N is the map size and n = [N/2]. Within the optimisation, the costly convolution is replaced with a direct multiplication. For the largest test cases, this amounts to a speedup of one or two orders of magnitude. If the observational noise in S is Gaussian, the errors in the Fourier amplitudes are also normally distributed, justifying the use of the least-squares formula. One important restriction in Eq. (8) is that the error estimates dS_k are now constant values for each frequency, not being able to be varied pixel by pixel. Compared to the alternative hypothesis of constant relative uncertainties over the maps, this leads to differences when a single beam contains widely different intensity values (changing locally the relative weight of the different pixels) or when the shape of the observed spectra changes significantly (changing locally the relative weighting of the different bands). It is also more difficult to implement general priors or even to constraint optical depths to positive values (i.e. to enforce W_i,k ≥ 0). In the following, the fits performed in Fourier space are referred to as N-TMPL-F fits.

3 Results

We compared the performance of the SED models for different simulated observations. Results are presented for simple toy models (Sect. 3.1), a Bonnor-Ebert sphere, (Sect. 3.2), a filament that includes the emission from stochastically heated grains (Sect. 3.3), and finally for a more realistic large-scale magnetohydrodynamic (MHD) simulation (Sect. 3.4).

3.1 Toy models

3.1.1 Unconstrained fits

To get basic understanding of the performance of the N-TMPL and N-MBB models, we started with simple toy models. In Fig. 1a, the observations consist of two emission components, one at a T₁ = 15 K and another at a temperature that changes linearly over the 128 × 128 pixel maps from T₂ = 10 K to T₂ = 20 K. The data include 160, 250, 350, and 500 μm maps with 3% relative observational noise. The beam size was set in all bands equal to half of the pixel size, to eliminate in this test the effect of the beam convolutions. This enables more direct comparison between the 1-MBB and multi-component models, when only the latter have the beam convolutions as part of the fitted model. The 2-TMPL model uses MBBs at T = 13 K and T = 17 K as the SED templates. The data were initially analysed assuming β = 1.9, the same value that was used to generate the synthetic observations. Figure 1 plots the quantity τ/τ₀, the ratio between the estimated and the true optical depths.

The 1-MBB model behaves in Fig. 1a as expected. The correct optical depth τ is recovered when the emission consists of a single temperature (T₂ ~ T₁ = 15 K), and τ is increasingly underestimated as the difference between the two temperature components increases. Since the model has only two free parameters, the scatter of the optical depth estimates is small.

The 2-TMPL model also has two free parameters (W₁ and W₂) and the τ scatter is similar to that of the 1-MBB model. However, the errors increase rapidly away from the two selected T_C values and are positive between the two T_C values and negative outside the T_C range.

The construction of the simulated data in Fig. 1a matches the assumptions of the 2-MBB model, which also provides practically unbiased results. However, the number of free parameters is four, and the statistical noise of the τ estimates is two times larger than in the 1-MBB and 2-TMPL cases. The absolute accuracy of the individual τ estimates is therefore on average not better than for the 1-MBB model.

Figure 1a quotes the χ² values of the fits. In Fig. 1 and later, we always quote the χ² values normalised by the number of fitted data points.

The values are not scaled with the degrees of freedom (which is not well defined for non-linear models or models containing constraints; Andrae 2010) and thus directly compares how well the models match the observed intensities. A smaller χ² value (a good match to intensities) does not necessarily mean more accurate τ estimates. The 2-TMPL fit was done with MCMC calculations, and the plotted values are averages over the MCMC samples. We also plot the 1-σ confidence region (containing 68% of the MCMC samples), which follows very closely the actual dispersion of the τ estimates. However, these do not reflect the true error, which is dominated by systematic errors.

Figure 1b uses a different set of synthetic observations. The optical depths are distributed in each pixel at different temperatures according to a Gaussian N(T₀,σ_T). In the plot, σ_T is 1K for the first half of the x-axis and 5 K for the second half. The mean temperature 〈T〉 increases in both halves linearly from 12 K to 20 K. The temperature distributions are also truncated to T > 8 K. The 1-MBB model performs well, with 20% bias reached only at the lowest temperatures and with the wider σ_T 5 K temperature distribution. Most 1-MBB estimates remain within the 1-σ error band of the 2-MBB model that has larger statistical noise. Like in Fig. 1a, the systematic error of the 2-TMPL model are significant. This is true especially at the highest temperatures, where the dust temperature distribution extends clearly above the T_C values.

Figure 1c uses the same observations as frame b but the fits assume a lower value of β = 1.7, which thus leads to lower τ estimates. The drop affects the 2-MBB fit slightly more than the 1-MBB fit, making the latter the most accurate in terms of the τ estimates. The wrong β value also exacerbates the 2-TMPL errors at high temperatures, where the estimates now fall even below zero. The unphysical values are the result of the model trying to match the SEDs (with temperatures outside its T_C range) using a negative weight for the colder component and a larger positive weight for the warmer component. It is thus clear than the over 10 K range of dust temperatures cannot be satisfactorily fitted with just two T_C components.

Fig. 1 Ratio of estimated and true optical depths τ/τ0 for 1-MBB, 2-TMPL (TC=13 K and 17 K), and 2-MBB models. The observations consist of 160, 250, 350, and 500 μm intensities with 3% relative noise. The dots correspond to individual pixels (every 10th pixel), the solid lines show moving averages, and the shaded regions the 1-σ dispersion of the τ estimates. The χ1 values of the ML fits (normalised by the number of fitted data points) are quoted in the frames. In frame a, observations consist of one emission component at T1 = 15 K and one component varied in the range T2 = 10–20 K (on x-axis), both with the same optical depth. The additional solid blue lines correspond to the 1-σ error estimates from the MCMC fit of the 2-TMPL model. In frame b, the observations correspond to a continuous temperature distribution, T ~ N(〈T〉, σT). The x-axis shows the 〈T〉 value, with σT = 1 K for the first and σT = 5 K for the second half. Frame c is the same as frame b but the SED fits assuming β = 1.7 instead of the correct value of β = 1.9.

3.1.2 Fits constrained with priors

Figure 2 shows the corresponding results, where we include a prior T ~ N(18 K, 5K) for the 2-MBB fits. The intensities of the 2-MBB fits (parameters I_i,c) and the 2-TMPL fits (parameters W_i,c) are restricted to non-negative values by optimising the logarithm of the relevant variables (cf. Appendix A). With these constraints, both models could be successfully fitted using MCMC. There is little difference between the 1-MBB and 2-MBB fits, although the latter has a smaller bias in the case of the σ(T) = 5 K observations. The non-negativity has a clear effect on the 2-TMPL fits, as the optical depths are now overestimated rather than underestimated for temperatures above the T_C range. This shows that the Fig. 1b fit already contained negative components, even when the total τ estimates were positive. Although the unphysical negative values have been eliminated, the χ² values of the fits are now much larger and the systematic errors of the τ values still increase rapidly outside the adopted T_C range.

Figure 3 shows fits for 3-TMPL and 3-MBB models. With T_C = 13, 15.5, and 19 K, the bias in 3-TMPL results can be constrained to ~20%. With a prior of T ~ N(15 K, 3 K), the results of 3-MBB fits are unbiased in the σ_T = 1 K part of the map (left half of the plot) but still underestimate the σ(T) = 5 K observations by ~10%. Even if the systematic errors are partly smaller than in the 1-MBB fit, there is no clear net improvement because of the larger statistical noise.

Fig. 2 Same as Fig. 1, but forcing each component of the 2-MBB and 2-TMPL models to positive optical depths. In the case of 2-MBB. the model also includes a prior on the temperatures, Τ ~ Ν(18 Κ, 5 Κ). The plotted 2-TMPL and 2-MBB data are from MCMC calculations while the χ2 values are the average values of the maximum-likelihood solutions.

Fig. 3 Same as Fig. 2b, but including results for 3-TMPL (MCMC fit) and 3-MBB (ML fit) models. The 3-TMPL model uses TC = 13, 15.5, and 19 K, and the 3-MBB model includes a temperature prior Τ ~ N(15 K, 3 K).

3.1.3 Fits including beam convolution

Figure 4 shows the first comparison where the beam size is different for different bands. The synthetic observations correspond to the sum of two temperature components. The first one is at 15 K, apart from two Gaussian sources where the temperature changes to 10 Κ or 20 K. The second temperature component increases linearly from 10 Κ to 20 Κ over the map. The pixel size is 10″, and the full width at half maximum (FWHM) values of the Gaussian beams 20, 20, 25, and 36″ at 160, 250, 350, and 500 μm, respectively. The observations contain 3% relative noise at the resolution of the each map, and all results correspond to the ML solutions.

In Fig. 4, the 1-MBB model has again the lowest statistical errors, partly because all data have been convolved to a 36″ resolution. The systematic errors exceed 20% when the 15 Κ component is combined with the extended second component at ~10 K. (leftmost part of Fig. 4c), when the warm Gaussian source is combined with cold extended emission (around pixels index ~4500), and even more significantly when the warm 17–18 Κ extended emission is combined with the cold Gaussian source (pixel indices ~ 12 000).

The 2-TMPL model uses values T_C = 13 Κ and 18 K. The statistical noise is larger than for the 1-MBB model but now also corresponds to a factor of two higher angular resolution. The angular resolution of the model is nominally 20″, but because it uses all bands up to the 500 μm map with 36″ resolution, the solution might exhibit weak correlations up to that scale. The τ values are again significantly underestimated, when the emission contains emission from dust below the T_C values. At the location of the cold Gaussian source (10 Κ against the ~ 17–18 Κ background), the bias reaches 50%, the same as for the 1-MBB model. Towards the map centre, the emission corresponds to the single 15 Κ temperature and τ is overestimated by up to 20%.

The results of the 2-MBB model are practically unbiased. However, apart from the effect of convolutions, the synthetic observations also precisely match the assumptions of the 2-MBB model. The dispersion of the τ estimates is larger than for the other SED models and also larger than in the previous examples. Thus, the overall accuracy is only slightly better than for the 1-MBB model, which however also corresponds to a lower resolution. When the convolutions were performed using FFTs, the fits took of the order of 1 min for the 128 × 128 pixel maps.

Fig. 4 SED fits considering the effect of telescope beams. The 2-TMPL and 2-MBB cases include beam convolution as part of the fit while the 1-MBB fit is done at the lowest common resolution. The dust emission consists of the two temperature components shown in frames a and b. The map size is 128 × 128 pixels, the pixel size 10″, and FWHM beam sizes 20″, 20″, 25″, and 36″ at 160, 250, 350, and 500 μm, respectively.

3.2 Bonnor-Ebert sphere

We calculated synthetic surface brightness maps for a 2.6 M_⊙ critically stable Bonnor-Ebert (BE) sphere (Ebert 1955; Bonnor 1956), with an assumed gas temperature of 15 K. However, in the subsequent radiative transfer modelling the cloud is illuminated from the outside, with the Mathis et al. (1983) interstellar radiation field (ISRF), and optionally from the inside by a 10 L_⊙ point source that is modelled as a 10 000 Κ black body. These result in continuous radial dust temperature gradients and each LOS has emission from a different range of temperatures. In the outer parts, the colour temperature approaches 18 K. In the inner parts, the temperature is below 14 K for the externally heated and over 30 K for the internally heated case. The pixel size is set to 6″, with beam sizes 18″, 18″, 25″, and 36″ at 160, 250, 350, and 500 μm, respectively. The subsequent SED fits assume a constant β, which corresponds to the 160–500 μm spectral index of the dust model used in the radiative transfer calculations (Compiègne et al. 2011).

Figure 5 shows results for the externally heated model. The 1-MBB model provides even surprisingly accurate results (at 36″ resolution). The 2-TMPL model uses T_C values of 14.5 K and 17 K. These provide accurate estimates towards the cloud centre but overestimate τ in the outer parts, where the dust temperatures are above the Τ_C values. The 3-TMPL model has Τ_C values of 14, 15, and 17 K. Because the highest T_C is the same as in the 2-TMPL fit, the results are similar in the outer parts of the cloud, but there is minor improvement at the centre (corresponding to the added Τ_C value). The width of the LOS temperature distribution appears to be sufficiently narrow so that the 1-MBB model, by adapting to the LOS mean temperature, still performs well compared to the 3-TMPL model. The 2-TMPL and 3-TMPL fits enforced the positivity of the component weights but used no other priors. Without the positivity requirement, the χ² values would be lower but especially the 2-TMPL result would include negative weights at some radial distances.

The 2-MBB fits in Fig. 5 were also calculated without priors. The estimates are mostly unbiased but show a relatively high scatter. The corresponding noise in the temperature being visible in Fig. 5d,e, where temperatures are not well constrained for individual pixels, at scales below the beam size. The low χ² value also indicates the possibility of some overfitting.

Figure 6 shows the corresponding results for the internally heated BE model. The errors of the 1-MBB fit increase to 20% towards the cloud centre. The selected T_C values are 18 and 25 K for the 2-TMPL model and 18, 23, and 29 K for the 3-TMPL model. The 2-TMPL model is unable to follow the rapid temperature changes in the centre, and the oscillating errors result in some artefacts in the estimated radial density profiles. The 3-TMPL model performs better at the very centre, but the overall accuracy is not much better than for the 1-MBB model.

The 2-MBB model shows increasing bias towards the cloud centre and the τ/τ₀ curve actually follows that of the 1-MBB model predictions. This is partly a coincidence, because the latter correspond to a factor of two lower angular resolution. At the centre of the projected cloud, the 2-MBB fit uses two temperatures that are both close to 30 K and separated from each other by some 4 K. This is still quite small compared to the full range of the actual LOS temperatures, which extend from ~17 K to over 50 K, (the highest temperatures at scales below the beam sizes).

3.3 Filament with stochastically heated grains

We examined the optical depth estimates also as a function of the wavelengths. This has already been studied extensively from the point of view of observational noise and the LOS temperature variations (Shetty et al. 2009a,b; Malinen et al. 2011; Juvela & Ysard 2012b,a; Juvela et al. 2013a,b). Below the main questions are related not only to the relative performance of the N-MBB and N-TMPL methods but also to the potential effect of stochastically heated grains.

We used the filament simulations of Juvela & Mannfors (2023). The model is discretised onto a Cartesian grid of 200³ cells with a cell size of 0.0116 pc. The model consists of a single filament with a Plummer-like column density profile $$N=N_0{\left[1+{\left(r/R_0\right)}^2\right]}^{\left(1-p\right)/2}$$(9)

with R₀=0.0696 pc and p = 3. The column density perpendicular to the filament major axis is set to N(H₂) = 10¹¹ cm⁻¹. The filament is illuminated by a normal ISRF (Mathis et al. 1983). There can be an additional 15 700 K point source with a 590 L_⊙ luminosity, which is is placed 0.93 pc in front of or to one side of the filament. The source is 0.93 pc from the map centre, in a direction parallel to the filament main axis. The dust properties, including those of stochastically heated grains, are taken from Compiègne et al. (2011).

We calculated 200 × 200 pixel maps, where the pixel size corresponds to the 3D discretisation and is set equal to 5″ (corresponding to a distance of 479 pc). Synthetic surface brightness maps were computed at 70, 100, 160, 250, 350, and 500 μm with the assumed angular resolutions of 10″, 10″, 12″, 18″, 25″, and 36″, respectively. The wavelengths and resolutions are similar as in Herschel observations, except for the lower angular resolutions adopted for the 70 and 100 μm bands.

In the absence of the point source, the 1-MBB fits of the 160–500 μm observations give a narrow range of colour temperatures from 16 K to 17.7 K. Based on this, we selected temperatures T_C=16 Κ and 17 Κ for the 2-TMPL and T_C = 14, 16, and 18 Κ for a 3-TMPL fit. When the point source is added on the LOS towards the filament (0.93 pc in front of the filament and 0.93 pc from the projected map centre), the range of 1-MBB colour temperatures becomes 19.5–26 K. If the point source is to one side of the filament, the maximum exceeds 40 Κ but only towards the source position, in a region of relatively low column density. In these cases, we used T_C = 21 and 25 Κ for the 2-TMPL fits and T_C = 20, 23, and 26 Κ for the 3-TMPL fits. All surface-brightness observations contain 3% relative noise.

Figure 7a shows SEDs towards the map centre. These correspond to a single pixel and models without the point source or with the point source on the LOS towards the filament centre (but offset from the 3D centre of the filament). The fits show the systematic increase of the colour temperatures as the 1-MBB fits are extended to higher frequencies. The lower frame shows how the SED fits translate to τ estimates. The optical depths tend to be underestimated, and this also applies to the 2-TMPL and 3-TMPL models. The systematic error is largest for the 1-MBB model, roughly equal for the 2-TMPL and 3-TMPL models, and smallest for the 3-TMPL model. Like the average colour temperature, the bias of the τ estimates tends to increase as shorter wavelengths are included in the fit. The exception are the 2-MBB fits, where the statistical errors dominate. The error bars in Fig. 7b correspond to the τ estimates over a 0.7 pc long section of the filament spine. In models containing a point source, this section also includes the projected location of the point source.

Figure 7b shows that the χ² values per data point of the 2-MBB fits are of the order of one for all band combinations. In the 1-MBB fits, the increase in the number of bands is associated with a rapid increase of χ² values (as indicated by the SEDs in frame a). A similar trend exists for the 2-TMPL model and to a lesser extent for the 3-TMPL model, these of course depending on the chosen T_C values.

Fig. 5 Fits on synthetic observations of an externally heated Bonnor-Ebert sphere. Frames a–c show the mass-weighted average LOS temperature of the model cloud and the calculated 160 μm and 500 μm surface brightness maps. Frames d-e show the temperature maps and frame f shows the 500 μm predictions from the 2-MBB model. Frame g shows the ratio of the estimated and true optical depth (τ/τ0) as the function of τ0 for the four cases listed in the legend. The solid lines are moving averages, the shaded regions indicate the 1-σ interval in the τ/τ0 values, and the individual pixels outside this interval are shown as dots.

Fig. 6 Same as Fig. 5, but for a BE sphere with an additional internal radiation source. The TC values of the 2-TMPL and 3-TMPL models are given in the text. The resolution is 36″ for the 1-MBB model and 18″ for the other fits (for a nominal pixel size of 6″).

Fig. 7 SED fits of the filament models. Frame a shows the SEDs towards the map centre and the 1-MBB fits to the 4–6 longest wavelength bands. The lower set of symbols and curves corresponds to a model without a point source, and the upper ones to a model with a point source on the LOS towards the filament but at 0.93 pc in front of the filament. The colour temperatures of the 1-MBB fits are quoted in the frame. Frame b shows the ratio of recovered and true peak optical depths (left axis, filled symbols) and the average χ2 per data point over the whole map (right axis, open symbols). The x-axis indicates the four models. The square and star symbols are, respectively, for the externally heated model and the model with a point source. The black, blue, and red symbols correspond to fits that extend to shortest wavelengths of 160, 100, and 70 μm, respectively. The error bars show the standard deviation of the optical depth estimates along the filament spine.

3.4 MHD cloud simulation

As the final test we examined a MHD cloud simulation. The full model volume is (4 pc)³ in size, but we used only the sub-volume of (1.84 pc)³ that contains the densest structures. The column densities reach values above N(H₂) ~ 10²³ cm⁻², while the average volume density is still below 10³ cm⁻³. The simulation includes a few hundred embedded sources that increase the temperature variations. Details of the MHD simulation can be found in Haugbølle et al. (2018) and of the radiative transfer modelling in Juvela et al. (2022). Because the model includes a filamentary structure reminiscent of a small infrared dark cloud (IRDC), we refer to the model as the IRDC cloud.

The synthetic observations consist of maps at 160 250 350, and 500 μm. The effect of stochastically heated grains is not included. For the SED ftting, the main challenges are the large map sizes and the temperature variations associated with the cores. The 3D dust temperatures vary from less than 10 K in the cold cores to up to ~100 K in the hot cores, but the range of the observed colour temperatures is naturally smaller. The surface brightness maps have 1875 × 1875 pixels, and the assumed beam sizes are 2.3, 2.3, 3.1, and 4.5 pixels at 160, 250, 350, and 500 μm, respectively. We use the same wavelengths as in previous tests, although the above scaling would correspond to a source distance of only 25 pc for a telescope of the Herschel size. The observational noise is 3% of the observations.

Below we compare the 1-MBB results to 3-TMPL and 3-MBB fits and then, with additional priors, to N-TMPL and N-MBB flts with N > 3. We also look at the N-TMPL flts made in Fourier space separately, because those adopt a different noise model and are more limited in the use of priors. We are also interested in the relative run times.

3.4.1 IRDC cloud and 2-MBB and 4-TMPL fits

Figure 8 compares the 1-MBB ft to the results from 2-MBB and 4-TMPL fts for a 256 × 256 pixel area of the full IRDC maps. Even this small region includes more than two orders of magnitude differences in surface brightness, with colour temperatures ranging from ~ 14 K to over 30 K. The analysis uses a β value that is correct for the 160-500 μm interval. The 4-TMPL fit uses T_C values of 15, 18, 22, and 27 K. The models were fitted without direct priors, except for the temperatures of the 2-MBB fit being limited to values above 8 K.

In the 1-MBB fits, the average ratio between the estimated and true optical depths is 〈τ/τ₀〉 = 0.81 and the standard deviation σ(τ/τ₀) = 15%. The optical depths are always underestimated, and the errors reach ~60% in a region that contains the superposition of warm extended emission and a more compact and colder filamentary structure. The χ² values are mostly well above one and in this case also correlated with the magnitude of the systematic τ errors.

For the 2-MBB and 4-TMPL flts, the morphology of the τ and χ² maps is similar to the 1-MBB case. The bias of the τ estimates is somewhat lower, and also the χ² values are lower, especially considering that they refer to up to a factor of two higher angular resolution. The maximum errors are still almost 50% but they are limited to smaller regions. The T_C values of the 4-TMPL model were chosen so that they result in accurate τ estimates. However, the results are not very sensitive to the precise T_C values, as long as they cover the range of colour temperatures in the 1-MBB fit.

3.4.2 IRDC cloud and 11-TMPL and 4-MBB fits

When N is increased further, the N-TMPL and N-MBB fits require additional constraints. When β is kept constant, a 4-MBB model has eight free parameters and an 11-TMPL model has 11 free parameters, both thus much larger than the number of the observed wavelengths.

Figure 9 shows results for the same area as in Fig. 8. The TC values of the 11-TMPL model are placed logarithmically between 12 K and 40 K. The temperature prior is T ~ N(25 K, 3 K) for both the 4-MBB and 11-TMPL fits. The large χ² values in the upper right hand part show that the fits do not yet correctly describe the complex temperature distributions of this area. Elsewhere in the maps the χ² values are at or below ~1, but low χ² values do not translate to more accurate τ estimates. The τ/τ₀ ratios are now more constant but there is also higher positive bias. The use of different apparently reasonable priors result in 〈τ/τ₀〉 varying by more than 10%. In particular, weaker priors tend to lead to larger temperature fluctuations and larger positive bias in τ. In all cases, the dispersion in the τ/τ₀ ratios is equal or larger than in the 1-MBB case (not accounting for the difference in the angular resolutions).

The same fits for the whole 1875 × 1875 pixel maps are shown in Appendix B. The average statistics are similar to the previous image, although the lower average column density results in lower average bias, especially in the case of the 1-MBB fit. However, there are also regions of stronger internal heating, where all the fitted SED models can underestimate the optical depths by up to 60%.

Fig. 8 IRDC observations fitted with 1-MBB, 2-MBB, and 4-TMPL models. Frames a-b show the 250 μm intensities and colour temperatures from the 1-MBB fit. The second row shows the 1-MBB τ estimates relative to the correct values (frame c) and the χ2 values (frame d). The third and fourth rows show the corresponding data for the 2-MBB and 4-TMPL flts. The plots cover a 256 × 256 pixel area of the IRDC cloud model.

Fig. 9 Fits of IRDC model maps with 4-MBB and 11-TMPL models. The corresponding data for the 1-MBB fit are shown in Fig. 8.

3.4.3 IRDC cloud and N-TMPL fits in Fourier space

Fits of the N-TMPL model in Fourier space (i.e. N-TMPL-F fits) are attractive because of the faster computation, even if it only allows limited use of priors. The error estimates are also required to be constant over the map, and we set them equal to the mean of the 3% relative uncertainties. We also add to the minimised χ² function a penalty $$P=\xi {\displaystyle \sum _i{W}_i^2}.$$(10)

This biases the solution against individual high W_i and, more importantly, makes it unlikely that an observation would be fitted as a sum of positive and negative terms instead of a smaller number of smaller positive components.

Figure 10 compares the 1-MBB results to the 4-TMPL-F and 12-TMPL-F fits. The T_C values were placed logarithmically in the range 11–45 K. The fits were done without direct priors on the temperature. However, the 12-TMPL-F fit included a small penalty according to Eq. (10), with a ξ values such that the effect on the final χ² was ≲ 1%.

The 1-MBB results are identical to the previous tests, apart from the difference in the adopted error estimates, and still underestimate the true optical depth by up to 60%. Although the differences are not very large, the 4-TMPL and 12-TMPL fits are here more accurate, with a more symmetric error distribution around the correct value. The average 〈τ/τ₀〉 value of the 12-TMPL-F fit varies by ±0.1 units if the number of components is changed between 10 and 20, the lower limit of the T_C values between 10 K and 14 K, and the upper limit between 30 K and 40 K. The 4-TMPL-F fit is also relatively robust with respect to selected T_C grid. The 10% and 90% percentiles of the τ/τ₀ ratios (shown in Fig. 10) are very similar for the 4-TMPL and 12-TMPL fits, in spite of the factor of two difference in the χ² values.

The upper right corner includes a region where the 4-TMPL-F and 12-TMPL-F fits overestimate τ by almost -60%. If the T_C grid is extended to 55 K, the maximum errors of the 12-TMPL-F fit are confined to a more symmetric range of [−35%, +31%]. This is in agreement with the maximum LOS temperatures being higher than the maximum 1-MBB colour temperatures of ~35 K. When the highest T_C is 35 K, some components around T_C ~ 30 K have negative weights in the area of the largest τ errors, in spite of the use of Eq. (10) penalty. The solution is thus partly unphysical, although, to the extent that the emission is in the Rayleigh-Jeans regime, this does not translate to large errors in the τ estimates. When the T_C grid is extended to higher temperatures, practically all weights remain positive.

The fits performed with the same parameters for the whole 1875 × 1875 pixel maps are shown in Appendix B. Due to the lower average column density, the mean bias of the 1-MBB estimates is there closer to zero (Fig. B.2), and also the 10–90% intervals of the τ/τ₀ ratios are similar between the fits. There is a number of hot sources where τ is underestimated by all three methods.

4 Discussion

Most studies model the FIR dust emission with a single MBB function (1-MBB), and the observations (the noise level and number of observed wavelengths) also may not allow much more complex analysis. The main problem of the single-temperature model is the underestimation of column densities (Shetty et al. 2009b; Juvela & Ysard 2012a). The bias depends on the width of the temperature distribution, which makes SF regions with hot and cold cores particularly challenging. We have examined in this paper the use of multi-component SED models, using observations at their native resolution. We discuss the main findings below.

4.1 Basic tests

When the number of free parameters is less than the number of observed wavelengths, the SED fits can be done without additional constraints. The use of prior distributions may still be beneficial, especially when the data have low S/N.

In the tests of Sects. 3.1 and 3.2, the 2-MBB fits provided only small improvement over the 1-MBB fit. They resulted in lower bias in the τ estimates but with a larger statistical noise. The χ² values were not directly correlated with the precision of the τ estimates. The low number of SED components was more of a problem for the 2-TMPL model, because the T_C values remain fixed and, unlike in the N-MBB models, not adjusted for each LOS separately. The errors depended on the distance between the T_C values and the actual source temperatures and increased rapidly outside the T_C range. This is a problem, since the full range of temperatures may not be known a priori and is difficult to estimate from observations. On the other hand, the errors may have been exaggerated in Sect. 3.1, because of the discrete temperature values of the synthetic observations. In real sources the temperature distribution is always wider, leading to less sharp changes in the bias.

In Sect 3.2, the internally heated Bonnor-Ebert sphere could not be fitted properly with just two temperature components. The absolute errors of the 3-TMPL were below ~20% but showed rapid oscillations as the mean LOS temperature moved across the T_C grid. When the true solution is not available, such oscillation cannot be recognised, and this could have a significant effect on the derived density profiles. When beam convolutions are included in the SED models, the 2-MBB model also showed some problems with temperature oscillations at scales below the beam size. When the predictions of such a model are convolved to the observed resolution, the SED fits may be perfect but the predicted column densities will show positive bias that depends especially on pixels with the coldest temperatures. This might require additional constraints on the smoothness of the solution, provided that the temperature field in the source is indeed assumed to be smooth at the observed resolution. Such priors would also affect the effective angular resolution of the obtained τ maps.

4.2 Filament with stochastically heated grains

In Sect. 3.3, we examined the dependence of τ estimates on the observed wavelengths in the 70–500 μm range. It is well known that the inclusion of shorter wavelengths tends to increase the estimated temperatures (e.g. Shetty et al. 2009b). This is based on the warmer dust having a higher contribution (relative to its mass) to the observed intensities. We had additional effects from stochastically heated grains that contributed significantly only to the shortest wavelengths. In a strong radiation field, the 70 μm emission could be emitted by warm large grains, but in our model the stochastically heated grains dominated the 70 μm observations and had a clear contribution to the 100 μm emission.

When the fitted wavelength range was extended from 160– 500 μm to 70–500 μm, the optical depths were increasingly underestimated. The error was always more than 20% in the 1-MBB fits and mostly around 20% for the other fits. The 2-MBB fit was an exception, showing no clear bias (partly because of the larger statistical errors). There was no large difference between filaments that were illuminated by an isotropic radiation field or also by a nearby star. The main effect is that a higher radiation field results in higher temperatures, which always tends to decrease the systematic errors by moving the observations further into the Rayleigh-Jeans regime.

Fig. 10 IRDC observations fitted in Fourier space. Frame a shows the true optical depth and frame e the colour temperature from the 1-MBB fit. Frames b–d show the optical depths estimated with the 1-MBB, 4-TMPL, and 12-TMPL fits, respectively. The corresponding ratios of estimated and true optical depths are shown in frames f-h and the χ2 values in frames j-l. Frame i shows the standard deviation of the temperature values from the 12-TMPL fit. The values quoted in the frames are the 10%, 50%, and 90% percentiles of the plotted quantity. The maps cover the same 256 × 256 pixel area as in Fig. 8.

4.3 The IRDC model

The IRDC model represented a complex star-forming cloud with a large range of temperatures and column densities. Therefore, the 1-MBB fits showed significant errors with the optical depths being underestimated by up to 60%. The 2-MBB and 4-TMPL models resulted in some improvement, although this required a careful selection of the Τ_C values for the 4-TMPL model. However, as long as the number of free parameters is not larger than the number of observed wavelengths, the solution is more or less well constrained, and there is some benefit from using SED model that describe not only the average temperature but also the temperature variations within the beam.

We also tested SED models with a larger number of free parameters. These did not produce clear benefits (apart from lower χ² values). On the contrary, these require stronger priors that risk biasing the results. The noise in multi-component models tends to increase the effective temperature dispersion of the SED model, which always biases the τ estimates upwards. The run times scales for both N-TMPL and N-MBB models approximately linearly with the number of SED components. The N-TMPL fits can be done faster in Fourier space but, as discussed in Sect. 2, the method has its limitations. The error estimates are (at least in our implementation) limited to a single value per, instead of individual error estimates for each pixel and frequency. One also cannot use direct priors for the weights W_i, although it is possible to give priors for the Fourier amplitudes, as was done in Sect. 3.4.3 to bias against negative weights. In Fig. 10, both the 4-TMPL and 12-TMPL fits resulted in clear improvement over the 1-MBB results. However, when the same parameters were used to analyse the full 1875 × 1875 pixel map, the improvements were less clear. The differences between the different SED models can in itself be a useful diagnostic, giving indications of the model errors that are not covered by the normal fit error estimates.

4.4 Degeneracies in multi-component SED fits

One factor affecting the systematic and statistical noise of the multi-component fits is the beam convolution. The multi-component fits N-MBB and N-TMPL try to estimate Τ and τ at the highest observed resolution. Since the resolution varies with wavelength, the best intensity information is available on a higher resolution than the information on the SED shape. There are even weaker constraints for the SED model parameters in individual map pixel, because only the averages over the beam are directly constrained. This was most directly visible in the N-MBB fits as the small-scale noise.

Even without the added complication of beam convolution, multi-component fits suffer from degeneracies between the fitted parameters. This can be demonstrated easily even for observations of a single pixel. Figure 11 shows the SED for a LOS dust temperature distribution that consists of Gaussian components at 12 K and 17 K. Figure 12 shows the corresponding distributions of the estimated optical depths for 100 noise realisations. The average bias of the resulting τ estimates is 25% for the 1-MBB model, 20% for the 2-MBB model, and 10% for the 2-TMPL models. If 2-MBB fits are performed without priors, a significant fraction of fits include one component with very low temperature (fitting some noise feature), which then leads to τ being overestimated. The error can reach orders of magnitude, if temperatures T ≪ 10 Κ are allowed. In Fig. 11, we use a temperature prior T ~ Ν(15 Κ, 3 Κ), which is still not enough to prevent occasional 50% overestimation of τ. On the other hand, if we were observing a field that potentially includes cold cores, the above prior would bias their τ estimates in the other direction, leading the optical depth of cold cores to be underestimated.

The 2-TMPL model performs well in Fig. 11, because the T_C values were chosen to match the peaks of the true temperature distribution. The results change only little if the T_C values are both moved 1 Κ higher or lower. However, the ratio 〈τ/τ₀〉 drops to 0.80 if the T_C values are both moved 1 Κ further apart, and increases to 〈τ/τ₀〉=1.16 if the T_C values are both moved 1 K closer to each other (following the trends seen in Sect. 3.1). The most accurate solution also corresponds to the lowest average χ² value, 〈χ²〉 = 0.5 (compared to 〈χ²〉 ~ 1 for the wider and the narrower T_C separations). For all the T_C combinations mentioned above, the χ² distribution (over different noise realisations) extends from χ² ~ 0 to χ² ~ 2 and above. Therefore, for an individual measurement, the best model cannot be reliably selected based on the χ² value alone. While 2-TMPL and 2-MBB may have smaller bias, they are not a priori much more accurate in predicting correct τ values. For a collection of observations, the data may contain enough information for the selection of the most appropriate priors (manually or by using hierarchical statistical models). The observations should thus conform to a common, preferably narrow range of SED parameters. Otherwise the effect of priors will be diluted or, alternatively, they may bias the results for sources outside the main parameter distribution.

When the number of observed wavelengths is small, both the 2-TMPL and 2-MBB fits can result in χ² values below one. This suggests that the τ estimates may be significantly affected by the noise realisation (as well as any potential priors). The results of Sect. 3.1 showed that a low χ² value even of the 1-MBB model is not a guarantee of accurate column density estimates. When the number of free parameters is increased, values χ² < 1 can be obtained with different parameter combinations. This is illustrated by Fig. 13 that is based on one noise realisation from Fig. 12. This is worse than an average case but illustrates the potential for parameter degeneracies. Figure 13a shows the χ² values as a function of the temperatures of two components. The intensities of the components are optimised separately for each pixel in the figure (i.e. for every temperature combination). The plot is reminiscent of the degeneracy of Τ − β fits, with long curved valleys where the minimum χ² values remain almost constant (cf. Chastenet et al. 2021). Figure 13b shows that fits with χ² = 1 correspond to a range of solution, where the difference between the highest and lowest τ estimates is more than a factor of two. In 2-MBB fits, any of these temperature combinations could be chosen, depending on the priors. In 2-TMPL fits, the solution would be directly dictated by the pre-selected T_C values, which might also fall completely outside the χ² = 1 contour.

Based on the above, component degeneracies can be a problem already in a two-component model. In multi-component models, these risks are only likely to increase, and feasible solutions with χ² ~ 1 may be reached with many parameter combinations that correspond to different optical depths. In Fig. 14, we simulate 160–500 μm observations that correspond to a single Τ = 15 Κ MBB spectrum or the sum of two MBB components at 13 Κ and 17 K. We examined in this case the use of the 10-TMPL model where the T_C values are set logarithmically between 8 Κ to 27 K. We generated random weights for the ten components and registered the combinations that fitted all SED measurements to better than 5% accuracy. Because this is a condition for the maximum error in any of the four bands, it is more strict than the assumption of 3% relative uncertainty that was used in the previous tests. The detailed shapes of the τ/τ₀ distributions (shown in frames c and f) depend on the way the simulation was done. Nevertheless, they show that the recovered τ values are almost always above the true values. The optical dept can be underestimated only in frame f, when the fit relies on T_C components that fall between the true temperatures of 13 Κ and 17 K. This is another illustration of the obvious fact that four measurements cannot determine 12 parameters in a unique way. Conversely, Fig. 14 can be interpreted to show how the observed sources could have up to a factor of two differences in the optical depths, in spite of showing observationally identical SEDs.

Strict priors can be set only if the nature of the object is well known. If one uses a prior on temperature, such as Τ ~ N(T₀,δT) in some tests in this paper, the parameters T₀ and δΤ will clearly be more appropriate for some parts of the maps than for others. Even if T₀ is varied spatially (e.g. following the temperature from the 1-MBB fit), the parameter δΤ still has a clear effect on the solution (a small δΤ leading to a lower τ). As noted above, the significance of the model errors means that the absolute accuracy of a column density measurement (match to τ₀) cannot be judged based on the χ² values (match to observed intensities). Based on synthetic Herschel 70–500 μm observations of a cloud simulation, Koepferl et al. (2017) also noted that fits with low χ² values can correspond to very inaccurate estimates of temperature and column density and, in their case, the fits could entirely fail close to hot sources.

The PPMAP method (Marsh et al. 2015) is another multi-component SED fitting method that uses a combination of SED templates. These can correspond different temperatures (typically a number larger than the number of observed wavelengths) and optionally also to multiple β values (e.g. Howard et al. 2019). The model is thus in some respects similar to N-TMPL and should be subject to some similar problems of parameter degeneracy. The PPMAP model can use priors for example for the temperature. Furthermore, the program looks for a solution that fits the observations to a χ² ~ 1 accuracy with the smallest number of objects in the state space. In a field like the IRDC model it might at each position use only a couple of components that best match the local temperature distribution, also avoiding excessive fitting of minor noise structures in the SEDs. On the other hand, as discussed above, a low number of (SED) components is not a guarantee of unbiased results and the sanies χ² value may be reached with different parameter combinations with different predictions of the optical depth. All real source also have continuous temperature distributions, so the number of components needed (and thus also the resulting τ estimates) will depend on the observational noise.

Figure 15 shows examples of the τ distributions estimated with MCMC fits. The data are a subset of the IRDC model maps, and the fitted models are 4-TMPL and 2-MBB with the same parameters as in Sect. 3.4. All data are convolved to the same resolution (18″ for an assumed 6″ pixel size) to allow faster calculations where each map pixel is fitted separately. These MCMC fits (4-TMPL and 2-MBB models) of the 256 × 256 pixel maps took each less than 10 s (with 40 000 MCMC steps after the burn-in phase). Therefore, as long as convolution is not part of the model, MCMC remains feasible even for the largest maps. For the assumed 3% observational noise, the typical uncertainty of an individual τ measurement is less than 50%. However, in those map positions, where the fits underestimate or overestimate the true values the most, the true value is far in the tail of the probability distribution. This is natural since even MCMC gives error estimates only within the confines of the selected model. The uncertainties do not cover the model errors, whether the selected model is appropriate or unique in describing the given data.

In this paper, β was mostly fixed to the value that was used in generating the synthetic observations. If β were included as a free parameter in the SED fits, the formal uncertainties of the temperature estimates would increase significantly due to the well-known anti-correlation between the Τ and β parameters (Shetty et al. 2009a; Kelly et al. 2012; Juvela et al. 2013b). Low χ² values would be concentrated in curved regions in the (Τ, β) plane, where the observational noise could cause large changes in the best-fit parameters. The large uncertainty in Τ will directly translate into large uncertainty in τ. The χ² plane could even exhibit distinct minima that correspond to significantly different temperatures and column density estimates (Juvela & Ysard 2012b). The χ² values of the N-TMPL and N-MBB fits could also exhibit multiple minima, even if β we kept constant.

Fig. 11 Simulated single-pixel observation. Frame a shows the optical depth distribution as a function of temperature, and frame b shows the SED, with measurements at 160, 250, 350, and 500 μm.

Fig. 12 Results from the fitting different noise realisations of the Fig. 11 observations. Frame a shows a few examples of the observations (black crosses) and individual fits (solid curves). Frame b–d show the distributions of τ estimates (100 noise realisations) for the 1-MBB, 2-MBB, and 3-MBB models, respectively. Each frame quotes the mean value and the standard deviation for the estimates relative to the true value, τ/τ0.

Fig. 13 Values of χ2 (frame a) and τ/τ0 (frame b) as functions of the assumed temperatures T1 and T2 of the fitted SED templates. The observations correspond to one noise realisation selected from the simulations in Fig. 12. Frame b includes contours at χ2 = 0.1 and χ2 = 1.

Fig. 14 Fits to simulated observations with a 10-TMPL model. The observations consist of 160, 250, 350, and 500 μm intensities for a single Τ = 15 Κ MBB function (upper frames) or a sum of 13 Κ and 17 Κ components of equal optical depth (red circles in frames b and e). Frames a and d show random realisations of the combinations of weights W that result in less than 5% errors for all four predicted intensities (curves in frames b and e). The frames c and f show the distributions of the τ values predicted by those fits. The true optical depth is τ(250 μm) = 0.01.

Fig. 15 IRDC model observations fitted with MCMC methods. Frame a shows the true τ(250 μm) optical depths, and frame d the colour temperature from the 1-MBB fit. Frames b and c show the τ maps from the 4-TMPL and 2-MBB fits (averages over the MCMC samples) divided by the true optical depth τ0. Frames e and f show the posterior probability distributions for the four positions marked in the b and c frames, and the true values τ0 are indicated with vertical lines of the same colour (black, blue, red, and cyan).

Fig. 16 Empirical correction of the 1-MBB optical depth estimates of the IRDC cloud model. Frame a shows the ratio τ/τ0 between the original estimates and the true optical depths. Frame b shows the quantity Σ(ΔIv/Iv) that characterises the difference between each observed spectrum and the 1-MBB fit. Frame c shows the correlation between this quantity and the τ error, and Frame d shows the τ/τ0 ratio corrected for the trend in frame c. In frames a and d, the inserts show histograms of the τ/τ0 values (with linear y-axis). In frame c, the colour scale corresponds to logarithmic point density.

4.5 Correction to MBB fits

Instead of relying on possibly ill-constrained multi-component models, one possibility is to start with the 1-MBB fit and try to mitigate its bias. Juvela et al. (2015b) carried out radiative transfer (RT) modelling of dense clumps and analysed the resulting synthetic surface brightness maps with the 1-MBB method. The ratio between the true optical depth and the optical depth estimated from the synthetic observations was then used to correct the 1-MBB τ estimates. The accuracy of the correction depends on how well the RT model matches the observed target, and this approach is therefore mainly limited to sources with a simple structure. If the model matched observations perfectly, it would directly provide good estimates for the column density, without the need for any further SED analysis. However, it is more robust to use RT models only to derive fractional corrections that are then applied to 1-MBB results. If the target is spherical (or otherwise geometrically simple), methods like the Abell transformation (Roy et al. 2014) provide a more straightforward and computationally efficient way to estimate the temperature and column density structure of a source. However, a full RT model still has the advantage of presenting a physically self-consistent description of the source, where the temperatures cannot vary in an arbitrary fashion but must be consistent with the radiation sources, the density field, and the resulting attenuation of the radiation field. However, the increased complexity also means that it is more difficult to find a radiative transfer model that would be able to match all observations to the full precision of the observations (Juvela et al. 2013a).

One can ask, whether the observations themselves contain enough information to go beyond the simplest 1-MBB fits. If the measurements do not sufficiently constraint the solution, there will always be significant systematic errors that depend on the construction of the SED models.

The IRDC cloud of Sect. 3.4 is the most complex model examined in this paper, containing large temperature and density gradients. Figure 16a shows again the 1-MBB results where the τ estimates are in many regions only 50% of the correct value and the maximum errors are close to a factor of five. The main effect of temperature variations is to make the SED broader. With relative differences between the observed intensities and the 1-MBB model fit, $r\left(\lambda \right)=\left(I_v^{1-\text{MBB}}-I_v^{\text{observed}}\right)/I_v^{\text{observed}}$, we define a parameter $$\text{Σ}\left(\text{Δ}{I}_v/I_v\right)=r\left(160\text{\hspace{0.17em}μm}\right)-r\left(250\text{\hspace{0.17em}μm}\right)-r\left(350\text{\hspace{0.17em}μm}\right)+r\left(500\text{\hspace{0.17em}μm}\right),$$(11)

which becomes more negative for wider SEDs. Figure 16b,c shows that Σ(ΔI_v/I_v) is well correlated with the bias of the 1-MBB τ estimates. In Fig. 16d, the τ estimates have been corrected for this trend, using a third order polynomial fit to the data in Fig. 16c. This simple empirical correction is successful in eliminating most of the bias, although there is only small improvement towards hot sources. The correction is also more sensitive to noise (being a second-order correction to the 1-MBB fit), which reduces the improvement in the accuracy of individual measurements (as suggested by the dispersion in Fig. 16c). The example shows that the observations do contain more information than used of by the 1-MBB fit alone, but, with the assumed 3% observational noise, there also are limits to the potential improvement.

We tested also the use of more general machine learning methods to derive the corrections. A 246 × 246 pixel part of the IRDC maps was used as the training set for a small, four layers deep neural network. We omitted from the maps only 5-pixel wide border regions that are affected by some edge effects from the beam convolution. The training could be accomplished in a time that is comparable to the previous direct multi-component SED fits, and the trained network can then be used to quickly correct the 1-MBB estimates.

Figure 17 compares the accuracy between the original 1-MBB estimates and the neural-net corrected 1-MBB-NN estimates. For the 246 × 246 pixel training set, the correction naturally works well. However, it is not perfect, since Fig. 17b still shows faintly the pattern that dominated the errors in the original 1-MBB fits. The trained network was then used to correct data for the full 1875 × 1875 pixel map (Fig. 17c,d). The result is rather similar to the empirical correction in Fig. 16d. While in Fig. 16 the correction was estimated from the same map where it was then applied, in Fig. 17 the correction was derived from a very small subset (~1.5% of the pixels, the rectangle marked in Fig. 17c). The corrected τ values have a symmetric error distribution. Although there is some positive bias (τ/τ₀ ~ 1.07 + 0.09), the maximum errors have been significantly reduced. This is encouraging regarding the possibility of finding useful corrections, possibly tuned for a given type of sources, that are general enough to be applied also to new observations.

4.6 Computational cost

The computational cost of the different SED models varied by about four orders of magnitude. Once the data are convolved to the resolution of the longest wavelength observations, the 1-MBB fits of the 1875 × 1875 maps could be fitted in just a couple of seconds (Appendix A). Even the empirical corrections of Sect. 4.5 would add little to the run times. This can provide a more robust way to estimate optical depths, considering the potential degeneracies in the multi-component fits, although the angular resolution is limited to the lowest common resolution.

When the beam convolution is part of the fit, the run times increase by orders of magnitude. For the N-TMPL models only, the optimisation can be done faster in Fourier space (Appendix A). In the case of IRDC maps (1875 × 1875 pixels and four wavelengths), the run times were of the order of 1 min instead of the order of 1 h for the fits optimised in the real space. The run times also depend on the chosen convergence criteria. The χ² values tended to decrease initially fast, with only one or 2% further improvement over a larger number of further iterations.

The beam convolution is a major part of the computational cost. For small maps of up to ~256 × 256 pixels, the use of a GPU for the FFT computations provides no benefits, but for the largest maps (1875 × 1875 pixels) the speed-up was a factor of a few. The use of FFTs can cause artefacts at the map boundaries, due to the assumed data periodicity. It is possible to carry out the convolutions in real space, by direct summation of the pixel values over the beams. This gives more freedom in the use of measurement uncertainties (including missing values and potentially even error correlations) and in the handling of the map edges. This remains a viable option for smaller maps, especially when GPUs are used.

The MCMC fitting is typically more time consuming than directs χ² minimisation. On the other hand, it could provide more complete information of the uncertainties than just using the local χ² gradients. Some examples of the use of MCMC calculations were given in Figs. 1–3 and in Sect. 4.4, where each map pixel was still fitted separately. However, MCMC turns out to be viable for maps of moderate size even when the beam convolution is part of the model optimisation. Appendix C shows one example, where 320 × 320 pixel maps are fitted with the 3-TMPL model, the MCMC runs taking about 10 min on a GPU. The MCMC error estimates are consistent with the observed scatter in the τ estimates. However, we emphasise again that these formal error estimates, whether based on the χ² values or MCMC sampling, do not account for the systematic errors that were usually the dominant error component in our tests.

Fig. 17 Correction to 1-MBB τ estimates calculated with neural nets. Frame a shows the 246 × 246 pixel map that is used as the training set. Frame c shows the 1-MBB results for the whole IRDC map, where the location of the training map is indicated with a dashed box in the upper left corner. Frames b and d shows the corrected τ maps for the training set and for the full 1875 × 1875 pixel map, respectively. The lower frames include histograms of the τ/τ0 values with their mean and standard deviation values.

4.7 Further improvements

All fits in Sect. 3 were based on SED models that consist of discrete temperature components. This is not a fundamental limitation, and the basis functions could equally be based on parameterised temperature distributions. Appendix D shows one example, where the SED templates correspond to different Gaussian distributions of the optical depth versus dust temperature. The model parameters are thus the mean value 〈T〉 and the width δΤ of the temperature distribution and the optical depth. Although the run time was longer (about 1 h), the results were in this case also more accurate than those of the 3-TMPL model (Sect. C) that has the same number of free parameters.

The model used in Appendix D is different from both the N-TMPL and N-MBB ones, and it combines a low number of free parameters with the flexibility of adjusting the mean temperature in a continuous fashion. The MCMC runs used a lookup table of pre-calculated dust spectra. The same approach could be used for any SEDs, which could be based on empirical (e.g. log-normal or broken power laws) or model-based temperature distributions (e.g. Dale et al. 2001; Galliano et al. 2011; Galliano 2018). For the former, Désert (2022) advocated the use the logarithm of temperature as the more natural parameter. The basis functions could also directly encode prior information, such as the avoidance of unphysically low temperatures.

In this paper we have not considered the possible correlations between the input parameters. If the measurements of the different bands are correlated, the χ² values are calculated by replacing the sum of ${\left(X_i-{\widehat{X}}_i\right)}^2/\sigma {\left(X_i\right)}^2$ values with $\left(X-\widehat{X}\right){\text{Σ}}^{-1}{\left(X-\widehat{X}\right)}^T$, where X are the measured values, X the model prediction, and Σ the covariance matrix of the observations. With a handful of frequency bands, the effect on the run times would be at most a factor of a few. The spatial correlations of the input data have no effect on the 1-MBB fits, which are done independently for each resolution element. The correlations should nevertheless be taken into account when the data are initially convolved to a common resolution. If the convolution is part of the fitting procedure, the handling of the spatial correlations would have a large impact on the run times. Depending on the beam sizes, each pixel might have non-zero covariances with tens or even hundreds of other pixels. The additional cost would come from the increased complexity in the convolution (to be done now in real space) and in the evaluation of the goodness-of-fit function. The implementation would be more straightforward in the MCMC framework, where the derivatives of the probability function are not needed. However, there would still be significant practical limitations on the size of the maps that could be analysed.

5 Conclusions

We have examined the analysis of dust emission spectra with SED models that consist of N-fixed SED templates (models N-TMPL, N > 1) or N MBB functions (models N-MBB, N ≥ 1). Tests with synthetic observations have led to the following conclusions.

• The N-TMPL and N-MBB models can improve the standard single-component 1-MBB fits, but they require a careful model setup and, in the case of large N, the use of appropriate priors.

• The accuracy of the N-TMPL models critically depends on the selection of its SED templates. Errors increase rapidly when the temperature of the emission falls outside the temperature range covered by the templates.

• The N-MBB models adjust better to large temperature variations in the source maps, but they may need to be combined with priors to avoid very low temperatures that would cause large noise in the τ estimates. This is true especially if there are many pixels per beam and if the SED of individual pixels is thus not well constrained.

• The fitting of multi-component models is computationally feasible even for large maps, provided that the gradients of the χ² function are calculated analytically. In our tests, SED fits were performed for maps of up to 1875 × 1875 pixels.

• When the beam convolution is part of every optimisation step, the use of GPUs can significantly reduce the run times.

• We have presented alternative N-TMPL fits that were based on the optimisation in the Fourier space. This results in significant speed-up in the computations, but with some limitations in the use of priors and general error estimates.

• MCMC remains feasible for maps of at least up to ~500 × 500 pixels (run times of approximately 1 h or less), even when each MCMC step includes the convolution of the model-predicted maps. MCMC provides better estimates of the statistical errors, but these do not give any indication of the often larger systematic errors.

• Multi-component fits can contain fundamental degeneracies, where (in the studied cases) the same observed intensities could correspond to a factor of two differences in the true optical depths.

• The bias of the 1-MBB fits can be decreased based on observed data alone by analysing the widening of the SEDs. Neural networks provide a general way to derive such corrections.

The results call for caution when interpreting the results from multi-component SED fits. This is particularly true if the number of free parameters exceeds the number of direct observational constraints and the results thus clearly depend on the chosen priors. Having low χ² values for the fits does not mean that the corresponding τ estimates would be correct, but the comparison of alternative SED models can give some indication of the magnitude of the model errors. Our initial results suggest that models with parameterised continuous temperature distributions might perform better than the models (such as N-TMPL and N-MBB) where only discrete temperatures are used. This needs to be investigated further in a wider range of test cases. Radiative transfer modelling also remains a possibility in the analysis of emission from simple sources, and it is not much more time consuming than the fitting of a complex multi-component SED model.

Acknowledgements

This work was supported by the Academy of Finland grant no. 348342.

Appendix A Fitting algorithms

When the beam-convolved model predictions are compared to observations, the fitting is accomplished by minimising the function $${\chi }^2={\displaystyle \sum _{i,k}\frac{S_{i,k}-{\text{Σ}}_{\text{c}}{M}_{i,\text{c}}\otimes B_k}{\delta S_{i,k}^2}}={\displaystyle \sum _{i,k}\frac{S_{i,k}-{\widehat{S}}_{i,k}}{\delta S_{i,k}^2}},$$(A.1)

where i refers to the map position and k to the frequency band, and the model is the sum of multiple components indexed with c. The symbol ‘⊗’ stands for the convolution between a model map M_i,c and the telescope beam B_k, and Ŝ_k denotes the convolved model predictions. The summation goes over all map positions i and bands k.

The functions M_i,c are sums of SED components. The N-MBB fits use sums of MBB functions $$M_{i,c,k}=MBB{\left(I_{i,\text{c}},T_{i,\text{c}},{\beta }_{i,\text{c}}\right)}_{i,k},$$(A.2)

where I_i,c is the intensity at a reference wavelength (see below), and the result is also a function of the band k. In this paper, β was kept constant and only I_i,c and T_i,c are considered free parameters. The N-TMPL model is somewhat simpler, $$M_{i,c,k}=W_{i,\text{c}}{C}_{c,k},$$(A.3)

with free parameters W_i,c,k and the fixed SED templates C_c,k.

When the χ² functions is minimised with the conjugate-gradient method, only the function values and the first derivatives with respect to the optimised parameters X_i,c are needed. The derivatives with respect to $X_{i_0,j_0}$ can be written as $$\frac{\partial {\chi }_k^2}{\partial X_{i_0,c_0}}=-2{\displaystyle \sum _k\left[{\left(\frac{S_{i,k}-{\widehat{S}}_{i,k}}{d{S}_{i,k}}\right)}^2\otimes B_{i,k}\right]\frac{\partial M_{i,c,k}}{\partial X_{i_0,c_0}}}.$$(A.4)

Since the model prediction Ŝ itself is already the result of convolution, the gradient calculation thus includes two convolution steps. This can be shown by direct derivation of Eq. (A.1), by replacing the convolution there with an explicit summation over the beam.

The last term in Eq. (A.4) is the derivative of the model with respect to one of its parameter. In the sum, this is reduced to a single term (pixel i₀ and component c₀). If M is the MBB function, it may be numerically benefcial to express it using intensities relative to the value at a suitable reference frequency v₀, $$M{\left(I,T,\beta \right)}_{i,c,k}=I_{i,\text{c}}\frac{\exp \left(h{v}_0/\left(k{T}_{i,\text{c}}\right)\right)-1}{\exp \left(hv/\left(k{T}_{i,c}\right)\right)-1},$$(A.5)

where frequency ν₀ is the reference frequency and ν the frequency of the band k. The full M_i,k function is the sum of the individual components M_i,c,k. Using the values of the functions M_i,c itself, the final partial derivatives in Eq. (A.4) can be written as $$\frac{\partial M_{i,c,k}}{\partial I_{i,\text{c}}}=M_{i,k}/I_{i,\text{c}},$$(A.6) $$\frac{\partial M_{i,c,k}}{\partial T_{i,\text{c}}}=\frac{M_{i,k}hv}{T_{i,\text{c}}^2}\left(\frac{v}{1-\exp \left(-hv/\left(k{T}_{i,\text{c}}\right)\right)}-\frac{v_0}{1-\exp \left(-h{v}_0/\left(k{T}_{i,\text{c}}\right)\right)}\right)$$(A.7)

and $$\frac{\partial M_{i,c,k}}{\partial {\beta }_{i,\text{c}}}=M_{i,c,k}\log \left(v/v_0\right),$$(A.8)

thus replacing $X_{i_0,c_0}$ with I, T, or β. In the case of the N-TMPL flts, the derivatives ∂M/∂X are directly elements of the SED templates $$\frac{\partial M_{i,c,k}}{\partial W_{i,\text{c}}}=C_{i,k}.$$(A.9)

The non-negativity of the parameters I could be enforced with an additional penalty function or using an optimisation routine that directly support the use of constraints. We chose to optimise the parameter α = log(I), which also guarantees the positivity of the intensity values, I = e^α. For N-MBB flts, Eq. (A.6) becomes in this case $$\frac{\partial M_{i,\text{c}}}{\partial {\alpha }_{i,\text{c}}}=M_{i,\text{c}},$$(A.10)

where M itself is of course evaluated by replacing I with exp(α). In the N-TMPL flts, all ∂M/∂X derivatives are of the form $$\frac{\partial M_{i,c,k}}{\partial W_{i,\text{c}}}=e^{W_{i,\text{c}}}{C}_{c,k},$$(A.11)

if the original weights W_i,c are replaced with terms $e^{W_{i,\text{c}}}$.

When the optimisation of the N-TMPL model is done in Fourier space (models N-TMPL-F), observations S_i,k and the weight maps W_i,c are replaced in Eq. (A.1) by their Fourier counterparts and the convolution is replaced by direct multiplication with the Fourier-transformed beams. This is possible because the N-TMPL model is linear with respect to SED templates and $$S_k\approx {\displaystyle \sum _{\text{c}}{W}_{\text{c}}}{C}_{c,k}\Rightarrow ℱ\left(S_k\right)\approx {\displaystyle \sum _{\text{c}}{C}_{c,k}}ℱ\left(W_{\text{c}}\right),$$(A.12)

as used in Eq. (8).

The relative run times depend on many factors that include the map size, the number of fitted components, having convolutions as part of the model, and the computational methods used for the convolutions (i.e. either in real space or via FFTs). Figure A.1 compares the run time of the 1-MBB fit to 2-MBB, 4-TMPL, and 12-TMPL-F models for maps of different sizes. The observations are for the IRDC model and consist of 160, 250, 350, and 500 μm maps. The 1-MBB fit is by far the fastest, because the optimisation is done pixel by pixel and model does not include beam convolution. In this case, all maps are frst convolved to the same resolution before the SED fitting. The quoted run times are for the optimisation only, and these remain below one second even for the largest maps (when performed on a GPU). The other methods are slower because of the direct additional cost of the convolutions (although that cost is small in case of TMPL-F) and because the convolution connects all parameters and makes the model ftting a more diffcult global optimisation problem. The 12-TMPL-F run times correspond to a full run to a converged solution, where all calculations (including the conjugate gradient optimiser) are run on a GPU. In the case of the 2-MBB and 4-TMPL flts, the run times correspond to 13000 function evaluations. This may be sufficient in some cases but a good convergence may require even 2-3 times longer runs. In the 2-MBB and 4-TMPL fts, the optimiser is run on the host machine but the FFT convolutions are performed on a GPU.

As mentioned in Sect. 4.6, FFT calculations were faster on a GPU only for map sizes larger than 256×256 pixels. Because the beam sizes correspond to less than ten pixels (FWΗΜ=36″ at 500 μm, compared to 6″ pixels), the convolutions could also be done with a GPU in real space, without a drastic increase in the run times¹. The cost of the real-space convolution also increases only directly proportionally to the map dimensions (∝ Ν × M) instead of the NM lοg(ΜΝ) scaling of 2D FFTs.

Fig A.1 Comparison of run times of different SED fits. The observations are from the IRDC model, and the map size is varied from 128 × 128 to 1875 × 1875 pixels. Apart from 1-MBB, beam convolution is always part of model that is being fitted. The number of fitted SED components is given in the legend, and TMPL – F stands for optimisation in the Fourier space.

The calculations for Fig. A.1 were performed on a laptop with a discrete GPU (which is nevertheless three generations behind the most recent hardware available). When the optimiser itself was run on the host computer, the run times were mostly bound by the host (the single-threaded program) rather than the FFT calculations on the GPU.

Appendix B Additional fits of the full IRDC maps

Figure B.1 shows 11-TMPL and 4-MBB fits for the whole 1875×1875 pixel maps of the IRDC cloud model. The fit parameters are the same as in Fig. 9 in Sect. 3.4.2, where the smaller 256×256 pixel maps correspond to the upper left corner of these full maps.

Section 3.4.3 discussed fits of the 4-TMPL and 12-TMPL models with optimisation in the Fourier space. Figure B.2 shows fits for the full 1875× 1875 pixel maps using the same parameters.

Fig. B.1 Fits of the 11-TMPL and 4-MBB SED models to full IRDC cloud maps (1875×1875 pixels). The parameters are the same as in Fig. 9.

Fig. B.2 Comparison of 1-MBB fits and 4-TMPL-F and 12-TMPL-F fits performed in Fourier space. The parameters of the fits are the same as in Fig. 9, which corresponded to the upper left corner of the full IRDC maps.

Appendix C MCMC calculations and estimated source sizes

Although MCMC calculations are expected to be more time consuming than direct ML calculations, they turned out to be feasible for maps of moderate size. In the following, we show one example where the fit of the 3 – TMPL model is done with a basic MCMC simulation. The results are used to examine the estimated source sizes in the case of sources with radial temperature gradients.

The synthetic observations consist of surface brightness maps at 100, 160, 250, 350, and 500 μm, with angular resolutions of 3, 3, 4.3, and 6.1 pixels, respectively. The background has a temperature of T₁=15 Κ and an optical depth of one unit. We add to the maps Gaussian sources that have a size of FWΗΜ=20 pixels and a peak optical depth of two units. Each source has a linear temperature gradient from 15 Κ at the distance of r = FWHM to a lower or a higher temperature T₂ at the source centre. The temperature T₂ of the individual sources is varied between 11 Κ and 19 K. We add to the data white noise that corresponds to 1% relative errors in the surface brightness measurements.

The data were fitted with the 1-MBB model (resolution 6.1 pixels) and with the 3-TMPL model (nominal resolution 3 pixels), using components with T_C=12.5, 15.0, and 17.5 K. No priors were used, except the requirement of non-negative weights in the 3-TMPL fit. The MCMC calculations of the 3-TMPL model included 200 000 MCMC steps, which was enough for reproducible results.

Figure C.1 shows a simulation of 320 × 320 pixel maps that contain 16 Gaussian sources. In the 1-MBB estimates, some residuals are visible for most sources (Fig. C.1c). The 3-TMPL fits have a higher resolution but clear errors are visible only a couple of the coldest and warmest sources (Fig. C.1d,f).

We fitted each of the 16 sources with Gaussians to study the errors in the estimated source sizes. To make the results comparable between the two fits, the 3-TMPL predictions were first convolved down to the resolution of the 1-MBB map. Figure C.2 shows the ratios of the estimated and expected beam-convolved FWHM sizes for the 16 sources in Fig. C.1. The x-axis corresponds to the temperature T₂ at the centre of each source. The 1-MBB fits are seen to overestimate the size of the coldest source by ~15% and underestimate the size of the warmest source by 10%. In the 3-TMPL fit the errors remain below ~5%.

We calculated for each source also mean optical depths over and area of 40×40 pixel (twice the source FWHM). Unsurprisingly, the bias of the 1-MBB estimates increases with the temperature difference |T₁ − T₂|, and the optical depth is underestimated by up to 15% for both the coldest and warmest sources. The 3-TMPL estimates show a positive bias of ~4%, with values dropping when T₂ goes outside the range of T_C values. Because of noise (and requirement of non-negative component weights) the 3-TMPL model is thus generally slightly overestimating the width of the temperature distribution. The results are on average more accurate than for the 1-MBB model. This is partly due to a better description of the temperature variations and partly due to the higher effective resolution, which reduces the mixing of temperatures within the beam.

Fig. C.1 Fits of 1-MBB and 3-TMPL models to maps containing 16 Gaussian sources with different radial temperature gradients. The upper frames show the true value τ0 (frame a), the 1-MBB estimates (frame b), and the ratio of the 1-MBB estimates divided by the true values (frame c). Data in frames b-c have a resolution of 6.1 pixels. The lower frames show results of the 3-TMPL fit: the fit residuals (frame d), the estimated τ values (frame e), and the ratio of the 3-TMPL estimates and the true values (frame f). The nominal resolution in the lower maps is equal to 3 pixels.

Fig. C.2 Test with sources with radial temperature gradients. The average optical depth estimates (circles) and the source FWHM estimates (squares) are plotted relative to their correct values. The symbols correspond to 16 sources that are plotted against their central temperature T2. The 1-MBB and the 3-TMPL results are plotted in black and red. respectively. The shading shows the MCMC error estimates for the τ/τ0 ratio in the 3-TMPL fits.

Figure C.2 also shows the MCMC error estimates for the optical depths. The shaded region corresponds to the interquartile range in the MCMC samples over the 40×40 pixels area around each source. In the plot, the raw error estimates (per resolution element) are scaled down by the ratio between 40 pixels (the averaged area) and nominal resolution of 3 pixels. The MCMC error estimates seem to be consistent with the scatter between the sources (of similar T₂) but of course does not cover the larger component of systematic errors.

For the above 3-TMPL fit of the 320×320 pixel map, the MCMC run took about 10 min and is thus not significantly slower than the ML solution. In our tests, the corresponding analysis of 576×576 pixel maps took about half an hour, in agreement with linear dependence on the number of pixels. However, all calculations were done using a GPU, and the run time on a CPU would be at least one order of magnitude longer.

Appendix D MCMC calculations with an alternative SED model

In this paper the N-TMPL and N-MBB models were based on a set of a few discrete temperatures. However, each template spectrum of the N-TMPL models and even the basis functions of the N-MBB models could equally be integrals over some temperature distributions.

As a basic example, we tested a three-parameter model where the parameters are the optical depth, the central temperature T_C and the width δΤ of a Gaussian temperature distribution. A grid of SED templates was first calculated as the function of T_C and δΤ. The fitted parameters were the indices to the SED grid and the scaling of the selected template spectrum. As in Sect. C, the model has three free parameters and the fit was performed with MCMC methods. We refer to this as the GRID model.

Figure D.1 shows the results for the same synthetic observations as in Sect. C (maps of 320×320 pixels). Because these consist of two discrete temperatures (with spatially varying T₂), the fitted model of a single Gaussian temperature distribution does not exactly match the simulated data. The results are improved compared to the 3-TMPL fits in Fig. C.2, and the τ values show now very little bias. On the other hand, the run time (on GPUs) is also longer, as the calculation of three MBB functions is replaced by the more time-consuming look-up from a SED table. After an initial burn-in phase, the calculation of 200 000 MCMC steps took one hour.

Fig. D.1 Alternative fit of sources with radial temperature gradients. The symbols are as in Fig. C.2, but the red colour and the shaded error regions are for the GRID model.

MCMC calculations provide complete information of the posterior probability of the fitted parameters: the optical depth τ, the mean temperature 〈T〉, and the width of the temperature distribution δΤ. Figure D.2 shows examples of these distributions. With the high S/N of the observations assumed in this example, there are no strong correlations between the model parameters.

Figure D.3 compares the power spectra of the different optical depth maps. The analysed region excludes the borders where τ(GRΓD) solution has not been calculated (cf. Fig. C.1f). When the GRID map is convolved to the resolution of the 1-MBB map, the power spectra are nearly identical up to k ~ 0.05 pixel⁻¹. We assumed here that, before this convolution, the GRΓD model would have the nominal resolution of 3 pixels compared to the 6.1 pixel resolution of the 1-MBB model. Figure D.3 is approximately consistent with this assumption. The 1-MBB and GRID power spectra separate only at high spatial frequencies, presumably because the GRID solution has a higher statistical noise, even when degraded to the resolution of the 1-MBB solution.

Fig. D.2 Corner plot showing parameter correlations for three areas in the fits of Fig. D.1. The parameters are the optical depth τ, the mean temperature TC, and the standard deviation of the temperature distribution δΤC. The data correspond to averages over 3 × 3 pixels towards the centre of the coldest source (blue), a position in the 15 Κ background (black), and the centre of the warmest source (red colour). The contours are drawn at 10%, 50%, and 90% of the peak probability density.

Fig. D.3 Power spectra of the optical depth maps for the fits in Fig. D.1. The spectra are shown for the true optical depth τ0 and the 1-MBB and GRID estimates. For the comparison, the other maps are convolved to the resolution of the 1-MBB map.

References

All Figures

Fig. 1 Ratio of estimated and true optical depths τ/τ0 for 1-MBB, 2-TMPL (TC=13 K and 17 K), and 2-MBB models. The observations consist of 160, 250, 350, and 500 μm intensities with 3% relative noise. The dots correspond to individual pixels (every 10th pixel), the solid lines show moving averages, and the shaded regions the 1-σ dispersion of the τ estimates. The χ1 values of the ML fits (normalised by the number of fitted data points) are quoted in the frames. In frame a, observations consist of one emission component at T1 = 15 K and one component varied in the range T2 = 10–20 K (on x-axis), both with the same optical depth. The additional solid blue lines correspond to the 1-σ error estimates from the MCMC fit of the 2-TMPL model. In frame b, the observations correspond to a continuous temperature distribution, T ~ N(〈T〉, σT). The x-axis shows the 〈T〉 value, with σT = 1 K for the first and σT = 5 K for the second half. Frame c is the same as frame b but the SED fits assuming β = 1.7 instead of the correct value of β = 1.9.

In the text

	Fig. 2 Same as Fig. 1, but forcing each component of the 2-MBB and 2-TMPL models to positive optical depths. In the case of 2-MBB. the model also includes a prior on the temperatures, Τ ~ Ν(18 Κ, 5 Κ). The plotted 2-TMPL and 2-MBB data are from MCMC calculations while the χ2 values are the average values of the maximum-likelihood solutions.
In the text

	Fig. 3 Same as Fig. 2b, but including results for 3-TMPL (MCMC fit) and 3-MBB (ML fit) models. The 3-TMPL model uses TC = 13, 15.5, and 19 K, and the 3-MBB model includes a temperature prior Τ ~ N(15 K, 3 K).
In the text

Fig. 4 SED fits considering the effect of telescope beams. The 2-TMPL and 2-MBB cases include beam convolution as part of the fit while the 1-MBB fit is done at the lowest common resolution. The dust emission consists of the two temperature components shown in frames a and b. The map size is 128 × 128 pixels, the pixel size 10″, and FWHM beam sizes 20″, 20″, 25″, and 36″ at 160, 250, 350, and 500 μm, respectively.

In the text

Fig. 5 Fits on synthetic observations of an externally heated Bonnor-Ebert sphere. Frames a–c show the mass-weighted average LOS temperature of the model cloud and the calculated 160 μm and 500 μm surface brightness maps. Frames d-e show the temperature maps and frame f shows the 500 μm predictions from the 2-MBB model. Frame g shows the ratio of the estimated and true optical depth (τ/τ0) as the function of τ0 for the four cases listed in the legend. The solid lines are moving averages, the shaded regions indicate the 1-σ interval in the τ/τ0 values, and the individual pixels outside this interval are shown as dots.

In the text

	Fig. 6 Same as Fig. 5, but for a BE sphere with an additional internal radiation source. The TC values of the 2-TMPL and 3-TMPL models are given in the text. The resolution is 36″ for the 1-MBB model and 18″ for the other fits (for a nominal pixel size of 6″).
In the text

Fig. 7 SED fits of the filament models. Frame a shows the SEDs towards the map centre and the 1-MBB fits to the 4–6 longest wavelength bands. The lower set of symbols and curves corresponds to a model without a point source, and the upper ones to a model with a point source on the LOS towards the filament but at 0.93 pc in front of the filament. The colour temperatures of the 1-MBB fits are quoted in the frame. Frame b shows the ratio of recovered and true peak optical depths (left axis, filled symbols) and the average χ2 per data point over the whole map (right axis, open symbols). The x-axis indicates the four models. The square and star symbols are, respectively, for the externally heated model and the model with a point source. The black, blue, and red symbols correspond to fits that extend to shortest wavelengths of 160, 100, and 70 μm, respectively. The error bars show the standard deviation of the optical depth estimates along the filament spine.

In the text

Fig. 8 IRDC observations fitted with 1-MBB, 2-MBB, and 4-TMPL models. Frames a-b show the 250 μm intensities and colour temperatures from the 1-MBB fit. The second row shows the 1-MBB τ estimates relative to the correct values (frame c) and the χ2 values (frame d). The third and fourth rows show the corresponding data for the 2-MBB and 4-TMPL flts. The plots cover a 256 × 256 pixel area of the IRDC cloud model.

In the text

	Fig. 9 Fits of IRDC model maps with 4-MBB and 11-TMPL models. The corresponding data for the 1-MBB fit are shown in Fig. 8.
In the text

Fig. 10 IRDC observations fitted in Fourier space. Frame a shows the true optical depth and frame e the colour temperature from the 1-MBB fit. Frames b–d show the optical depths estimated with the 1-MBB, 4-TMPL, and 12-TMPL fits, respectively. The corresponding ratios of estimated and true optical depths are shown in frames f-h and the χ2 values in frames j-l. Frame i shows the standard deviation of the temperature values from the 12-TMPL fit. The values quoted in the frames are the 10%, 50%, and 90% percentiles of the plotted quantity. The maps cover the same 256 × 256 pixel area as in Fig. 8.

In the text

	Fig. 11 Simulated single-pixel observation. Frame a shows the optical depth distribution as a function of temperature, and frame b shows the SED, with measurements at 160, 250, 350, and 500 μm.
In the text

Fig. 12 Results from the fitting different noise realisations of the Fig. 11 observations. Frame a shows a few examples of the observations (black crosses) and individual fits (solid curves). Frame b–d show the distributions of τ estimates (100 noise realisations) for the 1-MBB, 2-MBB, and 3-MBB models, respectively. Each frame quotes the mean value and the standard deviation for the estimates relative to the true value, τ/τ0.

In the text

	Fig. 13 Values of χ2 (frame a) and τ/τ0 (frame b) as functions of the assumed temperatures T1 and T2 of the fitted SED templates. The observations correspond to one noise realisation selected from the simulations in Fig. 12. Frame b includes contours at χ2 = 0.1 and χ2 = 1.
In the text

Fig. 14 Fits to simulated observations with a 10-TMPL model. The observations consist of 160, 250, 350, and 500 μm intensities for a single Τ = 15 Κ MBB function (upper frames) or a sum of 13 Κ and 17 Κ components of equal optical depth (red circles in frames b and e). Frames a and d show random realisations of the combinations of weights W that result in less than 5% errors for all four predicted intensities (curves in frames b and e). The frames c and f show the distributions of the τ values predicted by those fits. The true optical depth is τ(250 μm) = 0.01.

In the text

Fig. 15 IRDC model observations fitted with MCMC methods. Frame a shows the true τ(250 μm) optical depths, and frame d the colour temperature from the 1-MBB fit. Frames b and c show the τ maps from the 4-TMPL and 2-MBB fits (averages over the MCMC samples) divided by the true optical depth τ0. Frames e and f show the posterior probability distributions for the four positions marked in the b and c frames, and the true values τ0 are indicated with vertical lines of the same colour (black, blue, red, and cyan).

In the text

Fig. 16 Empirical correction of the 1-MBB optical depth estimates of the IRDC cloud model. Frame a shows the ratio τ/τ0 between the original estimates and the true optical depths. Frame b shows the quantity Σ(ΔIv/Iv) that characterises the difference between each observed spectrum and the 1-MBB fit. Frame c shows the correlation between this quantity and the τ error, and Frame d shows the τ/τ0 ratio corrected for the trend in frame c. In frames a and d, the inserts show histograms of the τ/τ0 values (with linear y-axis). In frame c, the colour scale corresponds to logarithmic point density.

In the text

Fig. 17 Correction to 1-MBB τ estimates calculated with neural nets. Frame a shows the 246 × 246 pixel map that is used as the training set. Frame c shows the 1-MBB results for the whole IRDC map, where the location of the training map is indicated with a dashed box in the upper left corner. Frames b and d shows the corrected τ maps for the training set and for the full 1875 × 1875 pixel map, respectively. The lower frames include histograms of the τ/τ0 values with their mean and standard deviation values.

In the text

	Fig A.1 Comparison of run times of different SED fits. The observations are from the IRDC model, and the map size is varied from 128 × 128 to 1875 × 1875 pixels. Apart from 1-MBB, beam convolution is always part of model that is being fitted. The number of fitted SED components is given in the legend, and TMPL – F stands for optimisation in the Fourier space.
In the text

	Fig. B.1 Fits of the 11-TMPL and 4-MBB SED models to full IRDC cloud maps (1875×1875 pixels). The parameters are the same as in Fig. 9.
In the text

	Fig. B.2 Comparison of 1-MBB fits and 4-TMPL-F and 12-TMPL-F fits performed in Fourier space. The parameters of the fits are the same as in Fig. 9, which corresponded to the upper left corner of the full IRDC maps.
In the text

Fig. C.1 Fits of 1-MBB and 3-TMPL models to maps containing 16 Gaussian sources with different radial temperature gradients. The upper frames show the true value τ0 (frame a), the 1-MBB estimates (frame b), and the ratio of the 1-MBB estimates divided by the true values (frame c). Data in frames b-c have a resolution of 6.1 pixels. The lower frames show results of the 3-TMPL fit: the fit residuals (frame d), the estimated τ values (frame e), and the ratio of the 3-TMPL estimates and the true values (frame f). The nominal resolution in the lower maps is equal to 3 pixels.

In the text

Fig. C.2 Test with sources with radial temperature gradients. The average optical depth estimates (circles) and the source FWHM estimates (squares) are plotted relative to their correct values. The symbols correspond to 16 sources that are plotted against their central temperature T2. The 1-MBB and the 3-TMPL results are plotted in black and red. respectively. The shading shows the MCMC error estimates for the τ/τ0 ratio in the 3-TMPL fits.

In the text

	Fig. D.1 Alternative fit of sources with radial temperature gradients. The symbols are as in Fig. C.2, but the red colour and the shaded error regions are for the GRID model.
In the text

Fig. D.2 Corner plot showing parameter correlations for three areas in the fits of Fig. D.1. The parameters are the optical depth τ, the mean temperature TC, and the standard deviation of the temperature distribution δΤC. The data correspond to averages over 3 × 3 pixels towards the centre of the coldest source (blue), a position in the 15 Κ background (black), and the centre of the warmest source (red colour). The contours are drawn at 10%, 50%, and 90% of the peak probability density.

In the text

	Fig. D.3 Power spectra of the optical depth maps for the fits in Fig. D.1. The spectra are shown for the true optical depth τ0 and the 1-MBB and GRID estimates. For the comparison, the other maps are convolved to the resolution of the 1-MBB map.
In the text