Issue 
A&A
Volume 549, January 2013



Article Number  A85  
Number of page(s)  21  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201118190  
Published online  21 December 2012 
Full SED fitting with the KOSMAτ PDR code
I. Dust modelling
^{1} I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany
email: roellig@ph1.unikoeln.de
^{2} N. Copernicus Astronomical Center, Rabianska 8, 87100 Toruń, Poland
Received: 30 September 2011
Accepted: 17 October 2012
Aims. We revised the treatment of interstellar dust in the KOSMAτ photodissociation region (PDR) model code to achieve a consistent description of the dustrelated 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 multi component dust radiative transfer (MCDRT) code to to solve the frequencydependent 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 nonequilibrium 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_{2} formation on grain surfaces to include the EleyRideal effect, thus allowing for hightemperature H_{2} 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_{2} formation on grain surfaces that allows for chemisorption. This increases the total H_{2} formation significantly and the connected H_{2} 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_{2} formation provides a boost to the temperature of outer cloud layers, which is clearly traced by highJ 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.
Key words: astrochemistry / radiative transfer / methods: numerical / photondominated region (PDR) / ISM: molecules / infrared: ISM
© ESO, 2013
1. Introduction
Photodissociation 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 farultraviolet (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 timeconsuming 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 emission^{1}.
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 ^{13}C and ^{18}O. 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. Multicomponent dust radiative transfer model
The multicomponent dust radiative transfer (MCDRT) code allows one to solve the frequencydependent 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) i that may or may not be cospatial, and obey any dust size distribution. The grain size distribution is given by the general relation valid for grains with radius a_{i, − } ≤ a ≤ a_{i, + }; (1)where: n_{i}(a) is the number density of grains with size ≤ a and n is the number density of H nuclei (in both atoms and molecules: n = n_{H} + 2n_{H2}). The function f_{i}(a) 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, J_{λ}(r) and the dust temperature, T_{d,i}(r,a), which for now is computed from the assumed thermal equilibrium for each dust species and dust grain size at each radius r of the clump. These two quantities (mean radiation field and dust temperature) are now used in the KOSMAτ PDR model code to obtain a selfconsistent gasdust 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 (2)where i stands for sil and gra, A_{sil} = 10^{25.10}, and A_{gra} = 10^{25.13} (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^{0}) and ionized (PAH^{+}) particles. The grain size distributions for these dust components follow the general definition given by Eq. (1) and their functional forms (f_{i}(a) in Eq. (1)) are given by Eqs. (2)–(6) of WD01. The authors employ a powerlaw form for the large grains and a lognormal distribution for the very small grain population. Free parameters of the f_{i}(a) functions were determined by the best fit to the average extinction with R_{V} = A_{V}/E(B − V)) of 3.1, 4 and 5.5 (see WD01 for details). A_{V} is the absolute visual extinction at V = 5500 Å and E(B − V) = A_{B} − A_{V}. The three dust models from WD01 correspond to lines 7, 21, and 25 from Table 1 in Weingartner & Draine (2001b) and are denoted WD017 (R_{V} = 3.1), WD0121 (R_{V} = 4.0), and WD0125 (R_{V} = 5.5).
To compute the extinction efficiency Q_{ext} 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 Table 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 log (a). 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 M = 10 M_{⊙}, the total surface gas density n = 10^{5} cm^{3}, and the total FUV field χ = 1000 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 R_{tot}/5. This implies R_{tot} = 2.46 × 10^{17} cm = 0.08 pc. We assumed a total cosmic ray ionization rate of molecular hydrogen ζ_{H2} = 5 × 10^{17} s^{1}. The chemistry is based on the UMIST 2006 database for astrochemistry UDfA^{2} (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 τ_{FUV} = σ_{d,FUV}N_{H} (Sternberg & Dalgarno 1989), with σ_{d,FUV} being the effective dust grain extinction cross section in the FUV, and the total hydrogen column density N_{H}. 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 dusttoH ratio the total extinction cross section per H nucleus can be derived from (3)where is the extinction cross section of a single dust particle with size a at wavelength λ. Often the absolute extinction in magnitudes is used (4)By introducing the average mass extinction coefficient for each dust component, i.e., the mean extinction cross section per dust mass, (5)where ρ_{b,i} is the bulk dust grain density (i.e., the density of material of which dust of given sort is composed) and ρ_{dust,i} is volume density of the ith dust sort, the optical depth at any given point of the model clump can be written as (6)Relating this to the gas density provides (7)where n m_{H} = ρ_{gas}, and m_{H} is mass of a single H atom. Using Ψ = ρ_{dust}/(n m_{H}) as the total dusttoH mass ratio, and we can write (8)with (9)For a constant dusttoH mass ratio the integration over r in Eq. (8) can be performed giving an H column density N_{H} using (10)i.e., σ_{ext}(λ) = ⟨ K_{ext}(λ) ⟩ m_{H}Ψ.
Fig. 1
Extinction cross section per H nucleus for the WD017 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. 
3.1.2. Influence of the dust composition
Figure 1 shows the total extinction cross section per H nucleus computed for the WD017 (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 WD017 model (silicates – dashed, red in the electronic version; graphite – dashed, blue in the electronic version; and PAH^{0}+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^{0} 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.
Fig. 2
Relative extinction for the dust models of MRN (solid) , WD017 (longdashed; red in electronic version); WD0121 (shortdashed; green in the electronic version) and WD0125 (dotdotdashed; 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 k_{FUV} for the considered models. 
Figure 2 shows A_{λ}/A_{V} for the dust models tested in this paper: MRN (solid), WD017 (longdashed; red in the electronic version), WD0121 (shortdashed; green in the electronic version), and WD0125 (dotdotdashed; blue in the electronic version). The FUV range relevant for photoelectric heating and photodissociation reactions extends between 13.6 and 6 eV^{3}. Consequently, we define an FUVtoV color as k_{FUV} = ⟨ A(λ)/A(V) ⟩ _{λ} ≡ ⟨ A_{λ}/A_{V} ⟩ _{λ} where the averaging is performed over an energy from 6 to 13.6 eV and A_{V} is the visual extinction. The vertical dotted lines in Fig. 2 indicate this range of averaging and the horizontal lines show k_{FUV} for the considered models. Note that the MRN and WD017 dust models have the same value of k_{FUV} = 3.339. All values of the average relative dust extinction are collected in Table 4. From Fig. 2 it can be seen that k_{FUV} is significantly decreased in the WD0121 and WD0125 dust models. This is consistent with the higher values of R_{V}, i.e., a shallower gradient of A(λ)/A_{V} toward shorter wavelengths causing a lower UV extinction for the same A_{V} (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
As a result of the new radiative transfer computations we obtain the full, wavelengthdependent FUV radiation field at all clump radii, where I_{λ} is the specific intensity at wavelength λ. In Fig. 3 we plot J_{λ}(r) for radii in steps of R_{tot}/10 for the WD017 model. The effect of the prominent 2200 Å “bump” is plainly visible in the enhanced reduction of J_{λ} with decreasing radius.
The FUV optical depth (11)can be described by the effective dust extinction cross section (12)For adopted dust models, the obtained values of σ_{D,FUV} are collected in Table 4.
Fig. 3
Mean intensity J_{λ}(r) for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The lines show J_{λ}(r) for radii in steps of 10% of R_{tot}. 
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 σ_{D,FUV}. As long as the models use the same value of σ_{D,FUV}, the radiative transfer should always give the same scaling of the UV intensity with exp( − N_{H}σ_{D,FUV}) 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 σ_{D,FUV} = 1.76 × 10^{21} cm^{2}:

1.
The old KOSMAτ FUV calculations using pureforward scattering.

2.
The MCDRT result for a MRN dust distribution providing the same σ_{D,FUV}, 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 depthdependant mean FUV intensity scaled by the unattenuated mean FUV intensity J(r)/χ_{0} for these three models. J(r) is the total mean intensity, i.e., averaged over the full solid angle and integrated over the full wavelength range: . χ_{0} 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 4π solid angle. For r = R_{tot}, 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 J(r) 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 τ_{FUV} 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).
Equation (12) turns out to be a very good approximation when only forwardscaterring 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.
Fig. 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. 
Fig. 5
Mean intensity of the UV radiation J(r) scaled to the unshielded (empty space) Draine field χ_{0} vs. cloud depth. Different lines denote different dust compositions. 
The influence of dust composition on the resulting depthdependent 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: WD017 < WD0121 < WD0125.
3.1.4. Influence on photoreactions
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 photoprocess j with a wavelengthdependent cross section σ_{j}(λ) proceeds at a rate (Roberge et al. 1991) (13)J_{λ}(r) is the mean intensity of radiation at radius r in photons cm^{2}s^{1} nm^{1}. The integral runs from λ_{H} = 912 Å to the threshold wavelength λ_{j} for the process j^{4}. Astrochemical databases often provide parametrized reaction rate coefficients to allow the calculation of these photoreaction rates. For instance, UDfA gives the rate coefficients in the form (14)assuming an attenuation of the radiation provided by a standard MRN dust model in a planeparallel configuration using the radiation transfer results from Flannery et al. (1980) parametrized in terms of the perpendicular visual extinction A_{V}. χ_{0} is the FUV field strength at the edge of the cloud and the unshielded, freespace rate coefficient α_{j} 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.
The unshielded rate coefficient α_{j} 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 wavelengthdependant FUV attenuation, or the spectral distribution of the FUV field leads to a different result for Γ_{j} and consequently to a different fit parameter γ_{j}. However, performing the wavelengthdependent 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 timeconsuming or because neither dust properties nor ionization cross sections σ_{j}(λ) are known accurately enough to justify the computational effort of a full integration of the wavelengthdependent radiative transfer. Furthermore, knowledge of the wavelengthdependent σ_{j} 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 γ_{j} values provided by a data base like UDfA and corresponding fit values for a particular dust model to calculate the rate coefficient for species j and dust sort D. For a direct comparison with UDfA values we started from a model cloud with a mass of 10^{3} 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 photodissociation and photoionization rates of 72 species provided by van Dishoeck et al. (2006); van Hemert & van Dishoeck (2008) for all 25 dust models D from Weingartner & Draine (2001b) according to Eq. (13).
Fig. 6
Fit results for for HCO (black), CH_{2} (red), CN (green), and N_{2} (blue) for all dust models from Weingartner & Draine (2001b,their Table 1). The abscissa gives the number of the dust model. The large points show the original values of γ_{j} for the respective molecules. 
Dust model parametrization of all WD01 dust models.
To this end, we performed leastsquares fits to up to an A_{V} = 10 to determine new coefficients for all D and j. Figure 6 shows the new for four example molecules, HCO, CH_{2}, CN, and N_{2} for all 25 dust models.
The big dots on the left show the original value γ_{j} from UDfA^{5}. 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 R_{V} 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 λ = 760 Å for which (15)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. Using this parameter, we can approximate the dependence of on γ_{j} and the dust model by (16)Equation (16) reproduces the explicitly calculated attenuation of the rate coefficients within an accuracy of σ = 0.05, and a maximum absolute residual of 0.22. This is comparable to the uncertainties of due to fitting Γ_{j}(A_{V})/Γ_{j}(0) up to a different maximum A_{V}. Note that for the fitting we removed seven species (photodissociation of CH^{+}, SH^{+}, OH^{+}, HCO^{+}, CO, O, and SiO ) 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. (16) to them, we had to use corrected coefficients instead of the UDfA vales when computing . The values are given in Table 3.
Fig. 7
Comparison of explicitly fitted for HCO (black), CH^{2} (red), CN (green), and N_{2} (blue) for all dust models from Weingartner & Draine (2001b,their Table 1) with scaling relationship Eq. (16). The large points show the original values of γ_{j} from UDfA, the open diamonds show the values calculated using an MRN dust model, and the open squares denote the explicitly calculated values assuming . 
Corrected for the outlying species.
It turns out that the original MRN distribution with does not exactly follow the parametrization of Eq. (16), 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 HCO, CH_{2}, CN, and N_{2} with the results of Eq. (16). In addition we show the MRN results in the plot: the original γ_{j}, 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. (16)).
It is unclear why the γ_{j} from UDfA, i.e., the MRN values, do not follow the parametrization of Eq. (16). 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 10 000 K blackbody 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 A_{V,max} up to which the γ_{j} is fitted give different values of γ_{j}( 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 g and dust albedo ω will change the radiation field and accordingly affect the γ_{j} (see Appendix A for a discussion on the effect of different values of g).
We conclude that the scaling relation Eq. (16) is able to rescale the tabulated γ_{j} from UDfA to the dustspecific (fitted to the fully calculated Γ_{j}(A_{V})/Γ_{j}(0)) with an accuracy of about 10%.
3.2. Dust temperatures
Dust temperatures are calculated independently for each dust component and size. Figures 8–10 show depthdependent 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.
Fig. 8
Dust temperatures of silicate grains for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000 using the WD017 dust model. The different lines denote different grain sizes. 
Fig. 9
Dust temperatures of carbonaceous grains for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000 using the WD017 dust model. The different lines denote different grain sizes. 
Fig. 10
Dust temperatures of neutral PAHs for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000 using the WD017 dust model. The different lines denote different PAH sizes. 
Figure 11 compares grain temperatures of the components of the MRN and WD0107 dust model at the surface of a model clump as a function of grain radius a_{i}. 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 nonequilibrium 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 a_{i} 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.
Fig. 11
Grain temperature for the components of the MRN and WD0107 dust models as a function of the radius a_{i} in nm. The grain temperatures are calculated at the surface of a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote different dust components. Dashed lines show the MRN components and solid lines show the WD0107 components. 
The detailed knowledge of the dust temperature as a function of radius, grain size, and grain type enables us to treat all dustrelated processes separately for all grain sizes and types. However, since most relevant processes are grain surface reactions, e.g., formation of the H_{2} or photoelectric heating, it is usually sufficient to calculate surfaceaveraged quantities. For the ith dust component we can define the relevant mean dust temperature as (17)Averaging over all dust sorts is done according to (18)with the geometrical cross section of the ith sort σ_{d,i} and the total dust cross section σ_{d}. We will use the shorthand notation T_{d} = ⟨ T_{d} ⟩ _{NDS} and T_{d,i} = ⟨ T_{d,i} ⟩ _{a} from here on.
Figure 12 shows the mean dust temperature T_{d,i} for the WD7 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 T_{d,i}. Neutral PAHs tend to be marginally hotter than the ionized PAHs. The black line denotes the average dust temperature T_{d} for the WD017 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 T_{d} ≈ 22 K. Accounting for the full wavelength dependence and the detailed grain size distribution gives a significantly lower T_{d} ≈ 13 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 T_{d}. Figure 11 shows that the T_{d,i} for MRN and WD017 are very similar for equal grain sizes. However, the temperature of the PAHs in the WD017 model contributes significantly to the average dust temperatures in Eq. (18) and leads to higher T_{d} 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 gasgrain collisions are becoming the dominant cooling process.
Fig. 12
Mean dust temperature for the WD017 dust model for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote the average temperature of the individual dust components. The black line shows the overall mean dust temperature. 
Fig. 13
Mean dust temperature as a function of A_{v} for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote different dust compositions. 
Parameter of the WD01 and MRN dust models.
3.3. H_{2} formation
Molecular hydrogen is the most abundant molecule in the universe. The formation of H_{2} most effectively takes place on the surface of dust grains. The proper treatment of the H_{2} 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_{2} formation efficiency R_{H2} 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 atomsurface binding resulting from a Van der Waals force (dipoledipole interaction). Chemisorption is a strong atomsurface 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 (T ≤ 100 K), H_{2} formation involves the migration of physisorbed H atoms. At higher temperatures (T ≥ 100 K), H_{2} formation results from chemisorbed H recombination. The presence of these two types of binding sites allows H_{2} formation to proceed relatively efficiently at low and elevated temperatures.
For the dust component i, the formation rate can be written as (19)where n_{H} and v_{H} are the number density and thermal velocity of H atoms in the gas phase, n_{d,i}σ_{d,i} is the total cross section of the grain component i, ϵ_{H2,i} is the formation efficiency (Cazaux & Tielens 2004) on surfaces of the ith component, and S_{H} is the sticking coefficient of the H atoms (20)which depends on the gas and dust temperature, T_{g} and T_{d} (Hollenbach & McKee 1979). The thermal velocity is 1.45 × 10^{5}(T_{g}/100)^{1/2} cm s^{1}. The ith dust cross section per H nucleus is given by (21)Table 4 lists σ_{d,i} for the different dust models. Cazaux & Tielens (2002) expressed the general expression for the ith H_{2} formation efficiency: (22)ϵ_{H2,i} is limited by three terms: the first term A 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_{2}. The third term B governs the hightemperature regime, where evaporation of physisorbed atoms removes H from the surface before it can recombine to H_{2}. It competes with the hopping of atoms into chemisorbed binding sites where the atoms can remain on the surface to form H_{2} at significantly higher temperatures. The numerator ξ_{i}, describing the chemisorption, terminates the H_{2} formation at much higher temperatures, above several hundred Kelvin, when even chemisorbed atoms start to evaporate.
Fig. 14
H_{2} formation efficiency for carbon and silicate surfaces. A and B denote the temperature regimes of the corresponding terms in Eq. (22). The dashdotted line denotes the standard formation efficiency given by Hollenbach & McKee (1979). 
Fig. 15
H_{2} formation rate R_{H2} as a function of dust temperature T_{d} for different formation models and dust compositions. The gas temperature T_{g} is set to 50 K. The dashed lines indicate thelow temperature formation in the original CT formation rates with A ≠ 0. SD89 is the formation rate from Sternberg & Dalgarno (1989). 
Figure 14 shows ϵ_{H2} as a function of T_{d} for silicate (blue line) and graphite (red line) surfaces. The old formation efficiency from Hollenbach & McKee (1979) is shown as a dashdotted line for comparison. When applying this formation in the framework of the PDR model, see Fig. 14, the term A turned out to be problematic because it shuts down the H_{2} release into the gas phase once T < 10 K, accumulating all molecules on the grain surfaces. We question this term because Eq. (22) does not account for the H_{2} 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 A = 0 to ensure efficient H_{2} formation even on cold dust surfaces (Cazaux 2011, priv. comm.)^{6}. More details on how to calculate the H_{2} formation efficiency are given in Appendix C. The total H_{2} formation rate can then be written as (23)Figure 15 compares the new and old H_{2} formation rates. The formation efficiency given by SD89 follows the treatment of Hollenbach & McKee (1979), which adopted a critical dust temperature of T_{cr} ≈ 65 K above which the H_{2} formation efficiency drops below 0.5. H_{2} formation becomes inefficient above ≈ 100 K. The H_{2} formation in the CT dust model remains efficient up to 1000 K. The model by Cazaux & Tielens allows a much more efficient H_{2} formation at higher temperatures. Additionally, the significantly larger dust surface in the MRN and WD017 model lead to a much higher H_{2} formation rate than SD89.
Fig. 16
H_{2} formation rate R_{H2} as a function of dust temperature T_{d} for different dust compositions. The gas temperature T_{g} is set to 50 K. The dashed lines show the formation rate with suppressed H_{2} formation on PAHs. The solid lines show the formation rate allowing for H_{2} formation on PAH surfaces. 
While accounting for chemisorbed H atoms significantly increases the H_{2} 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_{2} formation rate. The resulting H_{2} formation rate for one H atom cm^{3} is shown in Fig. 16. Solid lines show the formation rate as a function of T_{d} for a constant T_{g} = 50 K for the WD017, WD0121, and WD0125 models allowing for H_{2} formation on the surface of PAHs. The dashed lines show the formation rates if the PAH surfaces are excluded from the H_{2} formation calculation. Consequently, the dust model with the largest effective cross section, WD017 shows the highest formation rate in both cases.
Fig. 17
H_{2} formation heating rate for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models have been computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 
The strong increase in the total H_{2} formation rate has an important effect in addition to the easier H_{2} production. Each formed H_{2} 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 rotationalvibrational excitation. If we assume a uniform distribution, 1.5 eV are available to heat the gas. Here we can measure the H_{2} formation rate across the cloud by equivalently looking at the H_{2} formation heating.
By including H_{2} formation according to Eq. (23) we also increase the formation rate significantly at low values of A_{V}. This is shown in Fig. 17. All models in the plot were computed using the WD017 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_{2} formation treatment from CT (with A = 0), suppressing and allowing H_{2} formation on PAH surfaces, are shown by the dotted and dashed line. The much larger surface area of the WD017 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_{2} formation heating.
It is important to note that the H_{2} formation heating is extensive at cloud depths where most of the hydrogen is in the form of H. This leads to the curious situation of a strong H_{2} formation heating in the absence of H_{2} molecules. The gas temperatures at these depths can reach more than thousand Kelvin. Under these conditions two more H_{2} destruction reaction become important: 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 (24)where k_{(C 1)} and k_{(C 2)} are the temperaturedependent 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_{2} 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 H_{2} formation can take place on PAH surfaces the formation rate from Cazaux & Tielens predicts rates of a few times 10^{16} s^{1} 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 R ∽ 2−4 × 10^{17} cm^{3} s^{1} 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 H_{2} formation formalism in the Meudon PDR code and also found a strong increase of the local H_{2} formation rates. If we prevent H_{2} 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 H_{2} line excitation (Habart et al. 2011). H_{2} formation on PAHs has also been proposed by Habart et al. (2003) to explain the H_{2} 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 T is the gas temperature in Kelvin and ψ is in units of K^{1/2} cm^{3}. To convert the FUV field from Habing G to Draine χ units use G = 1.71χ. These parametrizations are fairly accurate when 10^{3} K ≤ T ≤ 10^{4} K and 10^{2} K^{1/2} cm^{3} ≤ ψ ≤ 10^{6} K^{1/2} cm^{3}. However, the parameter range in PDR models allows ψ to reach values as high as 10^{9} K^{1/2} cm^{3} and as low as 10^{2} K^{1/2} cm^{3}. Weingartner (2009) provided updated calculations of Γ_{pc} and Λ_{gr} for the dust distributions WD017, WD0121, and WD0125 up to ψ = 10^{9} K^{1/2} cm^{3}. These rates are fairly well reproduced by the new analytic fits: (25)and (26)Numerical values for C_{i} and D_{i} are tabulated in the appendix in Tables D.1 and D.2.
Fig. 18
Photoelectric heating rate for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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 Λ_{gr} > Γ_{pc}. 
Fig. 19
Gas temperature T_{g} for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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_{2} formation. Right panel: the different lines show the effect of different H_{2} formation treatment. All models in the right panel were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 
Figure 18 shows the PE heating rates for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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 (26) for the different dust models. Dotted lines indicate negative values, i.e., a transition from heating to cooling because Λ_{gr} > Γ_{pc}. For A_{V} > 1 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, WD0125, shows the strongest PEH rate while the model WD0107 drops quickest. However, it is interesting to note that for A_{V} < 0.1 the PEH rate for the WD0107 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
The changes in computing the PE heating and the H_{2} 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_{2} formation treatment (right panel). Among the different dust models presented by Weingartner & Draine (2001b) the WD017 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 0.1 < A_{V} < 2, 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_{2} 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.
Fig. 20
H and H_{2} abundances for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote different treatments of the H_{2} formation. The solid line shows a model using the standard H_{2} formation rate as given by Sternberg & Dalgarno (1989). The dashed and dotted lines show model results for H_{2} formation according to CT where the formation on PAH surfaces is switched on and off, respectively. All three models assume a WD017 dust distribution. The green dashdotted line assumes an MRN dust model and uses the CT H_{2} formation. 
Overall, we can distinguish two regions: For optical depths A_{V} ≳ 0.5, 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_{2} 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 photoionization or photodissociation processes are affected. Secondly, the stronger H_{2} formation efficiency leads to an increased H_{2} formation heating contribution, which in turn produces an increase in gas temperature at low A_{V}. 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.
Fig. 21
CO and C^{+} abundances for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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_{2} formation. Right panel: the different lines show the effect of different H_{2} formation treatment. All models in the right panel have been computed using the WD017 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 HH_{2} transition zone on the applied H_{2} formation rate. We compare three models assuming a WD017 dust model with full radiative transfer and photoelectric heating treatment. The only difference between the model calculations is the treatment of H_{2} formation. We apply (1) the standard H_{2} formation rate as given by Sternberg & Dalgarno (1989) (black, solid line); (2) the H_{2} formation given by CT with H_{2} formation taking place on the PAH surface (blue, dashed line); (3) and the H_{2} formation given by CT with H_{2} formation on PAHs suppressed (red, dotted line). H_{2} 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_{2} formation influences the details of the HH_{2} transition. In that particular model clump, the dust temperature is about 50 K and almost identical in all models, leading to a H_{2} 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 (σ_{d} = 4.14 × 10^{22} cm^{2}) followed by CT without PAHs and CT including PAHs (see Table 4). This agrees with the HH_{2} 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 WD017 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 WD017 model comes from silicates.
In contrast to H_{2} , 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_{2} formation treatment influences the chemistry. All models in the left panel were calculated using the SD89 H_{2} formation. The transition from C^{+} to CO occurs at different A_{V} 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 A_{V}, while models with a weak dust attenuation, e.g. WD0125, exhibit the same transition deeper in the cloud. Additionally, we note an enhanced CO density at low A_{V} in the WD017 model. This is a result of the higher gas temperature due to the stronger PE heating efficiencies of that dust model.
The right panel in Fig. 21 shows the isolated effect of the different H_{2} formation treatment. All models in the right panel were computed using the WD017 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 photodissociation of CO, is the same across all models. However, the CO abundance at lower A_{V} shows very large differences. The models using the CT dust formation scheme produce more CO at the cloud edge compared to the SD89 H_{2} formation. The CT models provide a stronger H_{2} formation heating contribution. The additional heating term produces a significant population of hot CO in the outer layers of the model clump.
Fig. 22
Total line integrated clumpaveraged intensities of ^{12}CO transitions for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 
Fig. 23
CH and CH^{+} abundances (left) and column densities (right) for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 
Fig. 24
OH and OH^{+} abundances (left) and column densities (right) for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 
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 lineintegrated intensities start to show significant differences for the three models. In this particular case, the strength of the highJ CO transitions is proportional to the total H_{2}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 highJ CO emission lines compared to nonclumpy model attempts.
As an additional example of chemical species that predominantly form in the outer regions of molecular clouds and thus are affected by the H_{2} formation treatment we show in Fig. 23 the CH and CH^{+} abundances and column densities and in Fig. 24 the OH and OH^{+} abundances and column densities. All these species are significantly affected by the strong H_{2} 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 A_{V} ≈ 1. 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_{2} formation modeling with the accompanying formation heating effect may help to resolve the discrepancy between models and observations.
4.2.1. Full model SED
In Fig. 25 we show the full continuum and spectral line flux (in units of Jy) emitted by a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The two panels show the results for different dust models. The top panel shows the result for the dust model with the largest effective surface contributing to the gas heating, the WD017 model, equivalent to a R_{V} = 3.1. The strong gas heating is visible in the large population of intense highJ CO emission lines with J > 10 stemming from the outer layers of the cloud (see Table 5)^{7}. In the bottom panel we show the result for a MRN dust model, visible through the lack of PAH emission features in the NIR.
Fig. 25
Full continuum and spectral line flux emitted by a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The two panels show the results for different dust models. Top: WD017 dust model, equivalent to a R_{V} = 3.1. 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 WD017. The PE heating was calculated according to Bakes & Tielens (1994). 
In Table 5 we summarize the line emission of the WD017 and WD0125 (with and without H_{2} formation on PAHs), as well as the MRN models quantitatively. The differences between the WD017 and WD0125 models are limited to the line emission, the continuum emission is only marginally different. The WD0125 models show weaker highJ CO emission because of the smaller grain surface and correspondingly a weaker H_{2} formation heating. The highJ CO emission of the MRN model falls between the WD017 and the WD0125 dust models without H_{2} formation on PAHs. Allowing for H_{2} formation on PAHs will add a strong heating contribution from the formation processes on the large PAH surface, exciting the highJ CO transitions by a large factor.
5. Summary
We revised the treatment of the dust in the KOSMAτ PDR model code to achieve a consistent description of the dust related physics in the code. The computation of dust properties in numerical models of photodissociation regions afflicts the chemical and physical structure of a model cloud via multiple effects. The major areas where dust properties play an important role are

1.
the optical dust properties, i.e., their influence on the radiativetransfer

2.
the dust temperature (neglecting nonequilibrium heating for the very small grains), which influences surface chemistry, most importantly the formation of molecular hydrogen

3.
the heating and cooling capabilities of dust grains via photoelectric heating and gasgrain collisional cooling.
The effect of changes in these areas on the physical and chemical structure of a model can be profound and we described their respective influence in detail.
We notice two opposite effects: Increasing the fraction of small grains steepens the slope of the extinction curve so that the FUV intensity drops faster when going into the cloud, which leads to a temperature decrease at high optical depths. However, it also provides more surface for PE heating and H_{2} formation, thus increasing the heating efficiency. As a consequence, the gas becomes hotter in the outer layers, but colder inside the cloud when the fraction of small grains and PAHs is increased. Only tracers of the outer layers are able to reliably measure the gas heating efficiency. In contrast to common wisdom, this is not even true for ionized carbon [CII], which traces both regimes.
The most influential modification of the code is the treatment of H_{2} formation on the surface of dust grains. Allowing for chemisorption, i.e., including the EleyRideal effect into the formulation of the H_{2} formation efficiency, significantly increases the overall formation rate of molecular hydrogen. If a significant fraction of the dust consists of very small grains and PAHs, their very large contribution to the total grain surface leads to an additional enhancement of the total H_{2} formation. If the molecule formation is enhanced, the formation heating is enhanced as well, providing a profound heating contribution in the outer layers of the model clumps, which significantly affects the local chemistry.
Other PDR model codes are also able to perform wavelengthdependent UV radiative transfer. Compare for example Le Petit et al. (2006); Abel et al. (2005).
A database of photoionization and photodissociation 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).
van Dishoeck et al. (2006) and van Hemert & van Dishoeck (2008) provide their own fitted values of γ_{j,EvD} (up to an maximum A_{V} = 3) for all species for which they also give tables of σ_{j}(λ).
We ignore deviations for CO because photodissociation of CO is treated separately in any PDR code (van Dishoeck & Black 1988).
Acknowledgments
This work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Os 177/1–1. R.S. acknowledges support partly from the Polish NCN grant 2011/01/B/ST9/02229. We acknowledge the use of the UDfA (a.k.a. UMIST) chemical reaction databases^{8}. Some kinetic data we used were downloaded from the online database KIDA (KInetic Database for Astrochemistry^{9}. We thank the anonymous referee for many remarks that significantly improved this paper. We would like to thank J. Weingartner for the calculation of the photoelectric heating rates and collisional cooling rates for the extended parameter range and for the helpful remarks on the calculation of the total heating rate.
References
 Abel, N. P., Ferland, G. J., Shaw, G., & van Hoof, P. A. M. 2005, ApJS, 161, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Abel, N. P., Hoof, P. A. M. v., Shaw, G., Ferland, G. J., & Elwert, T. 2008, ApJ, 686, 1125 [Google Scholar]
 Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles (New York: John Wiley) [Google Scholar]
 Browning, M. K., Tumlinson, J., & Shull, J. M. 2003, ApJ, 582, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222 [NASA ADS] [CrossRef] [Google Scholar]
 Cazaux, S., & Tielens, A. G. G. M. 2010, ApJ, 715, 698 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T. 1978, ApJS, 36, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T. 2003a, ARA&A, 41, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T. 2003b, ApJ, 598, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Falgarone, E., Ossenkopf, V., Gerin, M., et al. 2010, A&A, 518, L118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flannery, B. P., Roberge, W., & Rybicki, G. B. 1980, ApJ, 236, 598 [NASA ADS] [CrossRef] [Google Scholar]
 Garrod, R. T., Weaver, S. L. W., & Herbst, E. 2008, ApJ, 682, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Goicoechea, J. R., & Le Bourlot, J. 2007, A&A, 467, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gry, C., Boulanger, F., Nehmé, C., et al. 2002, A&A, 391, 675 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Habart, E., Boulanger, F., Verstraete, L., et al. 2003, A&A, 397, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Habart, E., Boulanger, F., Verstraete, L., Walmsley, C. M., & Pineau des Forêts, G. 2004, A&A, 414, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Habart, E., Abergel, A., Boulanger, F., et al. 2011, A&A, 527, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hall, P., & Millar, T. J. 2010, A&A, 517, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Hollenbach, D. J., & Tielens, A. G. G. M. 1999, Rev. Mod. Phys., 71, 173 [Google Scholar]
 Hollenbach, D. J., Takahashi, T., & Tielens, A. G. G. M. 1991, ApJ, 377, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Jura, M. 1974, ApJ, 191, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Kirby, K., Roberge, W. G., Saxon, R. P., & Liu, B. 1980, ApJ, 239, 855 [NASA ADS] [CrossRef] [Google Scholar]
 Le Bourlot, J., Le Petit, F., Pinto, C., Roueff, E., & Roy, F. 2012, A&A, 541, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147 [NASA ADS] [Google Scholar]
 Lepp, S., & Shull, J. M. 1983, ApJ, 270, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Li, A., & Draine, B. T. 2001, ApJ, 554, 778 [Google Scholar]
 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Mennella, V., Hornekær, L., Thrower, J., & Accolla, M. 2012, ApJ, 745, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Qin, S.L., Schilke, P., Comito, C., et al. 2010, A&A, 521, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roberge, W. G. 1983, ApJ, 275, 292 [Google Scholar]
 Roberge, W. G., Jones, D., Lepp, S., & Dalgarno, A. 1991, ApJS, 77, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Roberts, J. F., Rawlings, J. M. C., Viti, S., & Williams, D. A. 2007, MNRAS, 382, 733 [NASA ADS] [CrossRef] [Google Scholar]
 Röllig, M., Ossenkopf, V., Jeyakumar, S., Stutzki, J., & Sternberg, A. 2006, A&A, 451, 917 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saxon, R. P., & Liu, B. 1986, J. Chem. Phys., 85, 2099 [NASA ADS] [CrossRef] [Google Scholar]
 Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Störzer, H., Stutzki, J., & Sternberg, A. 1996, A&A, 310, 592 [NASA ADS] [Google Scholar]
 Szczerba, R., Omont, A., Volk, K., Cox, P., & Kwok, S. 1997, A&A, 317, 859 [NASA ADS] [Google Scholar]
 Taquet, V., Ceccarelli, C., & Kahane, C. 2012, A&A, 538, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tielens, A. G. G. M. 2008, ARA&A, 46, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Dishoeck, E. F., & Black, J. H. 1982, ApJ, 258, 533 [NASA ADS] [CrossRef] [Google Scholar]
 van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771 [NASA ADS] [CrossRef] [Google Scholar]
 van Dishoeck, E. F., Jonkheid, B., & van Hemert, M. C. 2006, Faraday Discussions, 133, 231 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 van Hemert, M. C., & van Dishoeck, E. F. 2008, Chem. Phys., 343, 292 [NASA ADS] [CrossRef] [Google Scholar]
 Weingartner, J., & Draine, B. 2001a, ApJS, 134, 263 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Weingartner, J. C., & Draine, B. T. 2001b, ApJ, 548, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Welty, D. E., Hobbs, L. M., & Morton, D. C. 2003, ApJS, 147, 61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolfire, M. G., Tielens, A. G. G. M., Hollenbach, D., & Kaufman, M. J. 2008, ApJ, 680, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Woodall, J., Agúndez, M., MarkwickKemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yorke, H. W. 1980, A&A, 86, 286 [NASA ADS] [Google Scholar]
Appendix A: Isotropic vs. nonisotropic scattering
Fig. A.1
Left panel: mean intensity versus optical depth τ normalized to the incident specific intensity I_{0}. The solid black curve shows the depth dependence of the mean intensity assuming dust parameters ω = 0.4, g = 0.7, approximately describing the WD017 dust sort. The dashed and dotted curves denote mean intensities using the same albedo ω, but assuming isotropic and pure forwardscattering, respectively. Right panel: absolute error  log (J_{i}/I_{0}) − log (J_{k}/I_{0})  of assuming isotropic g = 0 and pure forwardscattering g = 1 with respect to more realistic values ω = 0.4, g = 0.7 (dashed and dotted lines, respectively). 
The MCDRT code assumes isotropic scattering, i.e., a mean scattering angle g = ⟨ cosΘ ⟩ = 0 when calculating the frequencydependent FUV radiative transfer. However, taking into account the detailed material properties and size distribution of interstellar dust results in different values of g. Reasonable values for interstellar dust in the FUV range fall around g ≈ 0.7 and ω ≈ 0.4 (Li & Draine 2001). In the following we show that even in this case of strong forwardscattering the isotropic scattering treatment poses a clear improvement compared to the pure forwardscattering cased used in our previous model.
To quantify the error we used the spherical harmonics method, presented by Flannery et al. (1980) to calculate the penetration of UV radiation of a spherical cloud. Applying their method, we calculated the attenuation of the mean intensity as a function of optical depth for a given asymmetry factor g and albedo ω. For more details on the method see also Roberge (1983), Le Petit et al. (2006), Goicoechea & Le Bourlot (2007).
In Fig. A.1, left panel, we plot the solution of the mean intensity for ω = 0.4, g = 0.7 normalized to the incident radiation field for a homogeneous spherical cloud with an optical depth to the center τ_{c} = 20 (black curve). The dashed blue curve shows the solution for the same albedo, but for isotropic scattering. This represents the solution of the MCDRT code. The dotted red curve gives the result for pure forwardscattering, i.e., the case exp( − (1 − ω)τ). This is equivalent to the previous approximation in KOSMAτ. The right panel shows the corresponding error when assuming isotropic or pure forwardscattering. It is apparent that the assumption of isotropic scattering is an improvement over our previous approach up to optical depths of about 16, only for very high τ the forwardscattering is the solution closer to the g = 0.7 case. However, the FUV influence at τ > 10 is negligible. Presently, we therefore sacrifice the dark cloud accuracy in favor of a more accurate description of the cloud surface, but in a future work we plan to update the MCDRT code to properly account for nonisotropic scattering.
Appendix B: Rescaling of photorates to account for different dust properties
Fig. B.1
Wavelength dependence of the assumed UV illumination and the photodissociation cross sections. The thick black line shows the unshielded emptyspace Draine FUV field, the dashed line shows an extension of the Draine field longward of 2000 Å suggested by van Dishoeck & Black (1982). The blue, green and red lines show the FUV radiation field inside the reference cloud at 95% R_{tot} for the WD017, WD0121, and WD0125 dust models, respectively. Left panel: the blue filled curve shows the continuum absorption cross section σ_{cont} of CH^{2} that leads to dissociation, the black dots with drop lines show the corresponding line absorption cross section σ_{line}. Right panel: same as left panel but for the SiH^{+} photodissociation cross sections. 
In the following we discuss the uncertainties of the proposed rescaling for different dust models D.
One important factor is the spectral shape of the assumed FUV radiation field. In Fig. B.1 we plot the photodissociation cross sections σ of two different molecules (CH^{2} in the left panel, SiH^{+} in the right panel) together with the FUV radiation field (thick black line: unshielded Draine field, colored lines: FUV field inside the clump at 95% of the total radius for the dust models WD017, 21, and 25). The total photodissociation cross section is calculated using Eq. (13). The different wavelength behavior of σ in the two panels demonstrates why the photodissociation rates of different species are attenuated differently when going deeper into the cloud. The right panel also emphasizes that the choice of the spectral shape of the illuminating FUV radiation field may influence the photodissociation rate of some species significantly; the choice between the two FUV fields results in an unshielded photodissociation rate difference of a factor two. However, differences in the attenuation behavior between different dust models will still be dominated by wavelengths shorter than 2000 Å because A(λ)/A(V) varies only weakly for λ > 3000 Å.
Fig. B.2
Comparison of the explicitly integrated photodissociation rates of CH^{2} for the WD017, WD0121, and WD0125 dust models (circles, squares, triangles, respectively) and the corresponding fits to the rates according to Eq. (14). The different line styles denote fits up to different maximum visual extinctions (solid: A_{V,max} = 3, dashed: A_{V,max} = 10, dotted: A_{V,max} = 15, and dotdashed: A_{V,max} = 20). 
Another uncertainty results from the choice of A_{V,max} up to which the γ_{j,fit} from Eq. (14) is fitted to the values Γ_{j}(A_{V})/Γ_{j}(0). In Fig. B.2 we show the Γ_{CH2}(A_{V})/Γ_{CH2}(0) for three dust models: WD017, WD0121, and WD0125. The lines shows the corresponding fits for different A_{V,max}. The range of is 1.92−2.24, 1.48−1.71, and 1.14−1.3 for the three dust models, equivalent to a 10−15% uncertainty. We calculated fits to for all available species and for all dust models to A_{V,max} = 10.
The third source of uncertainties is the functional dependence used in the fit through Eq. (16). We performed leastsquares fits to polynomials up to order 4 to derive a suitable scaling relation between the input parameters γ_{j}, describing the wavelengthdependent attenuation of the photodissociation rate of species j with A_{V}, , describing the wavelengthdependent attenuation of the FUV radiation field depending on the dust model D, and the output parameter . Figure B.3 visualizes Eq. (16) (black grid) and the explicitly calculated fits (colored dots) for the 25 WD01 dust models for all species with available cross sections. The data points are colored according to their deviation from Eq. (16).
Fig. B.3
Visualization of Eq. (16) (black grid) and the explicitly calculated fits (colored dots) for the 25 WD01 dust models for all species with available cross sections. The data points are colored according to their deviation from Eq. (16). Green and blue denotes deviations of less than 0.05 and 0.1. Red points signal deviations exceeding 0.1. All deviations are smaller than 0.22. 
Seven species, CH^{+}, SH^{+}, OH^{+}, HCO^{+}, CO, O2^{+}, and SiO, deviate significantly from our Eq. (16) with residual > 0.4. By shifting their γ_{j} to a corrected we are able to compensate for these deviations. Figure B.4 shows the fitted of the three strongest outliers and the corresponding from Eq. (16) using the original γ_{j} (dashed lines) and the corrected (solid lines).
The deviations of the seven outlying species are large compared to the rest of the 63 species that we calculated (compare Fig. B.3). In the following we will discuss the strongest outliers individually^{10}.
HCO^{+}
HCO^{+} showed the strongest maximum deviation of 1.12 from the plane defined by Eq. (16). UDfA gives γ = 2.0 referring to van Dishoeck et al. (2006) who give a value of γ = 3.32 in their Table 2. Unfortunately, all photodissociation rate coefficients from UDfA that refer to van Dishoeck et al. (2006) are inconsistent with their published numbers and should be replaced with more recent numbers. Their paper was still in preparation when Woodall et al. (2007) was published, so the mismatch is possibly due to prepublication updates of the final values.
OH^{+}
OH^{+} showed the secondstrongest deviations of up to 0.8 from Eq. (16). It possesses a strong absorption continuum only for wavelengths shorter than 1100 Å so that extensions to the Draine FUV field at longer wavelengths should not affect the final photodissociation rate. However, computing the unshielded photodissociation rate coefficient α for a Draine field yields 1.27 × 10^{11} s^{1} , comparable to 1.1 × 10^{11} s^{1} from van Dishoeck et al. (2006), while UDfA lists 1 × 10^{12} s^{1}. The reason for the mismatch is unclear, all sources refer to the same photodissociation cross section by Saxon & Liu (1986). We also find inconsistent values of γ published, between 1.8 and 3.5. We find a best fit to Eq. (16) with .
Fig. B.4
Comparison of the for the three strongest outliers CH^{+}, OH^{+}, and HCO^{+} and the corresponding scaling behavior of Eq. (16) using the original γ_{j} (dashed line) and the corrected . 
CH^{+}
For CH^{+} deviations from Eq. (16) of up to 0.6 were found. Woodall et al. (2007) give α = 2.5 × 10^{10} s^{1} and γ = 2.5 referring to Roberge et al. (1991), van Dishoeck et al. (2006) give α = 3.3 × 10^{10} s^{1} and γ = 2.94 using cross sections from Kirby et al. (1980). We find and . At A_{V} these lower values of γ lead to photodissociation rates 3 − 6 orders of magnitude stronger compared to rates calculated with the higher γ from above. The reason for the γ offsets is unclear.
SH^{+}
SH^{+} has comparable values of α, and γ in UDfA and in van Dishoeck et al. (2006). Small differences in α can be attributed to a strong absorption line at 3100 Å, which leads to a stronger unshielded photodissociation rate when applying a 10 000 K blackbody radiation field as used by Woodall et al. (2007). We find and . Differences of γ from the literature and our calculations most likely result from different assumed dust properties and fits to lower maximum A_{V,max}.
SiO and O
Both molecules show maximum deviations of 0.4 from Eq. (16). The unshielded SiO rate coefficients from UDfA and van Dishoeck et al. (2006) differ by a factor 16 (see comment on HCO^{+}) while their γ_{j} are the same.
In Fig. B.5 we show a histogram of the residuals after removing the outlying species. The standard deviation is σ = 0.05, the maximum absolute residual is 0.22. Photodissociation cross sections are available for a number of species that are not yet included in UDfA (van Hemert & van Dishoeck 2008). In Table B.1 we give α and γ from http://home.strw.leidenuniv.nl/~ewine/photo as well as our .
Fig. B.5
Histogram of the residuals after removing outliers. 
R ate coefficients and corrected for species not included in UDfA. PI and PD denotes photoionization and photodissociation.
Appendix C: H_{2} formation rate
We followed the model of H_{2} formation on interstellar dust grains via physisorption and chemisorption from Cazaux & Tielens (Cazaux & Tielens 2002, 2004, 2010).
For a given dust type, the total H_{2} formation rate is given by Eq. (19) (omitting the subscript index i in the following) where the formation efficiency is given by Eq. (22): (C.1)We set (μF)/(2β_{H2}) to zero to make sure that newly formed H_{2} molecules are able to leave very cold dust surfaces, equivalent to Eq. (13) in Cazaux & Tielens (2002). We can approximate ξ and β_{HP}/α_{PC} by (C.2)(C.3)Equation (C.3) takes into account the corrections from Cazaux & Tielens (2010). Table C.1 gives the applied model parameters for silicate and graphite surfaces.
Model parameter for silicate and graphite surfaces.
Appendix D: Photoelectric heating rates
Weingartner & Draine (2001a) showed that the photoelectric heating rate for the WD01 grain size distributions is fairly well reproduced by the parametrization (D.1)with , where T is the gas temperature in Kelvins and ψ is in units of K^{1/2} cm^{3}, and G = 1.71χ. Numerical values for the parameters C_{i} in Eq. (D.1) are given in Table 2 in Weingartner & Draine (2001a) for each dust distribution. The rate of cooling due to charged particle collisions can be approximated by (D.2)where . Values for the D_{i} for all WD01 grain distributions are given in Table 3 in Weingartner & Draine (2001a).
Fig. D.1
Fit results to the photoelectric heating rate, calculated for an extended range of . The different lines show results for gas temperatures 100 K, 1000 K, and 10 000 K. Γ_{pc} is given in units of 10^{26} erg s^{1}. 
Fig. D.2
Fit result to the net cooling rate due to collisions with charged particles, calculated for an extended range of . The different lines show results for gas temperatures 100 K, 1000 K, and 10 000 K. Γ_{pc} is given in units of 10^{26} erg s^{1}. 
Equations (D.1) and (D.2) are fairly accurate when 10^{3} K ≤ T ≤ 10^{4} K and 10^{2} K^{1/2} cm^{3} ≤ ψ ≤ 10^{6} K^{1/2} cm^{3}. However, the parameter range in PDR models allows ψ to reach values as high as 10^{9} K^{1/2} cm^{3} and as low as 10^{2} K^{1/2} cm^{3}. Weingartner (2009) provided updated calculations of Γ_{pc} and Λ_{gr} for the dust distributions WD017, WD0121, and WD0125 up to ψ = 10^{9} K^{1/2} cm^{3}.
We performed numerical fits to the updated photoelectric heating rates and to the cooling rate due to collisions with charged particles according to the parametrizations 25 and 26. The heating and cooling rates have been calculated by Weingartner (2009) with an extended parameter range of 10^{2} K^{1/2} cm^{3} ≤ ψ ≤ 10^{9} K^{1/2} cm^{3}. The details of the computation are described in Weingartner & Draine (2001a). To achieve a good fit over the full extended parameter range it was necessary to modify the algebraic form of the parametrizations relative to its original form in Weingartner & Draine (2001a). The final form is given in Eqs. (25) and (26). In Figs. D.1 and D.2 we show the fits to calculations of Γ_{pc} and Λ_{gr} for gas temperatures of 100 K, 1000 K and 10 000 K to demonstrate the quality of the parametrization. The C_{i} and D_{i} parameters are given in Tables D.1 and D.2 and will provide a fairly accurate reproduction of the total photoelectric heating rate Γ_{tot} between 10 K ≤ T ≤ 10 000 K and 10^{2} K^{1/2} cm^{3} ≤ ψ ≤ 10^{9} K^{1/2} cm^{3}.
All Tables
R ate coefficients and corrected for species not included in UDfA. PI and PD denotes photoionization and photodissociation.
All Figures
Fig. 1
Extinction cross section per H nucleus for the WD017 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. 

In the text 
Fig. 2
Relative extinction for the dust models of MRN (solid) , WD017 (longdashed; red in electronic version); WD0121 (shortdashed; green in the electronic version) and WD0125 (dotdotdashed; 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 k_{FUV} for the considered models. 

In the text 
Fig. 3
Mean intensity J_{λ}(r) for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The lines show J_{λ}(r) for radii in steps of 10% of R_{tot}. 

In the text 
Fig. 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. 

In the text 
Fig. 5
Mean intensity of the UV radiation J(r) scaled to the unshielded (empty space) Draine field χ_{0} vs. cloud depth. Different lines denote different dust compositions. 

In the text 
Fig. 6
Fit results for for HCO (black), CH_{2} (red), CN (green), and N_{2} (blue) for all dust models from Weingartner & Draine (2001b,their Table 1). The abscissa gives the number of the dust model. The large points show the original values of γ_{j} for the respective molecules. 

In the text 
Fig. 7
Comparison of explicitly fitted for HCO (black), CH^{2} (red), CN (green), and N_{2} (blue) for all dust models from Weingartner & Draine (2001b,their Table 1) with scaling relationship Eq. (16). The large points show the original values of γ_{j} from UDfA, the open diamonds show the values calculated using an MRN dust model, and the open squares denote the explicitly calculated values assuming . 

In the text 
Fig. 8
Dust temperatures of silicate grains for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000 using the WD017 dust model. The different lines denote different grain sizes. 

In the text 
Fig. 9
Dust temperatures of carbonaceous grains for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000 using the WD017 dust model. The different lines denote different grain sizes. 

In the text 
Fig. 10
Dust temperatures of neutral PAHs for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000 using the WD017 dust model. The different lines denote different PAH sizes. 

In the text 
Fig. 11
Grain temperature for the components of the MRN and WD0107 dust models as a function of the radius a_{i} in nm. The grain temperatures are calculated at the surface of a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote different dust components. Dashed lines show the MRN components and solid lines show the WD0107 components. 

In the text 
Fig. 12
Mean dust temperature for the WD017 dust model for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote the average temperature of the individual dust components. The black line shows the overall mean dust temperature. 

In the text 
Fig. 13
Mean dust temperature as a function of A_{v} for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote different dust compositions. 

In the text 
Fig. 14
H_{2} formation efficiency for carbon and silicate surfaces. A and B denote the temperature regimes of the corresponding terms in Eq. (22). The dashdotted line denotes the standard formation efficiency given by Hollenbach & McKee (1979). 

In the text 
Fig. 15
H_{2} formation rate R_{H2} as a function of dust temperature T_{d} for different formation models and dust compositions. The gas temperature T_{g} is set to 50 K. The dashed lines indicate thelow temperature formation in the original CT formation rates with A ≠ 0. SD89 is the formation rate from Sternberg & Dalgarno (1989). 

In the text 
Fig. 16
H_{2} formation rate R_{H2} as a function of dust temperature T_{d} for different dust compositions. The gas temperature T_{g} is set to 50 K. The dashed lines show the formation rate with suppressed H_{2} formation on PAHs. The solid lines show the formation rate allowing for H_{2} formation on PAH surfaces. 

In the text 
Fig. 17
H_{2} formation heating rate for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models have been computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 

In the text 
Fig. 18
Photoelectric heating rate for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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 Λ_{gr} > Γ_{pc}. 

In the text 
Fig. 19
Gas temperature T_{g} for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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_{2} formation. Right panel: the different lines show the effect of different H_{2} formation treatment. All models in the right panel were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 

In the text 
Fig. 20
H and H_{2} abundances for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines denote different treatments of the H_{2} formation. The solid line shows a model using the standard H_{2} formation rate as given by Sternberg & Dalgarno (1989). The dashed and dotted lines show model results for H_{2} formation according to CT where the formation on PAH surfaces is switched on and off, respectively. All three models assume a WD017 dust distribution. The green dashdotted line assumes an MRN dust model and uses the CT H_{2} formation. 

In the text 
Fig. 21
CO and C^{+} abundances for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. 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_{2} formation. Right panel: the different lines show the effect of different H_{2} formation treatment. All models in the right panel have been computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 

In the text 
Fig. 22
Total line integrated clumpaveraged intensities of ^{12}CO transitions for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 

In the text 
Fig. 23
CH and CH^{+} abundances (left) and column densities (right) for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 

In the text 
Fig. 24
OH and OH^{+} abundances (left) and column densities (right) for a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The different lines show the effect of different H_{2} formation treatment. All models were computed using the WD017 dust properties regarding the radiative transfer and the PE heating. 

In the text 
Fig. 25
Full continuum and spectral line flux emitted by a model clump of M = 10 M_{⊙}, n = 10^{5} cm^{3}, and χ = 1000. The two panels show the results for different dust models. Top: WD017 dust model, equivalent to a R_{V} = 3.1. 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 WD017. The PE heating was calculated according to Bakes & Tielens (1994). 

In the text 
Fig. A.1
Left panel: mean intensity versus optical depth τ normalized to the incident specific intensity I_{0}. The solid black curve shows the depth dependence of the mean intensity assuming dust parameters ω = 0.4, g = 0.7, approximately describing the WD017 dust sort. The dashed and dotted curves denote mean intensities using the same albedo ω, but assuming isotropic and pure forwardscattering, respectively. Right panel: absolute error  log (J_{i}/I_{0}) − log (J_{k}/I_{0})  of assuming isotropic g = 0 and pure forwardscattering g = 1 with respect to more realistic values ω = 0.4, g = 0.7 (dashed and dotted lines, respectively). 

In the text 
Fig. B.1
Wavelength dependence of the assumed UV illumination and the photodissociation cross sections. The thick black line shows the unshielded emptyspace Draine FUV field, the dashed line shows an extension of the Draine field longward of 2000 Å suggested by van Dishoeck & Black (1982). The blue, green and red lines show the FUV radiation field inside the reference cloud at 95% R_{tot} for the WD017, WD0121, and WD0125 dust models, respectively. Left panel: the blue filled curve shows the continuum absorption cross section σ_{cont} of CH^{2} that leads to dissociation, the black dots with drop lines show the corresponding line absorption cross section σ_{line}. Right panel: same as left panel but for the SiH^{+} photodissociation cross sections. 

In the text 
Fig. B.2
Comparison of the explicitly integrated photodissociation rates of CH^{2} for the WD017, WD0121, and WD0125 dust models (circles, squares, triangles, respectively) and the corresponding fits to the rates according to Eq. (14). The different line styles denote fits up to different maximum visual extinctions (solid: A_{V,max} = 3, dashed: A_{V,max} = 10, dotted: A_{V,max} = 15, and dotdashed: A_{V,max} = 20). 

In the text 
Fig. B.3
Visualization of Eq. (16) (black grid) and the explicitly calculated fits (colored dots) for the 25 WD01 dust models for all species with available cross sections. The data points are colored according to their deviation from Eq. (16). Green and blue denotes deviations of less than 0.05 and 0.1. Red points signal deviations exceeding 0.1. All deviations are smaller than 0.22. 

In the text 
Fig. B.4
Comparison of the for the three strongest outliers CH^{+}, OH^{+}, and HCO^{+} and the corresponding scaling behavior of Eq. (16) using the original γ_{j} (dashed line) and the corrected . 

In the text 
Fig. B.5
Histogram of the residuals after removing outliers. 

In the text 
Fig. D.1
Fit results to the photoelectric heating rate, calculated for an extended range of . The different lines show results for gas temperatures 100 K, 1000 K, and 10 000 K. Γ_{pc} is given in units of 10^{26} erg s^{1}. 

In the text 
Fig. D.2
Fit result to the net cooling rate due to collisions with charged particles, calculated for an extended range of . The different lines show results for gas temperatures 100 K, 1000 K, and 10 000 K. Γ_{pc} is given in units of 10^{26} erg s^{1}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.