EDP Sciences
Free Access
Issue
A&A
Volume 604, August 2017
Article Number A58
Number of page(s) 17
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201629944
Published online 04 August 2017

© ESO, 2017

1. Introduction

Dust chemistry plays an important role during the evolution of interstellar clouds. The presence of dust is ubiquitous in the interstellar medium (ISM) and it is an important constituent of the Galaxy. Having a mass of only about 0.7% of the gas (Fisher et al. 2014), these microscopic particles greatly impact the chemistry and thermodynamics of gaseous clouds (Gerola & Glassgold 1978; Dopcke et al. 2013; Hocuk et al. 2014), along with dominating the continuum opacity. Gas-phase species can use grain surfaces as a third body to form more complex molecules, thereby catalyzing reactions which may otherwise be too slow to be significant (e.g., Gould & Salpeter 1963; Cazaux et al. 2005; Garrod et al. 2008; Gavilan et al. 2012; Ruaud et al. 2015). Atoms and molecules can, on the other hand, also be depleted from the gas phase when the dust temperature is too cold for species to overcome the thermal desorption energy (e.g., Jones & Williams 1984; Lippok et al. 2013). In this way, the dust temperature crucially controls whether gas-phase species freeze out onto dust, or are enriched from the chemistry occurring on dust grains (Tafalla et al. 2002; Garrod & Herbst 2006; Hocuk et al. 2016).

The dust temperature is also an important parameter in chemical reaction rates. An exponential dependence on the dust temperature lies at the heart of most surface reactions. At cold 10 K temperatures, for example, the difference of a single Kelvin can imply a variation of the reaction rates by orders of magnitude. Thus, a precise knowledge of the dust temperature is imperative when performing rate calculations. Fortunately, calculated abundances seem to be more sensitive to the relative reaction rates between those which compete with each other rather than the absolute reaction rates (see e.g., Tielens & Hagen 1982; Chang et al. 2007; Cazaux et al. 2016). Nonetheless, a study on the dependencies of the dust temperature is highly desirable.

A third and fundamental importance of the dust temperature follows from its impact on the gas temperature. The gas-dust collisional coupling is the single most important heat transfer mechanism for gas at number densities above a few ×104 cm-3 for typical Galactic conditions (Hollenbach et al. 1991; Spaans & Silk 2006). This process dominates the gas cooling as long as the dust temperature is lower than the gas temperature. At densities of roughly >106 cm-3, the temperatures of the two phases are irrevocably linked, such that the gas temperature is essentially set by the dust temperature. In regions with density below ~104.5 cm-3, where line emission regulates the gas temperature, the temperature of the dust still plays a role. Here, dust can influence the gas temperature because the molecular abundances of species such as CO and H2 (at T> 100 K), which depend on surface chemistry, control the amount of ro-vibrational line cooling (see e.g., Hocuk et al. 2016).

In this work, we have derived the dust temperature semi-analytically for various types of grain material through solving the energy balance. We compare our results to a collection of observed dust temperatures with the Herschel Space Observatory. In this way, by constraining our calculations with the observed dust temperatures, our study sheds light on the composition of dust in the ISM. The dust temperature is also explored for the presence of ices on dust surfaces. To substantiate our method, we test our semi-analytical solutions against a numerical one, i.e., radiative transfer calculations.

In Sect. 2, we report on a collection of observations of the interstellar dust temperature from the literature. In Sect. 3, we describe the analytical method that we use in order to derive the dust temperature, where we discuss the parameter details in Sect. 4. We present our semi-analytical solutions for the dust temperature for various grain materials as well as for ice coated grains in Sect. 5. In Sect. 6, we compare our semi-analytical solutions against computations with the Monte Carlo radiative transfer code radmc-3d1 (Dullemond 2012). We then introduce our parametric expression for the dust temperature in Sect. 7. In Sect. 8, we discuss our theoretical solutions with respect to the observed dust temperatures. Finally, in Sect. 9, we summarize our conclusions.

2. Observed dust temperatures

We first looked at observations reported in the literature to check if a simple expression for the dust temperature as a function of AV can be found. We show that this is not practical because one is probing different environments and different conditions among the various sources included and because of a missing general understanding of the physics. We will use these observations at a later stage to derive our parametric formalism from the inferred dust properties (in Sect. 8).

2.1. Herschel observations

Recent observations with the Herschel Space Telescope provided a large number of measurements of the dust temperature in various sources. These sources consist of dense filaments, clumps, starless and prestellar cores, and protostellar cores. We select eight independent studies that report the dust temperature for a variety of sources and adopt their values, but exclude the protostellar cores from our selection as these have fundamentally different environments due to protostellar feedback. We present a compilation of the published dust temperatures at the specified column densities, NH, or AV. To stay within the same units, we convert NH to AV using the conversion factor 2.2 × 1021 cm-2 mag-1 (Güver & Özel 2009; Valencic & Smith 2015).

Figure 1 displays the collection of observationally obtained dust temperatures, with the references given in the legend of the figure.

thumbnail Fig. 1

Observed dust temperatures from eight independent studies. The dust temperature is plotted as a function of visual extinction. A least squares semi-log linear line is fit through the data as given by the dashed black line.

Open with DEXTER

We have drawn the best fit semi-log linear line through the data points, but the functional form is arbitrary. The fitting function we obtain is Td(AV) = 14.4 − 3.73log 10(AV) K. However, there is no physical basis for this behavior and such a fit cannot show the dependence in the radiation field. Furthermore, by using this fit, one is not able to extrapolate with certainty beyond the bounds of the observed range (AV ≃ 0.2 − 70 mag).

The error bars are considered wherever they are available, though, only for the dust temperatures. The error bars for NH or AV are often not reported and, if given, can have large uncertainties. For example, in the case of the 21 cold clumps study by Parikka et al. (2015), two methods for calculating N(H2) are given, that is, from dust continuum and molecular lines, which diverge greatly in some cases and thereby influence the results. We adopted the one that is recommended (dust continuum) by these authors. For the study of starless and prestellar cores in L1495 of Taurus by Marsh et al. (2014) we recovered the column densities by computing this ourselves using the provided number densities, radii, and the Plummer-like density profile.

From the 14 low-mass molecular cloud cores (EPoS project, Launhardt et al. 2013), we took the 7 starless cores and excluded the protostellar cores. Also from the bok globule CB244 (Stutz et al. 2010), we only considered the starless core measurement. The only filament in our collection is the dense filament of the Taurus molecular complex L1506 (Ysard et al. 2013), which appears to have a density of nH> 103 cm-3, where a 3D radiative transfer model is used for estimating the emission and extinction of the dense filament.

In the studies of the isolated starless core B68 (Nielbock et al. 2012), the star-forming core CB17 (Schmalzl et al. 2014), the dense cores in the L1495 cloud of the Taurus star-forming region (Marsh et al. 2014), and the starless cores B68 and L1689B (Roy et al. 2014) various techniques have been used to remove the line-of-sight (LOS) contamination, always resulting in a lower dust temperature (by about 0 − 4 K) than the ones obtained from dust spectral energy distributions (SED) only2. The used LOS correction techniques are, in the order of the above listed sources, an employed ray-tracing model, a modified black body technique together with a ray-tracing technique, a radiative transfer model (corefit/modust), and an inverse-Abel transform-based technique.

Despite such mixed origins in Fig. 1, i.e., with and without LOS corrections, the difference is not directly obvious from the plot, except that data points corrected for LOS effects generally have a lower temperature toward higher AV. This can be perceived by looking at the red and black points.

2.2. Environmental differences

The external radiation field strength in many of these sources is not known. In units of the Habing field (Habing 1968), the ones that are known have the best estimates of G0 = 1 for L1506 and B68 (Ysard et al. 2013; Roy et al. 2014), G0 = 0.18 to 1.18 for the dense cores in L1495 (Marsh et al. 2014, their standard χISRF corresponds to G0 = 1.31), and G0 ≈ 2 for L1689B (Steinacker et al. 2016a), which was initially (over)estimated to be G0 ≈ 10 (Roy et al. 2014). For the cores CB244, B68, and CB17, Lippok et al. (2016) estimate an enhancement of factor 2.5, 2.2, and 3.0, respectively, relative to their interstellar radiation field (ISRF). The thermal dust temperature in the diffuse medium surrounding the remaining objects ranges between ~16 to 20 K. These values are typical for the Milky Way diffuse ISM and the ISRF is thus assumed to be close to the standard Galactic radiation field (G0 ≈ 1 − 2). Other differences, such as, number density, dust-to-gas mass ratio, type and size of grains, turbulence, and magnetic field are all important factors not taken into account. For these reasons, we targeted similar types of environments, the starless cold dense regions, where the above mentioned conditions are not expected to vary greatly. In spite of the unknown intrinsic differences, we still proceeded to overlay the various observations in a single plot (i.e., Fig. 1) to give us a general indication about the dust temperature in such environments.

3. Solving for the dust temperature

A simple way of obtaining the dust temperature is by solving the dust thermal balance for equilibrium. The underlying assumption is that an equilibrium is quickly reached and maintained. The primary heating and cooling processes are fast enough, on the order of 10-5 s (radiative cooling) to minutes and hours (heating by interstellar photons, e.g., Draine 2003a), to justify this approach. Presuming that the main heating of dust is caused by the ISRF and that the primary cooling comes from the isotropic modified3 black body emission, the energy balance for a single grain can be set up as follows (cf. Krügel 2008) (1)Here, Qν is the absorption efficiency (Sect. 4.2) that depends on a, the grain radius, Td is the dust temperature, Bν is the Planck function, Jν is the ISRF flux (Sect. 4.1), and Dν(AV) is the attenuation factor (Sect. 4.3), with AV the visual extinction. In this work, we assume a geometry that is spherically symmetric and that the interstellar radiation is coming from all directions, but that the cloud is large enough to shield the radiation from one side. Hence, we take 2π for the right-hand side. Our choice is best described by a semi-infinite slab. This represents the edges of a cloud very well, whereas the center of a cloud should tend to 4π. Assuming that a medium is in radiative equilibrium, for a distribution of dust grains, one may apply a second integral over grain sizes in Eq. (1). The integral equation becomes independent of a if a fixed size is adopted. In this work, due to the complex and evolving nature of dust grains, we present our solutions for the canonical size of a = 0.1 μm (e.g., Kruegel 2003).

This relatively simple concept is often adopted to obtain the dust temperature. This usually involves assumptions for certain aspects of the calculation (i.e., Qν, Jν, or Dν) or is simplified by limiting the solution to a desired range, which we discuss in the next section. Equilibrium solutions are, of course, always time independent and tend to consider simple geometries, like a slab or a sphere. The benefit is that the calculation is fast and stable, while the solutions are considered to be satisfactory.

It is advisable to note that for small dust grains (a ≲ 50 Å), the equilibrium solution will not hold, since there will be large temperature fluctuations following single-photon heating events (Draine & Li 2001). Although the mass in small grains is low, a substantial fraction of the emission from diffuse clouds may be coming from them (Draine & Li 2001; Li & Draine 2001). Non-equilibrium solutions should therefore be used for a better treatment of very small grains in diffuse regions.

One can expand the energy balance equation by adding more heating and cooling terms. One factor that may be important for the heating of dust grains are cosmic rays (CR). Upon hitting a dust grain or a molecule, cosmic ray particles can either heat the grain locally, i.e., impulsive spot heating (Leger et al. 1985; Ivlev et al. 2015), or globally, for example, due to secondary UV photons generated following H2 fluorescence in the Lyman and Werner bands. CRs are insensitive to a gas column density of up to NH ≳ 1023 cm-2 (Padovani et al. 2009; Indriolo et al. 2015) and have an attenuation length of NH ≃ 6 × 1025 cm-2 (Umebayashi & Nakano 1981). This renders the energetic secondary UV photons nearly independent on the cloud optical depth, which results in a constant contribution to the right-hand side of the energy balance equation. For heating the dust grains by secondary UV photons, adopting a CR induced UV photon (CRUV) flux of FUV = 2 × 104 s-1 cm-2 (the “low” proton model of Ivlev et al. 2015) and a photon energy of 13 eV per CRUV, the total intensity becomes ICR = 4.2 × 10-7 erg s-1 cm-2 sr-1. From these estimations, we can already report that the CR impact on the dust temperature for a standard Milky Way ionization rate, i.e., ζH2 = 5 × 10-17 s-1, turns out to be minimal. We find that due to CRs the dust temperature increases by ΔTd ≲ 0.1 K. Higher CR rates have been reported (e.g., Indriolo et al. 2015) and some CR attenuation is expected (e.g., Ivlev et al. 2015). The impact of higher cosmic ray rates for the considered simplified case (i.e., without attenuation) is explored in Appendix A.

4. Parameter details

4.1. Jν, the ISRF

Each parameter in the integral of Eq. (1) is a function of frequency that goes from 0 to infinity. In actuality, however, this limit is set by the ISRF, which affects the other parameters. Typically, the Galactic ISRF covers the wavelength regime between microwave (3000 μm) and far-ultraviolet (FUV, 0.1 μm), with contributions from stars (OB stars and late spectral classes), dust (cold, warm, and hot), and the cosmic microwave background (CMB). Starlight dominates the emission between the FUV and the near-infrared (NIR), while dust mainly emits reprocessed starlight between the NIR and the far-infrared (FIR). Beyond this range, there is little emission to be significant for the dust temperature. The shorter wavelengths are more important toward lower AV, while the longer wavelengths become more relevant at higher optical depths.

thumbnail Fig. 2

Parameters for the dust equilibrium. Top panel, the ISRF intensity as a function of frequency. The green line represents the adopted ISRF in this work. Bottom left panel, experimental and calculated absorption efficiencies for various grain materials. Scattering is not included in these. Bottom right panel, extinction curves from various studies. The filled black circles show the observed data from Mathis (1990). The black solid line is a fit to the data as given by Cardelli et al. (1989), which is adopted in this work. The sub panel zooms in at the lower wavenumbers where the adopted extinction curve below 0.15 μm-1, given by the black solid line, is interpolated from Mathis (1990).

Open with DEXTER

Because the ISRF has contributions from various sources, it results in a continuum spectrum that goes from the FUV to the CMB. The shape and power of this spectrum is described by Mathis et al. (1983) and Black (1994), though, we do note that the Galactic ISRF of most cores is anisotropic and dependent on the location within the Galaxy. In the work of Zucconi et al. (2001), henceforth to be referred as ZWG01, the ISRF is approximated by a parametric fit to the data of the above mentioned authors. However, the UV part of the spectrum is omitted, because the aim of ZWG01 was to model the temperature at AV> 10 mag. This part of the spectrum can be covered by adopting the UV background of Draine (1978). This approach is also reported by Glover & Clark (2012) and Bate & Keto (2015). Aside from adding the UV part of the spectrum down to λ = 0.091μm (13.6 eV), we also adapt the mid-infrared (MIR) part of the spectrum to a smooth modified black body instead of the power-law with a cut-off at 100 μm of ZWG014. Our modified function at the MIR (around ν = 2 × 1013 s-1) is (2)where h is the Planck constant, kB is the Boltzmann constant, c is the light speed, and Wi is the weighting factor. Our fitted values give us Wi = 3.4 × 10-9 and Ti = 250 K. This approach approximates the data better than the mentioned partial power-law (see the top panel of Fig. 2), albeit that this part of the spectrum is actually not smooth and largely dominated by polycyclic aromatic hydrocarbon (PAH) emission, see for example, Porter et al. (2006). The full expression of the adopted ISRF is given in Eqs. (B.1) and (B.2).

In the top panel of Fig. 2, we plot the ISRF as mean intensity Jν (erg s-1 cm-2 sr-1 Hz-1) versus frequency. The black filled circles in this figure represent the original data from Black (1994) and Mathis et al. (1983), which have a low UV contribution, whereas the black long-dashed line displays the Draine (1978) UV field. The red short-dashed line shows the fit of ZWG01. We have also added in this figure the ISRF adopted by Hollenbach et al. (1991), from here on HTT91, by two separate functions: the FIR/CMB and the UV part are both shown on the figure. The green solid line, that is adopted in the current work, shows the combined ZWG01 and Draine ISRFs, which covers the whole frequency range, in consensus with earlier works (Glover & Clark 2012; Bate & Keto 2015).

4.2. Qν, the absorption efficiency

Table 1

Considered dust material types adopted from literature.

Interstellar dust grains efficiently absorb photons with wavelengths smaller than their own size. Longer wavelength radiation is not entirely absorbed and there is an efficiency related to this, which is given by the frequency dependent parameter Qν. Scattering is not considered in the efficiency Qν since scattered radiation will not thermally affect a dust grain. Scattering may extend the path length of a photon which increases the probability of absorption of radiation. For this, one ideally needs to keep track of all the scattered radiation at each point in a cloud. This can be done numerically. We consider scattering through the attenuation of radiation, which is discussed in Sect. 4.3. Scattering is very inefficient at wavelengths much larger that the size of the scattering object (λ-4), but may become important at wavelengths around λ ≲ 1 μm.

We take Qν = 1 for λ ≪ 2πa, i.e., the geometric optics approach, and Qνσext/πa2, which is less than unity, in the Rayleigh limit λ ≫ 2πa, where σext (cm2) is the extinction cross section. The dust emissivity for cooling is an equally important aspect as dust extinction is for heating. Both are treated by the same efficiency parameter Qν. It is generally acceptable to assume that the emission efficiency equals the absorption efficiency, i.e., Qν,em = Qν, since a good absorber is also a good emitter.

The absorption efficiency can be theoretically constructed, where the simpler models assume a power-law dependence, or can be experimentally measured, with material-specific absorption features. When the function follows a power-law, since Qν is on both sides of Eq. (1), only the slope of the function, i.e., the power of ν, matters for the dust temperature. In reality, however, dust has material-specific absorption features. With detailed semi-analytical calculations or direct laboratory measurements it is possible to obtain the opacity coefficient κν (cm2 g-1), also known as the mass absorption coefficient, with high precision. The relation between Qν and κν for spherical grains is as(3)where ρd is the bulk mass density of dust, which is roughly around 3 g cm-3 for silicate grains, but actually depends on the dust refractory material composition.

In the present work, we adopt a realistic set of absorption efficiencies. The opacities for the considered dust materials are gathered from the references provided in Table 1. The adopted absorption efficiencies are either experimentally measured or theoretically calculated from the optical properties of the refractory materials. The obtained data is in units of Qν (LD93), optical constants n,k (HM97, Dr03), or κν (OH94, WD01, KYJ15, Or11). The Mie theorem allows the calculation of κν from n,k (e.g., Bohren & Huffman 1983), while κν is converted to Qν using Eq. (3). This shows that we could just as well integrate Eq. (1) over a3κν (mass) instead of a2Qν (surface). Since some of the models have an underlying size distribution, the conversion with a fixed grain size makes the assumption that the canonical value of grain radius 0.1 μm represents the mass weighted average over the grain size distribution (Mathis et al. 1977; Kruegel 2003). We expect that the size of the grains will not significantly change during the evolution of a diffuse cloud to a prestellar core (e.g., Hirashita & Omukai 2009; Schnee et al. 2014; Chacón-Tanarro et al. 2017).

In the bottom left panel of Fig. 2, we show the absorption efficiencies obtained from detailed calculations and laboratory experiments for the various types of dust material. To have a matching wavelength coverage, we extrapolate the data where there are no measurements and we limit Qν ≤ 1, i.e., remain within the geometric optics approach (not allowing more than 100% absorption efficiency), following Hollenbach et al. (1991).

In our selected list of materials Graphite is calculated for 0.1 μm grains at 25 K, Silicate A (quartz glass) is measured at 10 K, Silicate B (astronomical silicate) is composed for 0.1 μm grains at 20 K, and Mixed D is calculated at 10 K. For Mixed A we adopt the uncoagulated model, for Mixed B we adopt the Milky Way RV = 5.5 model, where RV is discussed in the next section, and for Mixed C we adopt the uncoagulated “CMM” model. Other details can be found from the original papers.

4.3. Dν (AV), the attenuation factor

The attenuation of radiation in the ISM mostly arises from dust and large molecules. Depending on frequency, these particles absorb and scatter light. For radiation traveling through a medium, when calculating the optical depth τν considering only absorption, one has to integrate the absorption coefficient αν (cm-1) along the path s, that is, (4)αν is related to the opacity κν through the relation αν = ρκν. Since attenuation scales as exp( − τν), one needs to know the optical depth at all frequencies to find the solution for Td.

At visible wavelengths (λ = 5500 Å) the relationship between τV and AV is straightforward and the attenuation factor simplifies to exp( − 0.921 AV). Taking this as a reference, the attenuation at different wavelengths is scaled by the wavelength-dependent attenuation coefficient γλ (Aλ/AV), such that the attenuation factor Dν(AV) becomes (5)The coefficient γλ, that is given by the extinction law, is parameterized by Cardelli et al. (1989) for the Milky Way. The extinction law accounts for both absorption and scattering, and its shape is the result of the contribution of three main components: PAHs, small carbon grains, and silicates. We use their 5-part function for the wavelength range 0.1 μm to 3.4 μm. For longer wavelengths we adopt the tabulated values of Mathis (1990).

The attenuation coefficient is now only a function of the optical parameter RV, which is the total-to-selective extinction ratio RVAV/E(BV). Two typical RV’s in the ISM are for diffuse clouds with RV = 3.1 and for dense clouds with RV ~ 5. For the present work, we adopt RV = 5, but also discuss RV = 3.1 briefly. In Table 2, we show γλ for a number of wavelengths.

Table 2

Attenuation coefficient γλ (RV = 5).

The impact of RV on the extinction curve is quite large at shorter wavelengths, however, above λ ≃ 0.55 μm (or γλ< 1) the differences in γλ are small (<16%, for a nice overview see Fig. 2 of Mathis 1990).

In Fig. 2 bottom right panel, we show γλ as a function of wavenumber (1 /λ). Here, one can notice the broad band feature at 2175 Å and the FUV rise toward shorter wavelengths. The origin of the prominent 2175 Å feature is not fully understood, but is believed that carbonaceous materials, such as PAHs (Draine 2003a; Xiang et al. 2011; Zafar et al. 2012) or amorphous hydrocarbons (Jones et al. 2013), are responsible. We also overlay the extinction curves used by HTT91, ZWG01 (assuming a reference frequency of c/ 5500 Å), and Garrod & Pauly (2011), henceforth GP11 (assuming RV = 5). In the sub panel of Fig. 2 (bottom right), we zoom in on the lower wavenumbers to highlight the differences, which are important for embedded regions.

5. Dust temperature: semi-analytical solutions

5.1. Bare grains

We display the dust temperatures obtained by solving Eq. (1) with our semi-analytical model (described in Sect. 3) for the materials graphite, silicate (SiO2, MgFeSiO4), and carbonaceous-silicate mixtures in the left panel of Fig. 3. The adopted opacities for the seven different dust material types used in this work are provided in Table 1.

thumbnail Fig. 3

Dust temperature solutions. The left panel displays the obtained dust temperatures for various grain materials, which have no ices. The right panel shows the impact of ice formation.

Open with DEXTER

All model solutions arrive at a cloud edge (AV = 0.04 mag) temperature that lies between 14 − 18 K. The temperature variation at high optical depths (AV = 150 mag) is quite small for all dust types, ranging between 5 and 6 K, where the silicate dust types are generally warmer. The silicate dust types have a smaller temperature difference between the minimum and maximum AV and have different dependencies with extinction compared to the carbonaceous and mixed dust types. The graphitic dust temperature profile, on the other hand, is hard to differentiate from the mixed dust profiles. The larger differences that arise at lower AV are a result of the greater variation in the absorption efficiencies at higher frequencies, see Sect. 4.2. The similar temperatures for graphite and mixed dust and the lower silicate temperature at low AV is expected, since the opacity of carbon grains in the optical and near-infrared is higher than that of silicate grains. While for mixtures, the carbon component dominates.

5.2. Icy grains

We also examine the dust temperature when the grain surface is covered by (pure water) ice. Ice formation changes the opacity of dust in a way that it may reduce the dust temperature at the cloud edge, though, no ice is expected in these regions, while ice formation increases the dust temperature at high optical depths. This is because of the prominent ice features in the infrared (IR), which increase the opacity at those wavelengths. As heating effectively occurs throughout the spectrum, cooling of dust is, on the other hand, temperature dependent. Therefore, due to the ice bands in especially the near-IR (NIR) and the far-IR (FIR), the dust grain is either cooled or heated more depending on its temperature (see for a review Woitke 2015).

OH94, KYJ15, and Or11 have also modeled grain opacities by coating dust surfaces with ices. Their method utilizes the optical properties of water ice and by applying the effective medium theory (see e.g., Min et al. 2008; Woitke et al. 2016). The OH94 modeled icy mantles range from having a few monolayers5 of ice to 100 monolayers.

Using the OH94 ice models, we display in the right panel of Fig. 3 the dust temperatures for ice covered dust grains with thick ice mantles (ice volume 4.5 × grain volume, i.e., all water is frozen), thin ice mantles (ice volume 0.5 × grain volume), and no ice mantles. While not drawn in this panel, the KYJ15 and the Or11 ice models also indicate the same trend, having only thin ice and bare surface models to compare. We note that both the KYJ15 and the Or11 base ice model opacities are not entirely separable from coagulation, hence we did not display them in Fig. 3, but we review and show them in Appendix C. We thus can conclude that ice formation has a clear and notable impact on the dust temperature, and that at AV = 150 mag, the ices make a difference of about 0.8 K. Where thick ices result in Td = 6.1 K at AV = 150 mag, bare grains cool to Td = 5.3 K.

It is, however, not expected to have thick ice mantles on dust grains at low AV (3 mag), since UV radiation can photodesorb ices. Moreover, the adsorption rates will be quite low at the cloud edge because of the lower densities. In dense cores, on the other hand, the expectation is that thick ice mantles will cover the cold dust surfaces where UV radiation no longer plays a significant role. In a realistic case, one should go from the black line (bare grain) in Fig. 3, left panel, to the blue triple-dot dashed line (thick ice) transitioning around AV ~ 3 − 6 mag. We emphasize that this results in a less curved, quasi linear, thermal profile. This is an important point which should be taken into account in simulations.

6. Dust temperature: numerical solutions

In addition to the semi-analytical solutions, we calculate the dust temperatures using the same opacity models (Sect. 4.2) with the Monte Carlo radiative transfer code radmc-3d (Dullemond 2012). The numerical approach for solving the radiative transfer problem allows us to calculate the dust temperature in arbitrary geometries and density distributions. Despite its flexibility, it is not always suitable for large scale hydrodynamical simulations that require time-dependent solutions, because the Monte Carlo approach is computationally intensive.

There is a noteworthy fundamental difference between the semi-analytical and the numerical method. In the former case, the efficiencies adopted for the heating and cooling (through Qν in Eq. (1)) differ from the opacity adopted for the attenuation of the ISRF (i.e., Dν in Eq. (1)). Where we choose to use the observationally obtained attenuation factor provided by Cardelli et al. (1989) and Mathis (1990) for all our semi-analytical solutions (see Table 2), and thus independent of the considered opacity, radmc-3d uses the given dust opacities to calculate the attenuation factor self-consistently. This means that the same absorption efficiency (or opacity) is responsible for the attenuation of the ISRF. Despite its self-consistency, the latter may not be a true representation of the conditions in space, e.g, when pure materials like silicate are considered, and especially at UV wavelengths when small grains or large carbon-chain molecules, such as PAHs, are neglected.

Due to this difference, we expect moderate deviations between the methods. Furthermore, our numerical model uses a one dimensional spherical coordinate system in contrast to the plane-parallel approximation applied in the semi-analytical approach. The visual extinctions measured in one system must be converted to the other before a comparison is made (see e.g., Flannery et al. 1980; Röllig et al. 2007).

6.1. The radmc-3d code

We utilize version 0.39 of radmc-3d and take advantage of the multi-threading mode of the code. We consider the same ISRF and dust opacity tables as discussed in the previous sections. The dust temperature is then calculated as a function of radius in 1D for a spherically symmetric idealized molecular cloud. The radial density distribution of the cloud, ρ(r), follows a power-law profile and is given according to ρ(r) = ρ0(r/Rref)-2, where r, the radial position, runs from 0.1 AU (cloud center) to 6 pc (cloud edge). To ensure that the model probes both low and high visual extinctions, the radial coordinate grid is set logarithmically with 2000 resolution elements and finer resolution at the cloud edge. The reference radius Rref is taken as 0.5 pc and ρ0, the density at the reference radius, is 4 × 10-21 g cm-3. This is equivalent to a gas number density of nH = 1 × 105 cm-3 with a gas-to-dust ratio of 0.01. The solutions for the dust temperature are largely independent of the parameter choices (Rref, ρ0) or the profile, but gives us the necessary resolution and the dynamic range in AV that is desired in this work.

The model cloud core is isotropically irradiated by an external radiation field as shown in the top panel of Fig. 2 (green line). No internal heating source is considered. Besides the tracking of absorption and re-emission events of photon packages, which is inherently different from the semi-analytical calculations, radmc-3d is capable of taking into account (an)isotropic scattering of photons. We turn this mode off for better consistency with the semi-analytic model, but discuss and quantify the effect of scattering in Sect. 6.3. To reduce the intrinsic statistical noise of the Monte Carlo radiative transfer method, we set the number of propagated photon packages to 107, which is a relatively high number for an effectively 1D model. The estimated error will be on the order of 1/, where the photons will be spread out among the resolution elements causing the highest error to be at the core.

The radmc-3d code takes the opacity κν instead of the absorption efficiency Qν as input parameter. We convert between the quantities according to Eq. (3) and fix the grain radius to 0.1 μm as we do for the semi-analytical calculations. The considered opacity model types are listed in Sect. 4.2.

The code gives the radial distribution of dust temperatures for the opacity models as result. The radial position is converted to visual extinction by calculating the visual optical depth at each location from the cloud edge. The visual extinction at each radial position of the model grid is defined according to (6)where is the optical depth at λ = 5500 Å and is given by Eq. (4). In our model, is independent of the position in the cloud and can, therefore, be brought out of the integral. The integral then simplifies to a summation of the dust column density, whereas the optical depth is given by the product of the column density and the opacity. In order to have a one-to-one comparison with the semi-analytical models, we need to rescale the AV to account for the geometry, i.e., spherical geometry to semi-infinite slab. The rescaling factor is given by Röllig et al. (2007): (7)where μ = cosθ is the cosine of the radiation direction and τUV is the optical depth at UV, evaluated at 0.3 μm in the present work.

6.2. Numerically obtained dust temperatures

We show in Fig. 4 the dust temperature solutions obtained with the radmc-3d code for the dust materials graphite, silicate, and carbonaceous-silicate mixtures.

thumbnail Fig. 4

Dust temperature solutions obtained with the radmc-3d code for various grain materials as given in Table 1.

Open with DEXTER

thumbnail Fig. 5

Comparisons between the solutions obtained with the numerical (dashed lines) and the semi-analytical models (solid lines) for three different types of grain. The differences in temperature originates mainly from the underlying attenuation coefficient.

Open with DEXTER

We find similar trends and draw the same conclusions for the different dust materials as we did with our semi-analytical solutions. This means that the silicate dust types have, in general, a lower dust temperature at AV< 2 mag and a higher one at AV ≥ 150 mag as compared to the carbonaceous materials, which results in a smaller difference in temperature between the two boundaries. The temperature variation between all models at the edge (AV = 0.04 mag) ranges from 14.8 to 18 K, similar to the semi-analytical solutions, and when drawn further (to AV ~ 0) they completely match the semi-analytical solutions. The agreement also holds at high depths, that is, AV ≥ 150 mag.

The main differences as compared to the semi-analytical models come from the silicate grains. The temperature profiles of the silicate grains are shaped differently from the semi-analytical solutions. The silicate dust temperatures from radmc-3d can differ by up to 2.4 K, though mostly is less than 1 K. We attribute the differences obtained with radmc-3d to the attenuation of radiation, which is self-consistent instead of using the observed extinction curve. Therefore, for pure silicates, this misses the important carbon opacity contribution at short waves. We compare the differences in extinction curves in Appendix D, which clarifies the found discrepancies. As expected, the results from the two different methods do not match for pure silicate dust opacities, whereas they do match well for the mixed dust materials, especially Mixed A (OH94). The radmc-3d dust temperature solution for the Mixed A dust is practically identical to the semi-analytical solution. For three out of the seven models, i.e., Graphite, Silicate B (worst match), and Mixed A (best match), we show and highlight the differences in Fig. 5.

6.3. Scattering

The numerical approach allows us to change the geometry and include the scattering of interstellar photons on dust grains. We test and discuss the results at higher spatial dimension geometries in Appendix E. While the choice of geometry does not greatly impact the results, the expectation is that scattering may play a more appreciable role, especially for the attenuation of energetic photons (at short wavelengths). Here we discuss the impact of scattering.

Because we do not have the scattering opacities for all the dust models, the model Silicate B is alone computed and displayed as an example. From Fig. 6 we can see that the impact of scattering is such that the dust temperature rises at the cloud edge, while the temperature decreases at higher AV.

thumbnail Fig. 6

Dust temperatures with scattering. The dashed lines show the radmc-3d solutions (2D in blue, 3D in yellow) when scattering is treated. The midplane (θ = 0) values are adopted and for 3D, the median of φ is taken.

Open with DEXTER

The differences are larger for 3D than for 2D. This is easily explained by the fact that due to the higher dimensions, the path length of a photon increases (the photon can travel in more dimensions) and, therefore, gets more extincted and absorbed at low AV. This causes the lower AV temperature to be higher and the higher AV temperature to be lower. By not considering scattering, the underestimation from the 1D model is about 0.3 K (~2%) at low AV, whereas the overestimation, peaking at AV = 5 mag, is 1.5 K (~14%). At maximum AV (150 mag), the difference between 1D and 3D is about 0.8 K (~13%). We point out, once again, that in the semi-analytical calculations, scattering is considered through the extinction curve (see Sect. 4.3).

7. New parametric expression for Td

From the theoretical models that best approximate the observational results (these are the models that include the mixed dust material compositions, i.e., OH94, WD01, KYJ15, Or11) we create a simple and useful expression for the dust temperature. This is achieved by finding a function to match the model results. We find that the correlation is best reproduced by fitting a hyperbolic curve through the semi-analytical solutions. Our expression is formulated as (8)where the expression is scalable with χuv, the Draine UV field strength (Draine 1978). This matches our FISRF, the interstellar radiation field flux, where the mean flux in the radiation field with energy between 6 − 13.6 eV is G0 times 1.6 × 10-3 erg cm-2 s-1. For the Draine field, G0 = 1.7. In Fig. 7, we compare our expression with the expressions from past studies (i.e., Hollenbach et al. 1991; Zucconi et al. 2001; Garrod & Pauly 2011), which we describe in detail in Appendix F. In greyscales we highlight the variation that we get from our expression if we consider only thick (instead of thin) ices at AV ≥ 20 mag, or when we exclude the ice model, i.e., Mixed D (the turquoise line as given in the left panel of Fig. 3).

thumbnail Fig. 7

The parametric expression as a function of AV. We present our best fit to the theoretical calculations in this work and compare it to other parametric expressions found in the literature. The grey band illustrates the variation we get from considering thin or thick ices or from using three models to fit our expression instead of four.

Open with DEXTER

It is interesting to note that we find a χuv (or FISRF) with a power of 1/5.9 to best match the solutions, in close agreement with the analytical prediction of 1/6 by Draine (2011, see his Eq. (24.19)). The given formula is tested and validated for the range AV = 0.01 − 400 mag and χuv = 0.1 − 105 erg cm-2 s-1. The latter is shown in Appendix G. To account for ices, we fitted our line through bare dust results in the range 0 ≤ AV ≤ 6 mag and icy dust for AV> 6 mag (taking thick ices from OH94). This expression is constructed for an RV of 5, but we find only small differences when compared against RV = 3.1. The fitting function for RV = 3.1 is given in Appendix H, where the expression is also provided for other parameter choices.

8. Discussion: observed versus theoretical Td

We compare the observationally derived dust temperatures against the theoretical solutions and show them in Fig. 8.

thumbnail Fig. 8

Dust temperature solutions compared against observations. Top panel displays the obtained dust temperatures from semi-analytical solutions (Sect. 5) for various grain materials, without ices. The bottom left panel shows the observed dust temperatures against three parametric expressions found in the literature (Appendix F). The bottom right panel compares the observed dust temperatures against the radmc-3d solutions (Sect. 6).

Open with DEXTER

Here, we discuss our main findings with respect to observations.

8.1. Observations versus semi-analytic calculations

Our solutions provide a good match to the observational data as can be seen from the top panel of Fig. 8. Both the curvature as well as the values are captured well. Nevertheless, there appears to be a spread in the observational dust temperatures, particularly around AV ~ 10 mag. This may hint toward a spread in the underlying dust material (affecting opacity) or grain size distribution among sources, however, the uncertainties in the physical conditions and, especially, the ISRF may also just be the cause. Furthermore, when ices start to cover the surface of the dust, the composition of the refractory components tend to become less relevant. Indeed, when only selecting the LOS corrected observations and excluding the higher G0 environments, i.e., L1689B and CB17 (Roy et al. 2014; Schmalzl et al. 2014), the match improves greatly, as shown in Fig. 9. This indicates that the uncertainties arising from Qν at AV> 10 mag, i.e., when not considering ices, does not influence the dust temperature significantly.

We find that all models match the data reasonably, but the best match, with the lowest χ2, is Mixed C (KYJ15), the amorphous carbon-silicate mix with a carbon volume fraction of 0.48. This model has the highest carbon fraction among the considered mixed dust types. We also find that if we take ices into account in our models above an AV of 6 mag, the correlation to the observations improves.

8.2. Observations versus literature expressions

The parametric expressions from the selected literature studies do not agree well with the temperatures acquired observationally, with the exception of ZWG01 at AV ≳ 10 mag (Fig. 8 bottom left). Where the considered literature expressions were derived for a fixed parameter space, within the bounds of their respective studies, the expression attained from this work has a broad range of validity in AV and G0. We find that the greatest difference among the literature expressions arises from the adopted absorption efficiencies, i.e., Qν, which in all cases are power-law functions of frequency in the mentioned studies. Our calculations are an improvement in this domain, since we are not committed to a strict power-law, but rather adopt the measured and the intricately calculated opacities.

8.3. Observations versus RADMC-3D solutions

The temperatures obtained with the radiative transfer code are slightly lower, yet similar to the semi-analytical ones. Compared to observations, the dust temperatures obtained with radmc-3d seem to agree well especially at low AV, but are generally below the observed Td’s at AV> 2 mag, that is, except for Silicate B (Dr03). This indicates that using the observed extinction curve for the attenuation of radiation may indeed provide a more accurate picture at AV> 2 mag than the self-consistent extinction from the adopted opacities as radmc-3d does. Moreover, along the LOS a mixture of dust types and sizes will eventually be responsible for the attenuation. Ideally, one should modify the opacity input for radiative transfer calculations to take the components of the extinction curve (PAHs, carbon grains and silicates) into account. It is, on the other hand, also true that the dust temperatures from observations have the difficulties of the LOS. Even with the corrections, these may cause an overestimation of the dust temperature, especially in embedded regions (cf. Schnee et al. 2006, 2008; Roy et al. 2014). Thus, a lower dust temperature as compared to the observations at higher AV is expected.

8.4. Active star-forming regions

Our expression for the dust temperature, i.e., , scaled by the ISRF (see Fig. G.1), is in agreement with the observed dust temperature of ρ Oph A. For an environment with a 103G0 and at AV = 20 mag, our expression gives a dust temperature of 22.3 K, similar to the results obtained from line emission of H2O2 (22 ± 3 K, Bergman et al. 2011b) and not far from HO2 measurements (16 ± 3 K, Parise et al. 2012). While for an AV = 100 mag, we obtain 17.6 K from our expression, similar to the results retrieved from D2CO measurements of ρ Oph A (17.4 K, Bergman et al. 2011a).

thumbnail Fig. 9

Semi-analytical solutions compared against the LOS corrected observed dust temperatures. Only sources consistent with G0 ~ 1 are considered. The match improves markedly.

Open with DEXTER

9. Conclusions

The motivation of this study was to find and provide a parametric expression for the dust temperature to use in numerical simulations and chemical models. To this end, we calculated the dust temperature from basic principles for dust in thermal equilibrium by considering in detail the ISRF, the attenuation of radiation, and the dust opacities. We did this for various grain material compositions, i.e., for graphite, silicates SiO2 and MgFeSiO4, and carbonaceous silicate mixtures. We compared our calculations against solutions obtained with the Monte Carlo radiative transfer code radmc-3d and against recent observational results from Herschel that we collected from the literature.

We find that our semi-analytical solutions as well as our numerical solutions match the range of observed dust temperatures well at low and at high optical depths and also captures the overall extinction dependence (barring uncertainties in the ISRF). Mixed carbonaceous silicate dust material compositions match the observed temperatures of starless regions better than pure materials and give a narrow range temperature solution between 15.2 and 17.7 K at the edge, that is, AV = 0.04 mag. This conforms to 5 − 6 K around an AV = 150 mag. However, the dust surface should be covered by ices at high AV making the composition of the refractory components less relevant in this regime.

Considering the impact of ices, we find that ice formation changes the opacity of dust significantly enough to reduce the net cooling at AV> 10 mag. This allows the dust to be slightly warmer (~15% for thick ice and around 8% for thin ice) in highly embedded regions, which may be crucial in avoiding the freeze-out of H2 molecules in models (typically below Td ≃ 7 K). Ice formation helps to raise the dust temperature at high optical depths (AV ≳ 20 mag) by about 0.5 to 1 K, depending on ice thickness. The ices also aid in flattening the temperature profile, which helps in explaining the near (semi-log) linear profile inferred from the observed dust temperatures. We find that with ices (at AV ≳ 6 mag) our models give a better match to the observed dust temperatures.

From our best matching lines, we provide an analytical expression as a function of AV, given by Eq. (8) (detailed further in Appendix H), which can be scaled by the ISRF.


2

When simply using SED fitting, the average dust temperature along the sight line is obtained.

3

Here we mean lower than unity emissivity, but with ν dependence.

4

The ZWG01 power-law best matches the data for the range 5 − 70μm.

5

A monolayer is when the whole dust surface is covered by one molecule thick (~3 Å) layer of ice.

Acknowledgments

We thank M. Köhler, N. Ysard, and C. Ormel for providing their opacity data tables. We thank A. Ivlev for his contribution on cosmic rays and M. Röllig for sharing his geometry conversion code. S.H. especially thanks J. Steinacker and W.F. Thi for their great interest and constructive comments about this work. P.C., M.S., and S.C. acknowledge the financial support of the European Research Council (ERC; project PALs 320620). S.C. is also supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.017).

References

Appendix A: Scaling with cosmic rays

As a test to the impact of higher CR rates, for example in CR-dominated regions (e.g., Papadopoulos et al. 2011; Bayet et al. 2011), the calculations are also performed with elevated CR rates. In Fig. A.1 the dust temperature solution for a 1000 times higher CR rate, i.e., ζH2 = 5 × 10-14 s-1 is presented. As noted before, only the heating by UV photons created from CRs are considered, not the direct CR impact or other CR processes.

With an increased CR rate, the resulting dust temperatures above an AV ≈ 5 mag are much higher, while at lower extinctions this is not the case. CRs essentially set a floor temperature for the dust grains, below which they cannot cool, though, this could be overcome by the magnetic mirroring process (Padovani & Galli 2013). It is interesting to note that the higher grain temperatures, that is seen from the observations around an AV ~ 10 mag, can be reproduced by assuming a 1000 times higher CR rate.

thumbnail Fig. A.1

The dust temperatures from semi-analytical solutions for a 1000 times higher cosmic ray rate.

Open with DEXTER

Appendix B: The ISRF

The full expression for the ISRF that is used in the calculations of this work is reported here. The five modified black bodies from the work by Zucconi et al. (2001) is adopted, except for the MIR range which is changed as was presented in Eq. (2) into a sixth modified black body. To this function the UV part of the spectrum is added. The part of the spectrum without UV is given by the six modified black bodies as (B.1)where the values for Wi and Ti are given in Table B.1. The included values for the MIR match the 10 micron emission coming from hot dust, but is highly smoothed out. The UV part of the spectrum is adopted from Draine (1978), rewritten in the current form to match the units, that is given as (B.2)The combined radiation field is given by in units of erg s-1 cm-2 sr-1 Hz-1.

Table B.1

Parameters for the ISRF.

thumbnail Fig. B.1

Same Figs. 2 and 3, ISRFs (top) and dust temperature solutions (bottom). Here, the Galprop ISRF (yellow dotted line, top) is used for calculating the dust temperatures.

Open with DEXTER

Appendix B.1: Testing with a different ISRF

In order to see how strong a dependence on the chosen ISRF there is, an online available, observationally constrained, ISRF is taken from Galprop (Strong et al. 2007; Ackermann et al. 2012) and used for the computations. The Galprop code calculates and extrapolates the ISRF for every part of the Milky Way. For this exercise, the Galactic plane value at a distance of 8.5 kpc is adopted. In Fig. B.1, we see the results for this ISRF.

The dust temperatures are slightly higher at the edge, but lower at high AV. The first is due to a higher radiation flux at optical and UV wavelengths from the Galprop ISRF, while the latter follows from the excluded flux at FIR and microwave wavelengths (notice the missing CMB part in the Galprop ISRF).

We point out that a recent study by Steinacker et al. (2016b) derives a factor 4 higher flux for the local MIR and FIR components of the ISRF of L1689B.

Appendix C: The temperature impact of ices

thumbnail Fig. C.1

Td as a function of AV with ices. Similar to the right panel of Fig. 3, the impact of ices on the dust temperature is displayed. Top panel shows the OH94 models, but coagulated for 105 yr, while the bottom panels show the KYJ15 (left) and the Or11 (right) models.

Open with DEXTER

Similar to what is shown in Fig. 3, right panel, the impact of ices on the dust temperature from two other studies, i.e., from Or11 and KYJ15, is displayed here. See Fig. C.1. The sole exception is that the Or11 model does not reproduce the suggested trend at AV< 0.8 mag, where the ice covered dust results in a higher Td. However, this is not relevant in a realistic case, because ices should not be present in this regime. While in his work Or11 does provide two other types of ice mixtures which yield more consistent results at low AV. It is advisable to note that, as stated before, the KYJ15 and the Or11 base ice models are not entirely separable from coagulation. The KYJ15 models include four big grain aggregates, while the Or11 models perform a minimum of 30 Kyr coagulation, which is a relatively short timescale. Both of these models demonstrate the same trend as indicated by the OH94 models in Fig. 3.

Appendix D: Comparing extinction curves

The attenuation coefficient γλ = Aλ/AV of the extinction curve from Mathis (1990) and the self-consistent opacities of radmc-3d are compared, see Fig. E.

thumbnail Fig. D.1

Extinction coefficients Aλ/AV as a function of wavenumber 1 /λ (μm-1). A comparison is made between Mathis (1990), black points, and self-consistent opacities, colored lines, used by radmc-3d.

Open with DEXTER

The Mixed A (OH94) model shows a very similar wavelength dependence in the range as the observed curve. For this model, the resulting temperature profiles from the semi-analytical and the radmc-3d methods are practically identical (Sect. 6.2). The extinction curves resulting from pure silicate materials are very different than the observed extinction curve. At long wavelengths the differences are more complicated and highly wavelength dependent.

Appendix E: Geometrical dependence

The higher dimension models are set up similarly to the 1D model of radmc-3d (Sect. 6), a spherically symmetric cloud core with a power-law radial density profile (see Sect. 6.1), but in spherical coordinates (r,θ,φ). In Fig. E.1 the dimensional dependence is shown for two separate models: Mixed A (OH94) and Silicate B (Dr03).

thumbnail Fig. E.1

Dust temperature in different geometries. The top panel shows the dimension impact for the Mixed A dust material type. The bottom panel shows this for the Silicate B dust material. Blue line shows the 1D model results, green line the 2D results, and black line the 3D results. The greyscales illustrate the variation from all the φ angles in the 3D model.

Open with DEXTER

With higher dimensions, due to the increased number of grid cells, and without changing the number of photon packages, one loses precision. This results in a higher noise with increasing AV. One can particularly notice this from the 3D curves. Since there are the angles, θ and φ, to be considered, for simplicity, the midplane value for θ at higher dimensional geometries are adopted. However, all φ angles for the 3D model is given in greyscales. The median radial temperature profile along the φ coordinate is drawn in black. Between all dimensions, and for both dust models, the difference is always less than 0.7 K, that is, when considering the median value for the 3D model. There seems to be some geometrical dependence, but this dependence is small. The 1D model is consistently higher in temperature than the higher spatial dimension solutions. Due to the introduced noise, and the finer details of the code, it is hard to state if there is any systematic variation between the 2D and the 3D results.

Appendix F: Analytical expressions from the literature

There are several other studies that describe a method to calculate the dust temperature from first principles, assuming thermal equilibrium, and provide simple parametric expressions that can be used in astrophysical models. We report here solutions by three different groups in order to have a basis for comparison. The parametric expressions that will be discussed are acquired from Hollenbach et al. (1991), Zucconi et al. (2001), and Garrod & Pauly (2011), which we had identified as HTT91, ZWG01, and GP11, respectively.

Appendix F.1: The HTT91 expression

The solution by Hollenbach et al. (1991) assumes a one-sided slab geometry and the expression for the dust temperature is formulated as (F.1)Here, ν0 represents the main absorbing frequency over the visual and UV wavelengths and is suggested as ν0 = 3 × 1015 s-1, G0 is the UV flux in terms of the Habing field (Habing 1968), τ100 is the emission optical depth at 100 μm, and T0 is the equilibrium dust temperature at the cloud edge due to the unattenuated incident FUV field alone. Given the condition for T0, this parameter equates to  K. HTT91 assume that the incident FUV flux equals the outgoing flux of dust radiation at T0, such that . Knowing T0 fixes τ100 to a value of 0.001.

thumbnail Fig. F.1

Comparison between three analytical expressions for Td as found in the literature.

Open with DEXTER

The HTT91 expression is a useful function which combines the heating by UV, CMB, and the re-processed IR. It is a function of G0 and AV, but the AV dependence is only for the attenuation of the incident UV field. The authors assume a fixed scaling of 1.8 for the UV attenuation with frequency, which is roughly a Planck averaged value over the extinction curve. The expression is scalable with the background UV field such that it can be applied to photodissociation regions. The assumed one-sided slab geometry will not affect the solution for the cloud edge, but will shape the transition region and center temperature, unless large enough column densities are considered. Qν is set to unity for every ν>ν0, and to ν/ν0 when this is not the case. This gives a linear scaling with frequency for the absorption efficiency. In reality, however, the extinction features from the type of grain material should make it highly non-linear. The HTT91 expression is often adopted and highly referenced in other studies.

Appendix F.2: The ZWG01 expression

The expression provided by Zucconi et al. (2001) is defined as (F.2)where VNIR, MIR, and FIR stand for contributions from the visual plus near-infrared (NIR), mid-infrared, and far-infrared, respectively. The individual dust temperatures are given as (F.3)The ZWG01 solution for the dust temperature is designed for the range 10 ≲ AV ≲ 400 mag. An additional expression is provided for AV = 1 − 2 mag, which we will not consider hereafter because of the narrow range.

The ZWG01 expression is based on the observed dust temperature of L1544 at various AV. These authors solve the thermal balance using observed quantities and provide a parametric fit to the solution. The dust temperature expression is solely a function of AV. The thermal balance is solved without considering the UV field thereby purely concentrating on the visual and IR part of the spectrum. This is justified, because the focus is on regions with high AV where UV light will be extincted. The ISRF spectrum is represented by a sum of five modified black bodies employing the data of Black (1994) and Mathis et al. (1983), and a power-law for the MIR. The dust opacities are adopted from OH94 for the thin ice, 105 yr coagulated case. ZWG01 do this by considering the opacities of a few peak frequencies and assume a power-law behavior rather than using the entire opacity table. They scale Qν as (ν/νpeak)η, where η is obtained from the OH94 models.

The authors note that can be scaled with the background radiation field by multiplying the dust temperature with a factor of the field intensity to the power 1/5.6 (Zucconi et al. 2001; Galli et al. 2002). This notion presumes, to an acceptable degree, that the whole spectral shape of the ISRF scales in the same manner, albeit the CMB should certainly not.

Appendix F.3: The GP11 expression

The expression by Garrod & Pauly (2011) is presented as (F.4)which is only valid in the range 0 ≤ AV ≤ 10 mag.

The GP11 formulation complements the low AV regions and is designed to be used in combination with . The authors adopt a scaling of Qνaν2 for the efficiency following Krügel (2008), and fix a, the grain radius, to 0.1 μm. This gives a steep scaling to Qν, which is only expected at long wavelengths. Their choice results in a high value of Qν below λ ≲ 30 μm and a low value in the submillimeter. This results in a higher dust temperature at low AV. A critical note is that GP11 adopt the right-hand side of Eq. (1) from ZWG01, which means that UV is not taken into account into the ISRF either. This is a rather important omission at low AV. GP11 also add a corrective value of 0.316 K to their solution to ensure continuity with ZWG01 at AV = 10 mag. For Dν(AV), see Eq. (1), the authors follow Cuppen et al. (2006) thereby using the tabulated values of Mathis (1990), times a factor 0.8, to evaluate Aν/AV. This is an improvement over fixing it to a single value as was the case in HTT91 and also over the three-part power-law method as was favored by ZWG01.

Appendix F.4: Comparison between expressions

The three dust temperatures derived from the three expressions in the previous sections are presented in Fig. F.1.

In this figure, we compare the three solutions for the ISRF as provided by the respective studies, which differ from each other.

It is clear that the differences between the models are quite large. The temperatures at cloud edge, which is the simplest condition because of no depth or geometry dependence, already vary between 12.2 K and 18.9 K. This disparity mainly arises from a different treatment of absorption efficiency (Qν, see Fig. 2). At high optical depths, that is, AV = 150 mag, the solutions of HTT91 and ZWG01 result in the dust temperatures of 4 and 6 K. Here, the difference mostly comes from the different treatment of the FIR part of the spectrum. Among the three expressions, the dust temperature profile from the ZWG01 solution is more in line with the observations.

In Fig. F.1, we also show one extra curve, which is given by the red dot-dashed line. This is taken from Van Borm (2016). Using a similar approach to ZWG01, the author derived analytically the dust temperature with a better treatment for the ISRF. It is worthwhile to note that this parametric expression reproduces the observed dust temperatures better than the ones from earlier studies. The full expression for can be found in Chapter 3.2.H of that work.

Appendix G: Scaling with background field

The thermal balance equation (Eq. (1)) is also solved for various different radiation field strengths. The assumption is that FISRF (not including CMB) scales proportionally to the UV field, which is characterized by χuv. Here, the semi-analytical solutions are compared against the parametric fitting function for various radiation field strengths, i.e., χuv = 0.1, 1, 10, 102, 103, 104, 105 (Eqs. (8) and (H.1)). The CMB part of the radiation field spectrum is not scaled with this factor for the semi-analytical calculations. The example given in Fig. G.1 is for mixed grains, bare dust. The parametric fits match the semi-analytical solutions excellently. We note that even at AV = 150 mag, heating by the ISRF dominates over the heating by the CMB.

thumbnail Fig. G.1

Dust temperatures from semi-analytical solutions (dashed lines) and from parametric fits from Eq. (8) (solid lines) for various ISRF strengths. The correspondence is noteworthy.

Open with DEXTER

Appendix H: Additional expressions

Table H.1

Three- and five-parameter fit coefficients

Additional parametric expressions constructed from the fits to the semi-analytical solutions are given here. The fits are provided for every type of grain material (graphite, silicate, mixed), with and without ices, and for RV = 3.1 (fiducial model is RV = 5).

A slightly more accurate and an extended version of the expression in Eq. (8) is provided here, i.e., (H.1)where X = log 10(AV). The coefficients α, β, γ, (δ, ϵ) are given in Table H.1 for the three-parameter and the five-parameter fits.

In Fig. H.1, the difference between the various parametric fits, considering the 5-parameter fits, are highlighted and plotted on top of the observed dust temperatures. The observed Td’s, consisting of only LOS corrected data here, are fit with the hyperbolic function of Eq. (H.1), rather than the linear fit given in Sect. 2.

thumbnail Fig. H.1

Comparison between selected (graphite, mixed, and with ices) hyperbolic fits overlaid on top of observations (red crosses). Observational data only include LOS corrected points. Including ices into the mix brings the solutions nearer to the observed dust temperatures.

Open with DEXTER

Appendix I: Scaling with redshift

In order to scale our expression with redshift z, the increase in CMB temperature has to be considered. The CMB temperature will rise according to the relation . This will change the dust temperature solutions in the following way (I.1)see, for example, da Cunha et al. (2013).

This simple prescription is not taking into account the changes in the UV field at high redshift due to different star formation rates and various other processes.

All Tables

Table 1

Considered dust material types adopted from literature.

Table 2

Attenuation coefficient γλ (RV = 5).

Table B.1

Parameters for the ISRF.

Table H.1

Three- and five-parameter fit coefficients

All Figures

thumbnail Fig. 1

Observed dust temperatures from eight independent studies. The dust temperature is plotted as a function of visual extinction. A least squares semi-log linear line is fit through the data as given by the dashed black line.

Open with DEXTER
In the text
thumbnail Fig. 2

Parameters for the dust equilibrium. Top panel, the ISRF intensity as a function of frequency. The green line represents the adopted ISRF in this work. Bottom left panel, experimental and calculated absorption efficiencies for various grain materials. Scattering is not included in these. Bottom right panel, extinction curves from various studies. The filled black circles show the observed data from Mathis (1990). The black solid line is a fit to the data as given by Cardelli et al. (1989), which is adopted in this work. The sub panel zooms in at the lower wavenumbers where the adopted extinction curve below 0.15 μm-1, given by the black solid line, is interpolated from Mathis (1990).

Open with DEXTER
In the text
thumbnail Fig. 3

Dust temperature solutions. The left panel displays the obtained dust temperatures for various grain materials, which have no ices. The right panel shows the impact of ice formation.

Open with DEXTER
In the text
thumbnail Fig. 4

Dust temperature solutions obtained with the radmc-3d code for various grain materials as given in Table 1.

Open with DEXTER
In the text
thumbnail Fig. 5

Comparisons between the solutions obtained with the numerical (dashed lines) and the semi-analytical models (solid lines) for three different types of grain. The differences in temperature originates mainly from the underlying attenuation coefficient.

Open with DEXTER
In the text
thumbnail Fig. 6

Dust temperatures with scattering. The dashed lines show the radmc-3d solutions (2D in blue, 3D in yellow) when scattering is treated. The midplane (θ = 0) values are adopted and for 3D, the median of φ is taken.

Open with DEXTER
In the text
thumbnail Fig. 7

The parametric expression as a function of AV. We present our best fit to the theoretical calculations in this work and compare it to other parametric expressions found in the literature. The grey band illustrates the variation we get from considering thin or thick ices or from using three models to fit our expression instead of four.

Open with DEXTER
In the text
thumbnail Fig. 8

Dust temperature solutions compared against observations. Top panel displays the obtained dust temperatures from semi-analytical solutions (Sect. 5) for various grain materials, without ices. The bottom left panel shows the observed dust temperatures against three parametric expressions found in the literature (Appendix F). The bottom right panel compares the observed dust temperatures against the radmc-3d solutions (Sect. 6).

Open with DEXTER
In the text
thumbnail Fig. 9

Semi-analytical solutions compared against the LOS corrected observed dust temperatures. Only sources consistent with G0 ~ 1 are considered. The match improves markedly.

Open with DEXTER
In the text
thumbnail Fig. A.1

The dust temperatures from semi-analytical solutions for a 1000 times higher cosmic ray rate.

Open with DEXTER
In the text
thumbnail Fig. B.1

Same Figs. 2 and 3, ISRFs (top) and dust temperature solutions (bottom). Here, the Galprop ISRF (yellow dotted line, top) is used for calculating the dust temperatures.

Open with DEXTER
In the text
thumbnail Fig. C.1

Td as a function of AV with ices. Similar to the right panel of Fig. 3, the impact of ices on the dust temperature is displayed. Top panel shows the OH94 models, but coagulated for 105 yr, while the bottom panels show the KYJ15 (left) and the Or11 (right) models.

Open with DEXTER
In the text
thumbnail Fig. D.1

Extinction coefficients Aλ/AV as a function of wavenumber 1 /λ (μm-1). A comparison is made between Mathis (1990), black points, and self-consistent opacities, colored lines, used by radmc-3d.

Open with DEXTER
In the text
thumbnail Fig. E.1

Dust temperature in different geometries. The top panel shows the dimension impact for the Mixed A dust material type. The bottom panel shows this for the Silicate B dust material. Blue line shows the 1D model results, green line the 2D results, and black line the 3D results. The greyscales illustrate the variation from all the φ angles in the 3D model.

Open with DEXTER
In the text
thumbnail Fig. F.1

Comparison between three analytical expressions for Td as found in the literature.

Open with DEXTER
In the text
thumbnail Fig. G.1

Dust temperatures from semi-analytical solutions (dashed lines) and from parametric fits from Eq. (8) (solid lines) for various ISRF strengths. The correspondence is noteworthy.

Open with DEXTER
In the text
thumbnail Fig. H.1

Comparison between selected (graphite, mixed, and with ices) hyperbolic fits overlaid on top of observations (red crosses). Observational data only include LOS corrected points. Including ices into the mix brings the solutions nearer to the observed dust temperatures.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.