Full SED fitting with the KOSMA-\mathrm{\tau} PDR code - I. Dust modelling

Full SED fitting with the KOSMA- PDR code - I. Dust modelling

M. Röllig I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany    R.Szczerba N. Copernicus Astronomical Center, Rabianska 8, 87-100 Toruń, Poland    V. Ossenkopf I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany    C.Glück I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany roellig@ph1.uni-koeln.de
Key Words.:
offprints: M. Röllig,


Aims:We revised the treatment of interstellar dust in the KOSMA- PDR (photo-dissociation region) model code to achieve a consistent description of the dust-related physics in the code. The detailed knowledge of the dust properties is then used to compute the dust continuum emission together with the line emission of chemical species.

Methods:We coupled the KOSMA- PDR code with the MCDRT (multi component dust radiative transfer) code to to solve the frequency-dependent radiative transfer equations and the thermal balance equation in a dusty clump under the assumption of spherical symmetry, assuming thermal equilibrium in calculating the dust temperatures, neglecting non-equilibrium effects. We updated the calculation of the photoelectric heating and extended the parametrization range for the photoelectric heating toward high densities and UV fields. We revised the computation of the H formation on grain surfaces to include the Eley-Rideal effect, thus allowing for high-temperature H formation.

Results:We demonstrate how the different optical properties, temperatures, and heating and cooling capabilities of the grains influence the physical and chemical structure of a model cloud. The most influential modification is the treatment of H formation on grain surfaces that allows for chemisorption. This increases the total H formation significantly and the connected H formation heating provides a profound heating contribution in the outer layers of the model clumps. The contribution of polycyclic aromatic hydrocarbons (PAH) surfaces to the photoelectric heating and H formation provides a boost to the temperature of outer cloud layers, which is clearly traced by high- CO lines. Increasing the fraction of small grains in the dust size distribution results in hotter gas in the outer cloud layers caused by more efficient heating and cooler cloud centers, which is in turn caused by the more efficient FUV extinction.


1 Introduction

Photo-dissociation regions (PDRs) are interstellar neighbors of HII regions that are already shielded from extreme ultraviolet (EUV) photons with an energy sufficient to ionize atomic hydrogen, but where the physical and chemical conditions of the atomic and molecular interstellar medium (ISM) are still governed by the remaining far-ultraviolet (FUV) radiation of nearby massive stars (Hollenbach & Tielens, 1999). Here, the transition from the atomic to the molecular ISM takes place, giving rise to a rich astrochemical environment. The combination of a chemically rich pool of species that includes a large reservoir of free electrons and strong energetic excitation conditions produces a wealth of spectroscopic emissions, carrying information on the internal chemical and physical structure of the PDR.

Comparing the results of theoretical and numerical calculations to the observed spectral line emission of PDRs is a frequently used method to infer their local physical conditions. Numerical models to describe the physical and chemical processes in molecular clouds have been successfully used for many years. They involve simultaneously solving the problem of a) radiative transfer for a cloud of gas and dust; b) the chemical structure by balancing formation and destruction processes for all included species; c) determining the thermal structure of the cloud by balancing all important heating and cooling processes.

Properly taking into account all relevant chemical and physical processes is extremely time-consuming in terms of computational power and usually leads to a large number of simplifying assumptions in the treatment of the physics and the chemistry. The exact choice of these simplifications depends on the field of application, on the expertise of the modeller, and on the available physical and chemical knowledge. Many processes are still not fully understood if not largely unknown. A prominent example is the description of interstellar dust (see Draine, 2003a, for a review on interstellar dust). Many important physical and chemical processes require a good knowledge of the properties of interstellar dust grains (Abel et al., 2008). Chemical reactions on the surface of dust grains appear to be the only efficient formation route for a number of important astrochemical species (Garrod et al., 2008; Hall & Millar, 2010). The specifics of their reaction kinematics depend to a large degree on the material properties and the structure of the dust grains. Porous, spongy grains provide a large surface for chemical reactions but might hinder the release of the newly formed species into the gas phase (Leger et al., 1985; Roberts et al., 2007; Taquet et al., 2012). Not only the material and shape, but also the size of the grains might be important. Very small grains and very large molecules, such as polycyclic aromatic hydrocarbons (PAHs), efficiently contribute to the gas heating by means of the photoelectric (PE) effect (Bakes & Tielens, 1994; Weingartner & Draine, 2001a) while larger grains are dominating the scattering and absorption of FUV photons. For a review on PAHs see Tielens (2008).

The rapid development in observational techniques over the last two decades left us with a vast amount of spectroscopic data that defies reproduction with simple numerical PDR models. We are far away from being able to really understand the conditions in any observed PDR (Röllig et al., 2007). In this paper we describe a number of steps that have been taken to increase the modeling power of the KOSMA- PDR model code. We tried to evolve the code toward a consistent treatment of dust physics and chemistry. By calculating the optical and energetic properties of a model cloud with a given dust composition we are now able to describe the full spectral emission characteristic of a model cloud including dust continuum and line emission111Other PDR model codes are also able to perform wavelength-dependent UV radiative transfer. Compare for example Le Petit et al. (2006); Abel et al. (2005)..

2 Code description

2.1 The KOSMA- PDR model code

The KOSMA- code (Röllig et al., 2006), a welltested and mature PDR model code (Röllig et al., 2007), applies spherical geometry to the problem of simultaneously solving the chemical and energy balance in an interstellar molecular cloud. KOSMA- is equipped with a modular chemical network, i.e., chemical species can easily be added or removed from the network and the network will rebuild dynamically, including isotopic chemistry of \ce^13C and \ce^18O. After obtaining the physical and chemical structure of the model cloud, a radiative transfer code is applied to calculate the resulting emissivities for the spectral line emission. Because of the spherical geometry of the model clump, radiation can reach a point inside the clump from any direction.

2.2 Multi-component dust radiative transfer model

The multi-component dust radiative transfer (MCDRT) code allows one to solve the frequency-dependent radiative transfer equations and the thermal balance equation in a dusty clump under the assumption of spherical symmetry (Yorke, 1980). A detailed description of the code’s main features is given in Szczerba et al. (1997). The code includes isotropic scattering, which is important at high optical depths. Further modifications to that version of the code have been made. First, we have added the possibility to solve the radiative transfer simultaneously for a number of dust sorts (NDS) that may or may not be co-spatial, and obey any dust size distribution. The grain size distribution is given by the general relation valid for grains with radius ;


where: is the number density of grains with size and is the number density of H nuclei (in both atoms and molecules: ). The function has a modular structure and can be changed easily to a different form. Second, because we adapted the code to the conditions typical for PDRs it was necessary to introduce an option for switching off the central source. In addition, because intense sources of an interstellar radiation field are ubiquitous in PDRs, we changed the spatial distribution of impact parameters (originally logarithmically spaced beginning from the inner shell of the clump) along which the ray equations are integrated. The end result of this code is the infrared radiation flux emitted by the clump, which can be compared to the radiation observed from real PDRs. However, as a byproduct, the code provides the mean intensity of the radiation field, and the dust temperature, , which for now is computed from the assumed thermal equilibrium for each dust species and dust grain size at each radius of the clump. These two quantities (mean radiation field and dust temperature) are now used in the KOSMA- PDR model code to obtain a self-consistent gas-dust model.

2.3 Dust models

For the purpose of this paper we selected three dust models of interstellar dust from Weingartner & Draine (2001b) (hereafter WD01), and the MRN dust model (Mathis et al., 1977) for comparison. The MRN interstellar dust model was constructed by fitting a dust model to the Galactic mean extinction curve. It consists of two separate dust populations, one for "astronomical silicates" (sil) and one for graphite (gra). In the MRN model the grain size distribution is identical for both dust components and is given by


where i stands for sil and gra, , and (WD01).

Recently gathered observational constraints led WD01 to construct new models for interstellar dust. The authors proposed that interstellar dust is composed of four main components: astronomical silicates, graphite, and populations of very small grains that consist of neutral (PAH) and ionized (PAH) particles. The grain size distributions for these dust components follow the general definition given by Eq. (1) and their functional forms ( in Eq. (1)) are given by Eqs. (2–6) of WD01. The authors employ a power-law form for the large grains and a log-normal distribution for the very small grain population. Free parameters of the functions were determined by the best fit to the average extinction with ) of 3.1, 4 and 5.5 (see WD01 for details). is the absolute visual extinction at V = 5500 Å and . The three dust models from WD01 correspond to lines 7, 21, and 25 from Table 1 in Weingartner & Draine (2001b) and are denoted WD01-7 (), WD01-21 (), and WD01-25 ().

To compute the extinction efficiency and albedo for each dust component, we assumed that grains are spherical and used the Mie theory (Bohren & Huffman, 1983). For silicates and graphite we used the dielectric constants from Draine (2003b), while for very small grains we followed the approach given by Li & Draine (2001). The minimal and maximal grain size for the size distribution of each dust component are given in Tab. 4. For very small grains we used 17 grain sizes, while for bigger grains we divided each dex of grain sizes into ten sizes, keeping equal distance in . The wavelength coverage used in the MCDRT code extends from 10 Å to 3000 m and is split into 333 wavelengths.

3 Impact on the PDR model

In the following we describe how the new, full radiative transfer (RT) computation compares to the old approximation used in the KOSMA- model. For this purpose we evaluated a reference model clump with varying FUV radiative transfer and different dust models. We kept the following model parameters constant for all models: the total clump mass 10 M, the total surface gas density  cm, and the total FUV field in units of the Draine field (Draine, 1978). We applied a power law density gradient with power index 1.5 with a constant central gas density for radii smaller than . This implies . We assumed a total cosmic ray ionization rate of molecular hydrogen s. The chemistry is based on the UMIST 2006 database for astrochemistry UDfA222www.udfa.net (Woodall et al., 2007), using 47 chemical species and 490 reactions in total.

3.1 UV continuum radiative transfer

3.1.1 Dust extinction

For any given line of sight, the optical depth at FUV wavelengths (Sternberg & Dalgarno, 1989), with being the effective dust grain extinction cross section in the FUV, and the total hydrogen column density . In the previous version of the KOSMA- code, the local radiation field was calculated by computing the optical depth as a function of angle and then performing an angular averaging of the resulting intensity over the full solid angle. A detailed description of the FUV transfer in the spherical model has been given by Störzer et al. (1996).

For a constant dust-to-H ratio the total extinction cross section per H nucleus can be derived from


where is the extinction cross section of a single dust particle with size at wavelength . Often the absolute extinction in magnitudes is used


By introducing the average mass extinction coefficient for each dust component, i.e., the mean extinction cross section per dust mass,


where is the bulk dust grain density (i.e., the density of material of which dust of given sort is composed) and is volume density of the -th dust sort, the optical depth at any given point of the model clump can be written as


Relating this to the gas density provides


where , and is mass of a single H atom. Using as the total dust-to-H mass ratio, and we can write




For a constant dust-to-H mass ratio the integration over in Eq. (8) can be performed giving an H column density using


i.e., .

3.1.2 Influence of the dust composition

Figure 1: Extinction cross section per H nucleus for the WD01-7 dust model (dashed lines) and for the MRN model (solid lines). The total cross sections are represented by thick lines. The contribution from each dust component is shown by red lines for silicates, blue lines for graphite, and green lines for very small grains.

Figure 1 shows the total extinction cross section per H nucleus computed for the WD01-7 (thick dashed line) and the MRN model (thick solid line) of interstellar dust. The other lines show contributions from each dust component incorporated into the WD01-7 model (silicates – dashed, red in the electronic version; graphite – dashed, blue in the electronic version; and PAH+PAH — dashed green line in the electronic version), and the solid lines present contribution from silicates (red) and graphite (blue) in the MRN model. One can see that in the WD01 model populations of very small grains (PAH and PAH) are responsible for a large part of the total extinction, especially for the 2200 Å band, while large grains dominate the total mean Galactic extinction in the MRN model.

Figure 2: Relative extinction for the dust models of MRN (solid) , WD01-7 (long-dashed; red in electronic version); WD01-21 (short-dashed; green in the electronic version) and WD01-25 (dot-dot-dashed; blue in the electronic version) of interstellar dust. The vertical dotted lines at wavelengths of 912 and 2066 Å, corresponding to 13.6 and 6 eV, show the range of averaging. The horizontal lines show for the considered models.

Fig. 2 shows for the dust models tested in this paper: MRN (solid), WD01-7 (long-dashed; red in the electronic version), WD01-21 (short-dashed; green in the electronic version), and WD01-25 (dot-dot-dashed; blue in the electronic version). The FUV range relevant for photoelectric heating and photo-dissociation reactions extends between 13.6 and 6 eV333H only absorbs FUV photons via Lyman and Werner electronic transitions in the 912-1100 Å range.. Consequently, we define an FUV-to-V color as  = where the averaging is performed over an energy from 6 to 13.6 eV and is the visual extinction. The vertical dotted lines in Fig. 2 indicate this range of averaging and the horizontal lines show for the considered models. Note that the MRN and WD01-7 dust models have the same value of . All values of the average relative dust extinction are collected in Tab. 4. From Fig. 2 it can be seen that is significantly decreased in the WD01-21 and WD01-25 dust models. This is consistent with the higher values of , i.e., a shallower gradient of toward shorter wavelengths causing a lower UV extinction for the same (see Fig. 2). This is the main source for the different UV attenuation properties of different dust types.

3.1.3 Influence of radiative transfer and scattering

Figure 3: Mean intensity for a model clump of 10 M,  cm, and . The lines show for radii in steps of 10% of .

As a result of the new radiative transfer computations we obtain the full, wavelength-dependent FUV radiation field at all clump radii, where is the specific intensity at wavelength . In Fig. 3 we plot for radii in steps of for the WD01-7 model. The effect of the prominent 2200 Å “bump” is plainly visible in the enhanced reduction of with decreasing radius.

The FUV optical depth


can be described by the effective dust extinction cross section


For adopted dust models, the obtained values of are collected in Tab. 4.

The original KOSMA- code did not include any angular dependence of the dust scattering but implicitly assumed pure forward scattering (Sternberg & Dalgarno, 1989) in terms of an extinction coefficient . As long as the models use the same value of , the radiative transfer should always give the same scaling of the UV intensity with exp( - ) at high optical depths because the radiation from scattering at other angles is quickly damped because of the longer optical path, but the scattering will increase the intensity close to the surface (Flannery et al., 1980).

To study the influence of FUV scattering we performed a set of model calculations with different scattering properties but the same effective value of  cm:

  1. The old KOSMA- FUV calculations using pure forward scattering.

  2. The MCDRT result for a MRN dust distribution providing the same , but where the albedo is artificially set to zero.

  3. The MCDRT result for the MRN distribution with a full treatment of scattering and absorption.

In Appendix A we discuss how the assumption of isotropic scattering compares to more realistic scattering properties and demonstrate that it poses a clear improvement compared to the pure forward scattering case used in our previous model.

Figure 4 shows a comparison of the depth-dependant mean FUV intensity scaled by the unattenuated mean FUV intensity for these three models. is the total mean intensity, i.e., averaged over the full solid angle and integrated over the full wavelength range: is the corresponding total mean intensity of the FUV field in the absence of the molecular cloud, i.e., the full unattenuated photon flux coming from solid angle. For , models 1) and 2) are slightly larger than 0.5, which would be the theoretical value for a fully opaque and very large clump. Higher values indicate the contribution from angles slightly larger than 90. By contrast model 3 shows a value of 0.58, indicating significant contributions to by scattered photons. Consequently, the mean intensity is higher for model 3) throughout the whole clump because of a large amount of scattered photons. The small differences between models 1) and 2) result from the full wavelength treatment compared to the attenuation by an average only. For high optical depths models 2) and 3) converge to the same intensities as all sidewards scattered photons turn insignificant there because of the longer optical path (Flannery et al., 1980).

Eq. (12) turns out to be a very good approximation when only forward-scaterring occurs, which is clearly reproduced by the detailed radiative transfer computations with the albedo set to zero. However, only the full radiative transfer treatment can produce the full spatial and spectral knowledge of the radiation field and the proper accounting for scattered photons.

Figure 4: Mean intensity of the UV radiation vs. cloud depth. The black line corresponds to the old KOSMA- result for pure absorption. The blue line corresponds to the MCDRT result for an MRN distribution with full treatment of scattering and absorption. The dashed line is the result for an MRN distribution with the albedo artificially set to zero.
Figure 5: Mean intensity of the UV radiation scaled to the unshielded (empty space) Draine field vs. cloud depth. Different lines denote different dust compositions.

The influence of dust composition on the resulting depth-dependent FUV intensity in the clump can be seen in Fig. 5. For the FUV intensity in the WD01 models, the following is true throughout the whole clump: WD01-7 < WD01-21 < WD01-25.

3.1.4 Influence on photo-reactions

The local UV intensity acts on the chemistry of the clump by heating the gas and dust and by inducing ionization and dissociation reactions. This photo-process with a wavelength-dependent cross section proceeds at a rate (Roberge et al., 1991)


is the mean intensity of radiation at radius in photons cms nm. The integral runs from  Å to the threshold wavelength for the process 444A database of photo-ionization and photo-dissociation cross sections for common astrophysical species can be found at http://home.strw.leidenuniv.nl/ewine/photo/. For more details see van Dishoeck et al. (2006) and van Hemert & van Dishoeck (2008). Astrochemical databases often provide parametrized reaction rate coefficients to allow the calculation of these photo-reaction rates. For instance, UDfA gives the rate coefficients in the form


assuming an attenuation of the radiation provided by a standard MRN dust model in a plane-parallel configuration using the radiation transfer results from Flannery et al. (1980) parametrized in terms of the perpendicular visual extinction . is the FUV field strength at the edge of the cloud and the unshielded, free-space rate coefficient assumes a mean FUV intensity of unity (in units of the Draine field) and FUV photons coming from all directions. At the edge of an optically thick molecular cloud the rate coefficient is about half of this value, depending on the dust scattering properties.

Figure 6: Fit results for for \ceHCO (black), \ceCH2 (red), \ceCN (green), and \ceN2 (blue) for all dust models from Weingartner & Draine (2001b, their Tab. 1). The abscissa gives the number of the dust model. The large points show the original values of for the respective molecules.
Figure 7: Comparison of explicitly fitted for \ceHCO (black), \ceCH2 (red), \ceCN (green), and \ceN2 (blue) for all dust models from Weingartner & Draine (2001b, their Tab. 1) with scaling relationship Eq. 3.1.4. The large points show the original values of from UDfA, the open diamonds show the values calculated using an MRN dust model, and the open squares denote the explicitly calculated values assuming .

The unshielded rate coefficient does not depend on the dust properties and we concentrated on the attenuation properties, parametrized in form of the . Changing either the dust properties, i.e., the wavelength-dependant FUV attenuation, or the spectral distribution of the FUV field leads to a different result for and consequently to a different fit parameter . However, performing the wavelength-dependent integration of Eq. 13 for every reaction and every radial point is often beyond the scope of PDR models, e.g., because it is too time-consuming or because neither dust properties nor ionization cross sections are known accurately enough to justify the computational effort of a full integration of the wavelength-dependent radiative transfer. Furthermore, knowledge of the wavelength-dependent is not available for all astrophysically relevant species and models have to apply Eq. 14 even if they can perform the full integration in principle.

Here we derive a scaling relation between the values provided by a data base like UDfA and corresponding fit values for a particular dust model to calculate the rate coefficient for species and dust sort . For a direct comparison with UDfA values we started from a model cloud with a mass of  M and a radius of 1.71 pc, large enough to neglect the effects of the spherical symmetry close to the surface, isotropically illuminated and calculated the photo-dissociation and photo-ionization rates of 72 species provided by van Dishoeck et al. (2006); van Hemert & van Dishoeck (2008) for all 25 dust models from Weingartner & Draine (2001b) according to Eq. 13.

To this end, we performed least-squares fits to up to an to determine new coefficients for all and . Figure 6 shows the new for four example molecules, \ceHCO, \ceCH2, \ceCN, and \ceN2 for all 25 dust models.

The big dots on the left show the original value from UDfA555van Dishoeck et al. (2006) and van Hemert & van Dishoeck (2008) provide their own fitted values of (up to an maximum ) for all species for which they also give tables of . . The figure shows that the new show the same behavior along the abscissa for all four example species. There are also three branches of , visible as prominent quantitative steps in Fig. 6 that correspond to the respective values.

The next step is to identify what dust model property is best correlated with . An extensive analysis of the dust properties reveals that we can find a wavelength for which


shows a monotonic ordering for the of the 25 dust models. Consequently, we chose as the parameter describing the dust. Table 1 gives for all 25 WD01 dust models.

Dust model , see Tab. 1 in WD01 , ratio of visual extinction to reddening. C abundance in double log-normal very small grain population. CaseCase A: variable grain volumes. Case B: grain volumes fixed at approximately the values found for . orderOrdering of the various dust models from lowest to highest value of .
1 3.1 0.0 A 4.13 19
2 3.1 1.0 A 4.31 20
3 3.1 2.0 A 4.44 21
4 3.1 3.0 A 4.56 22
5 3.1 4.0 A 5.02 23
6 3.1 5.0 A 5.22 24
7 3.1 6.0 A 5.62 25
8 4.0 0.0 A 2.79 9
9 4.0 1.0 A 2.98 10
10 4.0 2.0 A 3.10 12
11 4.0 3.0 A 3.26 15
12 4.0 4.0 A 3.54 17
13 5.5 0.0 A 1.66 1
14 5.5 1.0 A 1.79 2
15 5.5 2.0 A 1.97 6
16 5.5 3.0 A 2.21 8
17 4.0 0.0 B 3.01 11
18 4.0 1.0 B 3.11 13
19 4.0 2.0 B 3.16 14
20 4.0 3.0 B 3.30 16
21 4.0 4.0 B 3.58 18
22 5.5 0.0 B 1.80 3
23 5.5 1.0 B 1.88 4
24 5.5 2.0 B 1.96 5
25 5.5 3.0 B 2.10 7

Table 1: Dust model parametrization of all WD01 dust models. See Weingartner & Draine (2001b) for more details.

Using this parameter, we can approximate the dependence of on and the dust model by

fit parameter
Table 2: Fit parameters in Eq. 3.1.4.

Equation 3.1.4 reproduces the explicitly calculated attenuation of the rate coefficients within an accuracy of , and a maximum absolute residual of 0.22. This is comparable to the uncertainties of due to fitting up to a different maximum . Note that for the fitting we removed seven species (photo-dissociation of \ceCH+, \ceSH+, \ceOH+, \ceHCO+, \ceCO, \ceO2+, and \ceSiO ) from the data set because their represent strong outliers with respect to the remainder of the data set. We discuss possible reasons for these deviations in Appendix B. To apply Eq. 3.1.4 to them, we had to use corrected coefficients instead of the UDfA vales when computing . The values are given in Tab. 3.

\ceCH+ 2.5 2.11
\ceSH+ 1.8 1.39
\ceOH+ 1.8 2.76
\ceHCO+ 2. 3.036
\ceCO 2.5 3.03
\ceO2+ 2. 1.61
\ceSiO 2.3 1.87
Table 3: Corrected for the outlying species.

It turns out that the original MRN distribution with does not exactly follow the parametrization of Eq. 3.1.4, based on the WD01 dust models, but that it can be easily used in the same function when adjusting the parameter . Figure 7 compares for \ceHCO, \ceCH2, \ceCN, and \ceN2 with the results of Eq. 3.1.4. In addition we show the MRN results in the plot: the original , the explicitly calculated , and the assuming using large points, open diamonds, and open squares, respectively (the arrows in the figure demonstrate the shifting of the rescaled closer to Eq. 3.1.4).

It is unclear why the from UDfA, i.e., the MRN values, do not follow the parametrization of Eq. 3.1.4. There is a number of possible reasons for the deviation. Woodall et al. (2007) did not use a Draine FUV radiation field but assumed a 10000 K black-body radiation field. The different spectral shape usually accounts for approximately 10% difference but can in a few special cases result in a significant difference. Secondly, different up to which the is fitted give different values of ( see Appendix B for a more detailed discussion). The details of the radiation scattering is another possible reason for different results. Assuming different anisotropy factors and dust albedo will change the radiation field and accordingly affect the (see Appendix A for a discussion on the effect of different values of ).

We conclude that the scaling relation Eq. 3.1.4 is able to rescale the tabulated from UDfA to the dust-specific (fitted to the fully calculated ) with an accuracy of about 10%.

3.2 Dust temperatures

Dust temperatures are calculated independently for each dust component and size. Figures 8, 9, and 10 show depth-dependent grain temperatures for a number of differently sized grains for three of the four components. The plot for the ionized PAHs is omitted because their temperature behaves like that of neutral PAHs. The general behavior is that smaller grains exhibit higher grain temperatures.

Figure 8: Dust temperatures of silicate grains for a model clump of 10 M,  cm, and using the WD01-7 dust model. The different lines denote different grain sizes.
Figure 9: Dust temperatures of carbonaceous grains for a model clump of 10 M,  cm, and using the WD01-7 dust model. The different lines denote different grain sizes.
Figure 10: Dust temperatures of neutral PAHs for a model clump of 10 M,  cm, and using the WD01-7 dust model. The different lines denote different PAH sizes.

Figure 11 compares grain temperatures of the components of the MRN and WD01-07 dust model at the surface of a model clump as a function of grain radius . Graphites tend to have a higher temperature than silicate grains except for very large grain sizes. Small grains tend to be hotter than bigger grains. The PAHs are heated much more efficiently than the larger grains and have significantly higher grain temperatures. For simplicity, we deliberately neglected the non-equilibrium heating of very small grains and PAHs here and assumed equilibrium heating in this paper. In a subsequent paper, we will include the comparison with real observations accounting for the stochastic heating of the very small particles. The dependence on vanishes for sufficiently high extinctions i.e., once most of the remaining photons have wavelengths longer than the grains. The resulting equilibrium grain temperature then only depends on the absorption coefficient of the material.

Figure 11: Grain temperature for the components of the MRN and WD01-07 dust models as a function of the radius in nm. The grain temperatures are calculated at the surface of a model clump of 10 M,  cm, and . The different lines denote different dust components. Dashed lines show the MRN components and solid lines show the WD01-07 components.

The detailed knowledge of the dust temperature as a function of radius, grain size, and grain type enables us to treat all dust-related processes separately for all grain sizes and types. However, since most relevant processes are grain surface reactions, e.g., formation of the H or photoelectric heating, it is usually sufficient to calculate surface-averaged quantities. For the -th dust component we can define the relevant mean dust temperature as


Averaging over all dust sorts is done according to


with the geometrical cross section of the -th sort and the total dust cross section . We will use the short-hand notation and from here on.

Figure 12 shows the mean dust temperature for the WD-7 dust model as a function of cloud depth for all four dust components. As at the surface (Fig. 11), silicates, indicated by the solid blue line, are cooler than graphite grains (solid, red line). Graphite and silicate grain temperatures are very close, but the mean silicate temperature drops slightly faster for increased cloud depths than the graphites. The PAHs have a much higher . Neutral PAHs tend to be marginally hotter than the ionized PAHs. The black line denotes the average dust temperature for the WD01-7 dust.

The influence of the different dust models can be seen from Fig. 13. The black line shows the old dust temperature calculation in KOSMA- following Hollenbach et al. (1991), which gives a central dust temperatures  K. Accounting for the full wavelength dependence and the detailed grain size distribution gives a significantly lower  K at the center of the model clump. On the other hand, it also leads to warmer dust at the surface of the cloud compared to Hollenbach et al. (1991).

The dust temperature for the WD01 models follow their corresponding FUV intensities (see Fig. 5). The MRN model has a much lower . Figure 11 shows that the for MRN and WD01-7 are very similar for equal grain sizes. However, the temperature of the PAHs in the WD01-7 model contributes significantly to the average dust temperatures in Eq. (18) and leads to higher in models with PAH content.

The lower central dust temperatures of models with WD01 dust types also influences the gas temperatures. Models with WD01 dust show significantly lower central gas temperatures because gas-grain collisions are becoming the dominant cooling process.

Figure 12: Mean dust temperature for the WD01-7 dust model for a model clump of 10 M,  cm, and . The different lines denote the average temperature of the individual dust components. The black line shows the overall mean dust temperature.
Figure 13: Mean dust temperature as a function of  for a model clump of 10 M,  cm, and . The different lines denote different dust compositions.

3.3 H formation

model parameter silicates graphite PAH PAH total

[Å] 50 50
0.25 0.25
[%] 0.598 0.357 0.955

[Å] 3.5 100 3.5 3.5
1. 10. 0.01 0.01
[%] 0.828 0.0.228 0.038 0.038 1.130

[Å] 3.5 100 3.5 3.5
1. 10. 0.01 0.01
[%] 0.815 0.0.252 0.027 0.027 1.120

[Å] 3.5 100 3.5 3.5
1. 10. 0.01 0.01
[%] 0.812 0.0.268 0.019 0.019 1.118
Table 4: Parameter of the WD01 and MRN dust models

Molecular hydrogen is the most abundant molecule in the universe. The formation of H most effectively takes place on the surface of dust grains. The proper treatment of the H formation is crucial in any numerical PDR model (see e.g. Le Bourlot et al., 2012, for a recent update of the Meudon PDR code.). We updated the KOSMA- calculation of the H formation efficiency from the method described by SD89 to account for physisorption and chemisorption binding energies on silicate and graphite surfaces (Cazaux & Tielens, 2002, 2004, 2010, (hereafter CT)).

This model is based on two main points: (1) H atoms can stick to the grain surface in two different binding sites, a physisorption site and a chemisorption site. Physisorption is a weak atom-surface binding resulting from a Van der Waals force (dipole-dipole interaction). Chemisorption is a strong atom-surface binding interaction involving the valence electrons, known as covalent bond. These interactions determine whether atoms move on the surface, and thus, how quickly they can meet a recombination partner. (2) The mobility of a surface H atom results from the combination of two physical processes: tunneling and thermal diffusion; tunneling dominates at the lowest temperatures while thermal diffusion is most important at the highest temperatures. The calculation of the transmission coefficients, to go from site to site, are given in Cazaux & Tielens (2004, 2010).

At low temperature ( K), H formation involves the migration of physisorbed H atoms. At higher temperatures ( K), H formation results from chemisorbed H recombination. The presence of these two types of binding sites allows Hformation to proceed relatively efficiently at low and elevated temperatures.

For the dust component , the formation rate can be written as


where and are the number density and thermal velocity of H atoms in the gas phase, is the total cross section of the grain component , is the formation efficiency (Cazaux & Tielens, 2004) on surfaces of the -th component, and is the sticking coefficient of the H atoms


which depends on the gas and dust temperature, and (Hollenbach & McKee, 1979). The thermal velocity is  cm s. The -th dust cross section per H nucleus is given by


Table 4 lists for the different dust models. Cazaux & Tielens (2002) expressed the general expression for the -th H formation efficiency:

Figure 14: H formation efficiency for carbon and silicate surfaces. and denote the temperature regimes of the corresponding terms in Eq. (22). The dash-dotted line denotes the standard formation efficiency given by Hollenbach & McKee (1979).

is limited by three terms: the first term prohibits the evaporation of the newly formed molecules at low temperatures. The second term, equal to unity, dominates at higher grain temperatures; all incoming H atoms leave the surface as H. The third term governs the high-temperature regime, where evaporation of physisorbed atoms removes \ceH from the surface before it can recombine to H. It competes with the hopping of atoms into chemisorbed binding sites where the atoms can remain on the surface to form H at significantly higher temperatures. The numerator , describing the chemisorption, terminates the H formation at much higher temperatures, above several hundred Kelvin, when even chemisorbed atoms start to evaporate.

Figure 14 shows as a function of for silicate (blue line) and graphite (red line) surfaces. The old formation efficiency from Hollenbach & McKee (1979) is shown as a dash-dotted line for comparison. When applying this formation in the framework of the PDR model, see Fig. 14, the term turned out to be problematic because it shuts down the H release into the gas phase once K, accumulating all molecules on the grain surfaces. We question this term because Eq. (22) does not account for the H binding energy of 4.48 eV that upon formation is available to allow release into the gas phase even at very low temperatures. As a remedy we set to ensure efficient H formation even on cold dust surfaces (Cazaux (2011), priv. comm.)777All model results accounting for CT H formation shown in this paper assume .. More details on how to calculate the H formation efficiency are given in Appendix  C. The total H formation rate can then be written as

Figure 15: H formation rate as a function of dust temperature for different formation models and dust compositions. The gas temperature is set to 50 K. The dashed lines indicate the-low temperature formation in the original CT formation rates with . SD89 is the formation rate from Sternberg & Dalgarno (1989).
Figure 16: H formation rate as a function of dust temperature for different dust compositions. The gas temperature is set to 50 K. The dashed lines show the formation rate with suppressed H formation on PAHs. The solid lines show the formation rate allowing for H formation on PAH surfaces.

Fig. 15 compares the new and old H formation rates. The formation efficiency given by SD89 follows the treatment of Hollenbach & McKee (1979), which adopted a critical dust temperature of  K above which the H formation efficiency drops below 0.5. H formation becomes inefficient above  K. The H formation in the CT dust model remains efficient up to 1000 K. The model by Cazaux&Tielens allows a much more efficient H formation at higher temperatures. Additionally, the significantly larger dust surface in the MRN and WD01-7 model lead to a much higher H formation rate than SD89.

While accounting for chemisorbed H atoms significantly increases the H formation efficiency at higher surface temperatures, it is the total available grain cross section that distinguishes between different dust models. Table 4 also compares the total dust cross section per H nucleus for the different dust compositions. The PAH cross sections contribute the majority to the final H formation rate. The resulting H formation rate for one H atom cm is shown in Fig. 16. Solid lines show the formation rate as a function of for a constant  K for the WD01-7, WD01-21, and WD01-25 models allowing for H formation on the surface of PAHs. The dashed lines show the formation rates if the PAH surfaces are excluded from the H formation calculation. Consequently, the dust model with the largest effective cross section, WD01-7 shows the highest formation rate in both cases.

Figure 17: H formation heating rate for a model clump of 10 M,  cm, and . The different lines show the effect of different H formation treatment. All models have been computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating.

The strong increase in the total H formation rate has an important effect in addition to the easier H production. Each formed H molecule releases its binding energy of 4.48 eV. A common assumption (e.g., Habart et al., 2004) is that this energy is equally used to a) release the molecule into the gas phase, b) increase the kinetic energy of the molecule, and for c) internal rotational-vibrational excitation. If we assume a uniform distribution, eV are available to heat the gas. Here we can measure the H formation rate across the cloud by equivalently looking at the H formation heating.

By including H formation according to Eq. (23) we also increase the formation rate significantly at low values of . This is shown in Fig. 17. All models in the plot were computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating. The solid line shows the results for the formation rate from Sternberg & Dalgarno (1989). The result for the full H formation treatment from CT (with ), suppressing and allowing H formation on PAH surfaces, are shown by the dotted and dashed line. The much larger surface area of the WD01-7 dust model leads to a significant increase of the formation rate. Particularly the large PAH surface in that dust model contributes dominantly to the total grain surface and consequently to the total H formation heating.

It is important to note that the H formation heating is extensive at cloud depths where most of the hydrogen is in the form of \ceH. This leads to the curious situation of a strong H formation heating in the absence of H molecules. The gas temperatures at these depths can reach more than thousand Kelvin. Under these conditions two more H destruction reaction become important:

(C 1)
(C 2)

Each destruction process requires 4.48 eV to break the binding, in contrast to the release of 4.48 eV during the formation. The much more efficient formation of molecular hydrogen by Cazaux & Tielens makes it necessary to adopt an additional cooling term in the chemical network, accounting for the 4.48 eV energy consumption during the kinetic dissociation (Lepp & Shull, 1983). The cooling rate can be written as


where and are the temperature-dependent rate coefficients for reactions (C 1) and (C 2).

The role of interstellar dust in ISM physics and chemistry remains one of the key problems in modern astrophysics. In this section we demonstrate the strong effect of H formation on PAH surfaces on the chemistry and the thermal conditions in PDRs. However, it remains unknown whether PAHs really contribute to the formation of molecular hydrogen. Naturally, this is connected to the uncertainty about the detailed form and distribution of interstellar PAHs.

By assuming that \ceH2 formation can take place on PAH surfaces the formation rate from Cazaux & Tielens predicts rates of a few times  s when applied to the dust distributions presented by WD01. This is about a factor 10 higher than what has commonly been found for the diffuse medium (Jura, 1974; Browning et al., 2003; Gry et al., 2002; Welty et al., 2003; Wolfire et al., 2008). However, the usual values of  cm s are mean values along the line of sight while our results are local values and a direct comparison is difficult. Recently, Le Bourlot et al. (2012) presented an update to the \ceH2 formation formalism in the Meudon PDR code and also found a strong increase of the local \ceH2 formation rates. If we prevent \ceH2 formation on PAH surfaces, the total grain surface available to form molecular hydrogen is significantly reduced and we predict formation rates of the same order as found in the diffuse ISM. However, PDR models assuming diffuse gas formation rates are not able to reproduce the observed extent of \ceH2 line excitation (Habart et al., 2011). \ceH2 formation on PAHs has also been proposed by Habart et al. (2003) to explain the \ceH2 emission in  Oph, and recently Mennella et al. (2012) concluded from experimental studies that hydrogenated neutral polycyclic aromatic hydrocarbon molecules act as catalysts for the formation of molecular hydrogen.

3.4 Photoelectric heating

Weingartner & Draine (2001a) presented the photoelectric heating rate for the WD01 grain size distributions in a parametrized form as a function of , where is the gas temperature in Kelvin and is in units of Kcm. To convert the FUV field from Habing to Draine units use . These parametrizations are fairly accurate when and . However, the parameter range in PDR models allows to reach values as high as and as low as . Weingartner (2009) provided updated calculations of and for the dust distributions WD01-7, WD01-21, and WD01-25 up to . These rates are fairly well reproduced by the new analytic fits:



Numerical values for and are tabulated in the appendix in Table 8 and 9.

Figure 18: Photoelectric heating rate for a model clump of 10 M,  cm, and . The different lines denote different dust compositions. Models with an MRN dust composition apply the PE heating efficiency as given by Bakes & Tielens (1994). Dotted lines indicate a shift from heating to cooling because of .

Figure 18 shows the PE heating rates for a model clump of 10 M,  cm, and . The solid lines show the photoelectric heating (PEH) rates using the method described by Bakes & Tielens (1994). The PEH rate obtained for an MRN dust composition when solving the full FUV radiative transfer problem is enhanced compared to the old KOSMA- values because of the stronger FUV intensities due to the scattering. The dashed lines in Fig. 18 show the heating rates using the updated parametrizations from Eqs. (25) and (3.4) for the different dust models. Dotted lines indicate negative values, i.e., a transition from heating to cooling because . For the PEH rate of the WD01 dust models follow the FUV intensity of the respective models, i.e., the model with the highest mean FUV intensity, WD01-25, shows the strongest PEH rate while the model WD01-07 drops quickest. However, it is interesting to note that for the PEH rate for the WD01-07 dust produces the strongest PEH rate, leading to the highest gas temperature of all models (see Fig. 19).

4 Application of the single clump model

4.1 Impact on gas temperature

Figure 19: Gas temperature for a model clump of 10 M,  cm, and . Left panel: The different lines denote different dust models that affect the FUV radiative transfer and the PE heating. All models in the left panel werecalculated using the SD89 H formation. Right panel: The different lines show the effect of different H formation treatment. All models in the right panel were computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating.

The changes in computating the PE heating and the H formation, which depend on the applied dust model, have a profound effect on the thermal structure of the model clump. Figure 19 shows the separate effects of the changed FUV and resulting PE heating (left panel) and different H formation treatment (right panel). Among the different dust models presented by Weingartner & Draine (2001b) the WD01-7 dust models shows the highest total PE heating rate because it contains the largest population of very small grains (Weingartner & Draine, 2001a). Depending on , the PE heating rate from Bakes & Tielens (1994) can become higher than in the WD01 models despite its smaller population of very small grains. This is because they assume a higher electron sticking coefficient. In the chemically very active zone of , the new FUV treatment leads to systematically warmer gas temperatures. The net effect of the dust model on the PE heating strongly depends on the applied dust model, with the PE heating efficiency proportional to the size of the very small grain population. A second effect comes from the H formation heating. The right panel in Fig. 19 shows how the increase in the total available grain surface raises the gas temperature. However, this mostly affects the outer layers of a model clump where the high grain temperatures require accounting for chemisorbed H atoms to allow the formation of molecular hydrogen.

Overall, we can distinguish two regions: For optical depths , the different attenuation of the FUV field by different dust distributions dominates the gas temperature. The model with the lowest reddening, i.e., the highest FUV intensity at those depths, shows the highest gas temperature. At lower optical depths, the FUV intensity is in all cases high enough that the total dust surface available for PE heating and H formation determines the gas temperature. Here, the models with PAHs and many small grains, which also produce a higher reddening, show the highest gas temperature.

4.2 Impact on gas chemistry

The described changes to the model have a significant effect on the chemical structure of the cloud. The influence of the assumed dust composition and dust size distribution affects the FUV radiative transfer and the improved treatment of dust scattering increases the FUV intensity throughout the clump compared to the original approximation. Chemical species that are dominantly formed or destroyed via photo-ionization or photo-dissociation processes are affected. Secondly, the stronger H formation efficiency leads to an increased H formation heating contribution, which in turn produces an increase in gas temperature at low . Consequently, the abundance of the species that are dominantly formed in the outer regions of the molecular clump is changed because the higher gas temperature accelerates their respective formation and destruction reactions in the gas phase.

Figure 20: H and H abundances for a model clump of 10 M,  cm, and . The different lines denote different treatments of the H formation. The solid line shows a model using the standard H formation rate as given by Sternberg & Dalgarno (1989). The dashed and dotted lines show model results for H formation according to CT where the formation on PAH surfaces is switched on and off, respectively. All three models assume a WD01-7 dust distribution. The green dash-dotted line assumes an MRN dust model and uses the CT H formation.
Figure 21: CO and C abundances for a model clump of 10 M,  cm, and . Left panel: The different lines denote different dust models that affect the FUV radiative transfer and the PE heating. All models in the left panel have been calculated using the SD89 H formation. Right panel: The different lines show the effect of different H formation treatment. All models in the right panel have been computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating.

Both effects can go in opposite directions, but usually one effect dominates. In Fig. (20) we show the dependence of the H-H transition zone on the applied H formation rate. We compare three models assuming a WD01-7 dust model with full radiative transfer and photoelectric heating treatment. The only difference between the model calculations is the treatment of H formation. We apply (1) the standard H formation rate as given by Sternberg & Dalgarno (1989) (black, solid line), (2) the H formation given by CT with H formation taking place on the PAH surface (blue, dashed line), (3) and the H formation given by CT with H formation on PAHs suppressed (red, dotted line). H chemistry is unique in the sense that it only weakly depends on the gas temperature (see Eqs. (19, 20)). Figure 20 shows how the different treatment of the H formation influences the details of the H-H transition. In that particular model clump, the dust temperature is about 50 K and almost identical in all models, leading to a H formation efficiency of unity across the dust models shown. The only significant difference is the total dust surface of the various models. The SD89 approximation, case (1), has the smallest available surface ( cm) followed by CT without PAHs and CT including PAHs (see Table 4). This agrees with the H-H transition moving outwards with growing dust surface in Fig. 20. To emphasize this behavior, we added a fourth dust model to the plot, showing the results for the MRN model. The total surface in the MRN model is only slightly smaller than the WD01-7 model without PAHs and consequently, both models show a very similar behavior. Please note that the contribution of graphite and silicate to the total dust surface is still very different in the different dust models. The MRN model possesses about equal surfaces in both kinds, while the majority of the big grain surfaces in the WD01-7 model comes from silicates.

In contrast to H , all other chemical species in the gas phase are affected by variations of the local FUV intensity and the local gas temperature. This is apparent, for example, in the C to CO transition. Figure 21 shows the density profile of C and CO for different model calculations. In the left panel we show the influence of the different dust models in terms of radiative transfer and PE heating, while the right panel demonstrates how the different H formation treatment influences the chemistry. All models in the left panel were calculated using the SD89 H formation. The transition from C to CO occurs at different for all models. Consistent with the FUV intensities shown in Figs. 4 and 5, the models with the lowest FUV intensities show a C to CO transition at lowest , while models with a weak dust attenuation, e.g. WD01-25, exhibit the same transition deeper in the cloud. Additionally, we note an enhanced CO density at low in the WD01-7 model. This is a result of the higher gas temperature due to the stronger PE heating efficiencies of that dust model.

Figure 22: Total line integrated clump-averaged intensities of CO transitions for a model clump of M,  cm, and . The different lines show the effect of different H formation treatment. All models were computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating.

The right panel in Fig. 21 shows the isolated effect of the different H formation treatment. All models in the right panel were computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating. The C-CO transition is similar across all shown models, since the FUV intensity, which controlls the photo-dissociation of CO, is the same across all models. However, the CO abundance at lower shows very large differences. The models using the CT dust formation scheme produce more CO at the cloud edge compared to the SD89 H formation. The CT models provide a stronger H formation heating contribution. The additional heating term produces a significant population of hot CO in the outer layers of the model clump.

The production of hot, i.e., strongly excited, CO is visible in the spectral line emission of the model clumps. In Fig. 22 we show the corresponding total CO line emission of the model clumps shown in the right panel in Fig. 21. For transitions higher than J=8-7, the line-integrated intensities start to show significant differences for the three models. In this particular case, the strength of the high- CO transitions is proportional to the total H-forming grain surface. The outer layers of the model clumps contribute significantly to the total CO line emission of this model. As a consequence, any modeling with a large surface of the clouds, e.g., assuming a clumpy structure, will show much stronger high- CO emission lines compared to non-clumpy model attempts.

Figure 23: \ceCH and \ceCH+ abundances (left) and column densities (right) for a model clump of M,  cm, and .The different lines show the effect of different H formation treatment. All models were computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating.
Figure 24: \ceOH and \ceOH+ abundances (left) and column densities (right) for a model clump of M,  cm, and .The different lines show the effect of different H formation treatment. All models were computed using the WD01-7 dust properties regarding the radiative transfer and the PE heating.

As an additional example of chemical species that predominantly form in the outer regions of molecular clouds and thus are affected by the H formation treatment we show in Fig. 23 the \ceCH and \ceCH+ abundances and column densities and in Fig. 24 the \ceOH and \ceOH+ abundances and column densities. All these species are significantly affected by the strong H formation heating effect in the outer parts of the model clump. The column density effect is weakest for the CH because of the abundance peak at . The other three species show a stronger effect on the column density because of the weak formation in deeper clump regions. Recent Herschel observations show high column densities of these surface tracers (e.g. Qin et al., 2010; Falgarone et al., 2010), which were difficult to reproduce in existing PDR models. A more efficient H formation modeling with the accompanying formation heating effect may help to resolve the discrepancy between models and observations.

4.2.1 Full model SED

Figure 25: Full continuum and spectral line flux emitted by a model clump of M,  cm, and . The two panels show the results for different dust models. Top: WD01-7 dust model, equivalent to a . This is the dust model with the largest grain surface contributing to the heating. Bottom: MRN dust mode. This is a dust model with a slightly smaller grain surface compared to WD01-7. The PE heating was calculated according to Bakes & Tielens (1994).
line WD01-7 WD01-7 WD01-25 WD01-25 MRN
w/o PAHs with PAHS w/o PAHS with PAHS
[m] [erg/s/cm] [erg/s/cm] [erg/s/cm] [erg/s/cm] [erg/s/cm]
0.290 0.290 0.286 0.280 0.278
[C ii] 158
[O i] 63
[O i] 146
[C i] 609
[C i] 370
CO (1–0) 2601
CO (4–3) 650
CO (7–6) 372
CO (10–9) 260
CO (15–14) 174
CO (20–19) 130