© ESO, 2012

1. Introduction

At an important phase of a young low-mass star’s life, the natal envelope has dispersed, but the star is still surrounded by an optically thick protoplanetary disk, which is thought to be the location of planet formation (see Dullemond & Monnier 2010; Armitage 2010, for recent reviews). Our own solar system contains both gaseous giants and rocky planets. It is thus evident that both the dust and the gas content of the disk need to be studied in order to understand planet formation. The feedback of the forming star onto the disk can help decide the disk’s fate. The star irradiates the disk, heating both the gas and dust and changing the chemical composition of the gas by high-energy ultraviolet (UV) radiation and X-rays. Through heating of the disk’s upper atmosphere, a disk-wind is driven, which may eventually disperse the disk (Alexander et al. 2006; Ercolano et al. 2008; Gorti & Hollenbach 2009). The most abundant species tracing the ionizing radiation and warm temperature of a few 100 K are (ionized) atomic carbon (C, C⁺), atomic oxygen (O), and high-J lines of CO (with J_up ≳ 14). Apart from neutral carbon, which is accessible to ground-based telescopes, these species can now be observed using the PACS instrument (Poglitsch et al. 2010) onboard the Herschel Space Observatory (Pilbratt et al. 2010).

The main forms of carbon (C, C⁺, and CO) can thus be observed for the first time and through high-J CO lines also in warm regions. This allows us to probe the carbon chemistry and the physical structure of planet- and comet-forming zones, which were previously inaccessible to us. The carbon budget is important for planet formation, since the terrestrial planets in the inner solar system are known to have a deficit in carbon of three orders of magnitude relative to the solar abundance (Allègre et al. 2001). Indications of a carbon deficit are also found towards other planetary systems (Jura 2008). In the diffuse ISM, about 40−60% of the carbon is locked into amorphous carbon grains (Savage & Sembach 1996; Sofia et al. 2004, 2011). To be consistent with the Earth’s composition, these grains need to have been destroyed before planet formation. In comets, however, carbon grains have been detected, suggesting that grain destruction mechanisms play no role (Lee et al. 2010).

What region of the disk does the far-infrared (FIR) C, C⁺, CO and O lines trace? How does the partitioning between volatile and refractory carbon vary from clouds to disks? The simultaneous measurement of the species tracing a significant fraction of the disk may provide answers to these questions. In addition, besides the hard to observe H₂, these simple species form the most abundant constituents of the gas. Their lines are thus important tracers of the physical structure of the disk and allow us to finally test the predictions made by many existing physical-chemical models of disks. Neutral carbon [C i], for example, is predicted to be strong by models with a range of gas masses and dust settling (Jonkheid et al. 2007) but is not detected in observations of HD 100546 (Panić et al. 2010) and CQ Tau (Chapillon et al. 2010). In addition to the submillimeter lines discussed in our work, neutral carbon also has forbidden lines in the near-infrared (at 8729, 9826, 9852 Å), which are predicted to have detectable strengths (Ercolano et al. 2009b). To date, only a few detections towards pre-main-sequence stars have been reported (Looper et al. 2010).

The physics and chemistry of gas in protoplanetary disks, which is irradiated by the central star, is similar to the interface of general molecular clouds, where nearby massive stars heat and photoionize the gas. Modeling the infrared (IR) response of gas irradiated by UV and X-rays has been developed in the context of these photo-dominated or dissociation regions (PDRs) and X-ray dominated regions (XDRs) (e.g. Tielens & Hollenbach 1985; Sternberg & Dalgarno 1989; Kaufman et al. 1999; Meijerink & Spaans 2005). These models simulate the infrared line emission by solving the coupled system of chemical evolution and thermal balance. The coupling is the result of the cooling rates for the thermal balance depending on the abundance of the coolants, which is determined by the chemical network that itself also depends on the temperature. These physical-chemical models were applied in the context of protostars by Spaans et al. (1995) and Bruderer et al. (2009a, 2010) to the walls along outflows, irradiated by UV radiation arising from either accretion hot-spots or the protostar’s photosphere. For protoplanetary disks, physical-chemical models have been applied by e.g. Kamp & van Zadelhoff (2001), Glassgold et al. (2004), Kamp & Dullemond (2004), Gorti & Hollenbach (2004), Nomura & Millar (2005), Aikawa & Nomura (2006), Jonkheid et al. (2004), Jonkheid et al. (2006), Jonkheid et al. (2007), Gorti & Hollenbach (2008), Woitke et al. (2009a), Woods & Willacy (2009), and Ercolano et al. (2009a).

The well-studied intermediate-mass pre-main-sequence star HD 100546 is the focus of our present study. The B9.5Vne type star of 2.2 M_⊙, L_bol ~ 27   L_⊙, and T_eff = 10   500 K (van den Ancker et al. 1997) is at its distance of 103    ±    6 pc (Hipparcos) one of the closest Herbig Ae/Be sources. The star and disk have been extensively observed at all wavelengths in lines and the continuum. X-ray observations are reported in Stelzer et al. (2006). Both the Far Ultraviolet Spectroscopic Explorer (FUSE) and the International Ultraviolet Explorer (IUE) observed HD 100546 in the UV, detecting absorption in electronic transitions of H₂ probing the warm gas with a temperature of several 100 K (Martin-Zaïdi et al. 2008). Ardila et al. (2007) analyzed scattered-light images in the optical using ACS/Hubble finding evidence of minimal grain sizes larger than those of typical ISM grains. The scattered light observations also reveal structures resembling spiral arms possibly caused by either the perturbations of a companion (Quillen et al. 2005) or a warped disk structure (Quillen 2006). Features in the near-infrared (NIR) reveal a high crystalline-silicate dust fraction in the inner disk surface layers (Bouwman et al. 2003) and an abundant amount of PAHs whose emission is detected on extended scales (Geers et al. 2007). A wealth of continuum images in the near and mid-IR (e.g. Pantin et al. 2000; Grady et al. 2001; Liu et al. 2003; Grady et al. 2005) have helped to constrain the structure of the inner disk using SED fitting. They reveal a large inner hole of radius  ~10 AU (Malfait et al. 1998; Bouwman et al. 2003), which possibly has a Jupiter-sized planet orbiting inside (Acke & van den Ancker 2006). Interferometric observations (VLT-AMBER and MIDI, Leinert et al. 2004; Petrov et al. 2007) confirm the presence of the inner hole directly and set r_hole ~ 13 AU (Benisty et al. 2010). Some CO ro-vibrational emission at 4.7 μm was reported by van der Plas et al. (2009) and Brittain et al. (2009). The latter paper inferred that the CO emission originates from radii of 13−100 AU, which is consistent with the inner hole seen in the dust continuum. They found a hot rotational temperature of  ~1000 K and indications that both UV fluorescence and collisional excitation are important for the excitation. Hot gas in the disk surface layers is also detected by observing the H₂v = 1−0 S(1) line at 2.12 μm by Carmona et al. (2011).

The bulk of the disk mass is at cold temperatures, as traced by submillimeter lines and continuum. Submillimeter continuum was reported by Henning et al. (1994, 1998), who derived a dust mass of a few times 10^-4 M_⊙. The gas up to  ~100 K is traced by CO rotational lines up to J = 7−6, as reported by Panić et al. (2010). They derived a gas mass of  >10^-3 M_⊙ and concluded from the high ¹²CO J = 6−5/¹²CO J = 3−2 ratio that the disk atmosphere needs to be warmer than in T Tauri stars, probably owing to photoelectric heating of the gas by the strong UV field of the forming B9 star. Strong higher-J lines up to CO J = 30−29 were detected by Sturm et al. 2010 using the PACS instrument onboard Herschel. The detected lines have upper level energies of up to 2500 K.

In this work, we study the chemistry and excitation of simple species (C, C⁺, CO, and O) in the disk around a pre-main-sequence Herbig Ae/Be star. Our approach is twofold. We first want to interpret the PACS observations of Sturm et al. (2010) and the ground-based observations of Panić et al. (2010) taken towards HD 100546 and search for evidence of a warm disk atmosphere from the CO ladder, which probes a wide range of temperatures. We however also wish to study the dependence of the lines on certain parameters such as the carbon abundance in the gas phase, the gas/dust ratio, and others. In Sect. 2, we briefly summarize the modeling approach and describe the observations. In Sect. 3, we first present the results for one particular model that agrees reasonably well with the observations in Sect. 3 and then discuss the parameter dependences in Sects. 4 and 5. In Sect. 6, we summarize our results.

2. Model and observations

For this work, we use a new model based on Bruderer et al. (2009b,a, 2010), in which several changes and improvements have been implemented. The section provides a short summary of the model. Details and results of benchmark tests are reported in Appendix A.

Our physical-chemical model starts with (i) a dust radiative transfer calculation for a given dust and gas distribution to obtain the local UV radiation field and the dust temperature. Using an initial guess of the gas temperature; (ii) the chemical abundances are calculated, which serve as input for a molecular/atomic excitation calculation to obtain the cooling rates. Next; (iii) the thermal balance is calculated to obtain an improved guess of the gas temperature. Steps (ii) and (iii) are repeated until convergence is achieved. Raytracing and convolution to the telescope beam then yields line fluxes that can be observed with observations. These modeling steps are summarized in Fig. A.1.

Differences from the study of Bruderer et al. (2009b,a, 2010) are that we solve the chemical network for the steady-state solution instead of interpolating from a pre-calculated grid of abundances (Bruderer et al. 2009b). This new approach allows us to include detailed photodissociation and ionization rates based on the local wavelength-dependent UV field and the actual column densities for the self-shielding factors. The thermal balance is refined by obtaining atomic and molecular cooling rates from the excitation calculation used in Bruderer et al. (2010). This (approximate) method, which is similar to that of Poelman & Spaans (2005), allows us to calculate the molecular excitation quickly and with reasonable accuracy. We can now include the detailed velocity structure into the excitation calculation and account for heating of the gas by the IR pumping of molecules, which is then followed by collisional deexcitation. The method couples cells at different radii of the disk, unlike the “vertical direction only” escape probability approach employed in most previous disk modeling. The disadvantage of our method is that we need to iterate over the entire chemical structure to derive the local gas temperature, which increases the calculation time. However, the method allows us to calculate models with an arbitrary geometry including e.g. disks embedded in envelopes. Since our code contains all geometrical properties in one module, it can be extended to three-dimensional structures in the future. The chemical network in the new code has also been extended by incorporating hydrogenation reactions on grain surfaces, photodesorption of species frozen-out onto grains, and vibrationally excited molecular hydrogen (), which can be used to overcome activation barriers.

In this study, we concentrate on the rotational ladder of CO and do not consider vibration-rotation lines detected in the M-band at 4.7 μm, since fluorescent excitation by UV photons plays a significant role in their population. We note that the excitation of the rotational ladder is determined by collisions and the relative population of CO in the ν = 1 state is small at 4% for a thermalized population at 1000 K (Brittain et al. 2009). Thus, the results presented in this work would not change by including vibrationally excited levels of CO. For the same reasons, we also do not model the H₂ 2.12 μm line. In addition, this line is very sensitive to the depth at which the dust becomes optically thick and a quantitative comparison to the far-infrared lines is thus difficult.

2.1. Parameters for the HD 100546 model

To model HD 100546, we chose to use the dust density distribution obtained by Mulders et al. (2011). They fit the SED using a dust radiative transfer calculation with a vertical structure that obeys hydrostatic equilibrium assuming that T_gas = T_dust. Dropping this assumption results in a density structure with a puffed-up inner rim with about 10 AU (Kamp et al. 2010), but this region does not contribute significantly to the emission we probe here. For consistency, we also use the dust temperature of Mulders et al. (2011). The surface density power-law is taken to be Σ ∝ r^-1, which fits the VISIR images. Another dust density structure for the HD 100546 disk was proposed by Benisty et al. (2010), based on a analytical disk structure and fitted to the IR interferometric data. Their main conclusions about the inner disk were incorporated into the Mulders et al. (2011) model that we use. It is well-known that SED fitting yields degenerate solutions and the obtained dust structure does not necessarily reflect the gas structure, e.g. owing to dust settling effects. In addition, the structure of the disk is likely non-axisymmetric (Panić et al. 2010; Quanz et al. 2011) and with substructures (Ardila et al. 2007). However, this approximate approach is a reasonable one for studying the main characteristics of the gas emission, given the large uncertainties in the gas line modeling (Sect. 5.5).

We adopted a Hipparcos distance of 103 pc to HD 100546 (van den Ancker et al. 1998). The inclination and position angle were taken to be 42° and 145°, respectively (Pantin et al. 2000). We used a Keplerian velocity field around a 2.5 M_⊙ star, which was found by Panić et al. (2010) to closely reproduce the observed spectra. The outer radius is uncertain but estimated to be between 200 AU and 400 AU (Pantin et al. 2000; Panić et al. 2010).

Inside the large inner hole of the HD 100546 disk, a smaller inner disk extending from 0.25 AU to 4 AU was found (Benisty et al. 2010). We run models including the inner disk, but found that more than 99% of the flux of the lines discussed in this work emerges from the outer disk. We thus considered only the outer disk starting at 13 AU for this work. The dust opacity of the inner disk was however included to calculate the local UV field.

As input stellar radiation field (Fig. 1), we analyzed dereddened FUSE observations for wavelengths 911−1200 Å, which had been previously analyzed by Deleuil et al. (2004) and IUE observations (Valenti et al. 2003) for longer wavelengths (1200−3200 Å). The spectra were dereddened using the extinction of A_V = 0.36 (Ardila et al. 2007) and the extinction law of Draine (2003) for R_V = 3.1. The FUSE/IUE spectra were extended to longer wavelengths using the B9V template of Pickles (1998). The X-ray luminosity had been found to be log (L_X) = 28.9 erg s^-1 (Stelzer et al. 2006). Since the plasma temperature had not been well-constrained by observations, we assumed a rather high temperature of kT_X = 6 keV. This hard spectra allows X-rays to penetrate deep into the atmosphere. Nevertheless, the lines studied in this work change by less than  ~5% by switching on X-rays. Thus, the assumption about the shape (plasma temperature) of the X-ray spectra does not affect our results (Sect. 3.1). In addition to the stellar radiation field, we also account for the interstellar radiation field (Draine 1978). Another source of ionization are cosmic rays, for which a rate of ζ_c.r. = 5 × 10^-17 s^-1 is adopted.

The dust opacities control the penetration of the photons and thus both the photodissociation and ionization of molecules as well as the dust temperature. The dust grains used in the SED fitting by Mulders et al. (2011) are silicates with a 5% (in mass) composition of carbonaceous material. For gas/dust  = 100 and a present-day solar photosphere carbon abundance of 2.7 × 10^-4 relative to hydrogen (Asplund et al. 2009), this amounts to about 15% carbon bound in this population of grains. The grains have a size of 0.1–1.5 μm with a distribution proportional to a^-3.5, employing the Mathis, Rumpl and Nordsieck size distribution (Mathis et al. 1977). The optical constants are for irregularly shaped distributions of hollow spheres (DHS, Min et al. 2005). The DHS dust model has been adopted to be consistent with the SED fitting of Mulders et al. (2011). This dust population has a mass of 10^-4 M_⊙, which agrees with observations of small grains by Dominik et al. (2003). However, much of the mass can be in larger grains, and the gas/dust ratios given relative to this population of grains are thus upper limits. This population of grains does not contain the UV/visual attenuation by very small grains (VSGs) and polycyclic aromatic hydrocarbons (PAHs). We thus also use R_V = 3.1 and 5.5 opacities after Draine (2003) at wavelengths  ≲ 1 μm (Fig. 2) to include VSG/PAH absorption. Using these opacities at short wavelengths only affects (i.e. increases) the dust temperature in the uppermost layer, above the region from which the line emission discussed in this work emerges, as we checked. For simplicity, we thus use the same dust temperature for all models taken from Mulders et al. (2011).

The PAH abundance affects both the heating by means of the photoelectric effect, as well as the ionization balance, by means of recombination reactions. We adopt a 5% PAH-to-dust ratio within 100 AU and 1% outside 100 AU, following a simultaneous fitting of Spitzer and VISIR N-band observations using the N_C = 100 opacities of Draine & Li (2007) (Gijs Mulders, pers. comm.). The PAH-to-dust mass ratio of 5% corresponds approximately to the average Galactic PAH-to-dust mass ratio of 4.6%, which in turn represents about 50 ppm of C per H nucleon bound in PAHs (Draine & Li 2007) and 15% of the present-day solar photosphere carbon abundance (for gas/dust  = 100). The PAHs are only excited by UV in the uppermost layer of the disk (Visser et al. 2007) and the observed features thus trace this region. In addition to the large uncertainties in the molecular properties of the PAHs, the PAH abundance is only poorly determined in the region from which the lines discussed in this study emerge. We thus choose to scale the PAH abundance with the gas/dust ratio, such that gas with gas/dust  = 20 contains five times smaller amount of PAHs than for gas/dust  = 100 (Jonkheid et al. 2007). We examine the influence of the PAH abundance separately in Sect. 4.3. We assume PAHs to be photodestroyed in regions with integrated far ultraviolet (FUV) fluxes  >10⁶ ISRF, following the results of a more complete PAH chemistry by Visser et al. (2007). This affects however only the upper layer of the disk and we find that the lines discussed in this work are unaffected by this assumption. We note that the emission of all lines studied here emerges from z/r lower than the region of PAH destruction (Sects. 3.1, 3.4).

2.2. Observations

Spectrally resolved observations of the low-J¹²CO J = 7−6, J = 6−5, and J = 3−2 lines were obtained with the 12 m APEX telescope by Panić et al. (2010) and displayed a double-peaked line shape, that is consistent with a rotating disk. We note that the uncertainty in the J = 7−6 line is considerable and the J = 6−5 line is detected at a much higher signal-to-noise ratio. They have also searched for the 370 μm (809 GHz) [C i]  line, but were unable to detect it. We note that they have been using position switching with an off position 16′ from the source. Thus, it is unlikely that the emission at the source had been canceled out by extended emission. The 609 μm (492 GHz)  line of neutral carbon was not observed.

The [C ii] 158 μm , the [O i] 63 μm , and 145 μm lines and mid-J/high-J lines of CO with J_up = 14 to 30 were observed by Sturm et al. (2010) with the PACS instrument onboard Herschel as part of the “Dust, ice, and gas in time” (DIGIT) key program (PI: N. Evans). The CO J = 31−30 line was blended with the OH 3/2,7/2⁺−5/2⁻ line and CO J = 23−22 was potentially blended with H₂O 4₁₄−3₀₃. These two lines are thus taken as upper limits. We note that the apparent jump in CO line flux at J > 25 may have been introduced by uncertainties in the PACS flux calibration when the first PACS results were published, including those of Sturm et al. (2010). The larger uncertainty in these lines does however not affect any of our conclusions. The PACS data are spectrally unresolved. Except for [C ii], all lines are consistent with their emergence from within the central  pixel. To correct for extended emission in [C ii], we subtracted the pixels neighboring the central pixel. We also accounted for emission leaking from the center pixel owing to the extended point spread function. The origin of the [C ii] line is discussed in Sect. 5.3.

3. Results: a representative model

We now present the details of a representative model of HD 100546. By representative, we mean that it reproduces sufficiently well the observations enabling us to study the main physical and chemical effects. It is however not meant to be a best fit model. The dependence on the parameters assumed here is discussed in the next section. The outer radius of the representative model is assumed to be 400 AU. We define the depletion factor δ_C as the carbon abundance in the gas phase relative to the present-day solar photosphere abundance of 2.7 × 10^-4 relative to hydrogen (Asplund et al. 2009). In this section, we assume a carbon gas phase abundance of 1.3 × 10^-4 relative to hydrogen (δ_C ~ 0.5), which is characteristic of the diffuse ISM. Adopting gas/dust  = 20, we get a gas mass of  ~2 × 10^-3 M_⊙. We use the R_V = 5.5 dust extinction law for this model.

3.1. Physical structure

The physical conditions of the representative model are shown in Fig. 3. Only one quadrant of the disk is given and only regions with gas density  >10⁵ cm^-3 are shown. We present the gas density in Fig. 3a/b both in linear (r, z) as well as in (log (r),z/r)-space. The latter has the advantage that the inner part and upper disk can be seen clearly. Direct rays from the star to the disk correspond to vertical lines. The density reaches up to a few times 10¹⁰ cm^-3 in this model in the inner part of the disk midplane.

The X-ray ionization rate is given in Fig. 3c. Despite the low X-ray luminosity of this source of  ~8 × 10²⁸ erg s^-1, the ionization rate due to X-rays in the inner disk is still about three orders of magnitude higher than the standard cosmic ray ionization rate of 5 × 10^-17 s^-1. For the heating and cooling budget of the disk, X-rays can however be neglected, since the luminosity in the FUV band (6–13.6 eV) is at  ~4 × 10³⁴ erg s^-1 much higher than the X-ray luminosity. Nevertheless, the considerable ionization rate may affect the chemistry and for example lead to a destruction of water (e.g. Stäuber et al. 2006). The simple species considered in this work are however little affected by X-rays and we find that the CO, C, C⁺, and O fluxes change by less than  ~5% when X-rays are switched off.

The (integrated) FUV flux in units of the interstellar radiation field (Draine 1978; 1 ISRF  ~2.7 × 10^-3 erg s^-1 cm^-2) is shown in Fig. 3d. Very high FUV fluxes of up to 10⁸ ISRF are reached in the surface of the inner edge of the disk. In the upper atmosphere of the outer disk at r = 400 AU and z/r = 0.3 the FUV flux still reaches 10⁵ ISRF at a density of a few times 10⁵ cm^-3. We note that these are roughly the conditions of the Orion bar PDR (Hogerheijde et al. 1995; Jansen et al. 1995). Similar conditions are also predicted for the outflow walls, the “dense PDR” along the outflow of a high-mass star forming region (Bruderer et al. 2009a).

Figure 3e gives the FUV extinction in units of A_V, obtained by defining τ_FUV = −ln(I_att/I_unatt), with the attenuated FUV flux I_att and the FUV flux corrected only for geometrical dilution I_unatt (e.g. Bruderer et al. 2009a, Eq. (2)). Conversion to A_V is done by scaling with the dust opacities at 2070 Å  and 555 Å  and accounting for A_V ~ 1.086τ_V. The UV extinction is not used during the calculation, because the photoionization and dissociation rates are calculated from the wavelength-dependent UV spectra. The FUV extinction however shows that owing to scattering, FUV photons can penetrate much further into the disk, particularly in the outer parts (e.g. van Zadelhoff et al. 2003).

Figure 3f shows the adopted PAH abundance. The drop in the inner, upper disk is due to the photodissociation of PAHs. The drop at 100 AU is observationally constrained (Sect. 2.1).

The dust and gas temperature are presented in Fig. 3g/h. Figure 3i shows the ratio of the gas to dust temperature. The dust temperature does not exceed  ~250 K, but stays above  ~30 K across the entire model. Thus, CO does not freeze out onto the dust grains. The gas temperature exceeds the dust temperature by factors of up to 50 in the inner upper disk and reaches temperatures of up to  ~6000 K. In the outer upper disk, the gas temperature is still T_gas ~ 1000 K and thus about a factor of ten higher than T_dust. Deeper within the disk, below z/r = 0.15, collisions between gas and dust grains result in a coupling of the gas and dust temperature (T_gas = T_dust).

3.2. Chemical structure

Fractional abundances of C⁺, C, CO, O, H₂, and electrons relative to the total hydrogen density (n_H = n(H) + 2n(H₂)) are shown in Fig. 4. As found by previous chemical models of protoplanetary disks the chemical structure resembles a classical PDR in the vertical direction (e.g. Jonkheid et al. 2007). Below a top layer of atomic hydrogen, H₂ forms, as self-shielding becomes sufficient. In the outer disk, the inward column is also large enough to prevent H₂ from photodissociation and hydrogen is in molecular form.

Carbon in the upper disk (z/r ~ 0.2) is mainly in the form of C⁺, followed by neutral carbon and CO deeper within the disk. An exception is the “warm finger” of CO in the outer disk (r > 100 AU, z/r > 0.25). In this region, CO forms through the reaction of C⁺ with H₂ to CH⁺, followed by a reaction with O to form CO⁺ which then reacts with H₂ to form HCO⁺ and subsequently CO (see Jonkheid et al. 2007). The initiating reaction is strongly endothermic with an energy barrier corresponding to 4640 K, thus requires a high temperature to proceed. The OH formation at high temperature, by means of the reaction of O and H₂, leads to additional CO⁺ for this chain of reactions. The OH reacts with C⁺ to form CO⁺. At the inner edge of the disk at 13 AU below z/r ~ 0.1, CO is shielded from photodissociation owing to shadowing by the inner disk. The higher abundance of CO at r = 13 AU and z/r = 0.2 is due to the adopted H₂ formation rate: since the sticking coefficient decreases above  ~3000 K (Cuppen et al. 2010), the lower gas temperature at higher z/r leads to an increase in H₂ and subsequently in CO.

Atomic oxygen has a considerably high abundance throughout the disk, except for the midplane. This is because the elemental abundance of oxygen is higher than that of carbon and not all oxygen can be locked into CO, which is almost unaffected by photodissociation, by means of self-shielding. In addition, oxygen cannot be ionized directly by FUV photons since its ionization energy is slightly higher than 13.6 eV. In the inner disk with 13 AU  <r < 20 AU and z/r < 0.1, gas and dust are at a temperature  > 100 K and the gas is not exposed to strong FUV radiation. In this part of the disk, oxygen is locked into water. Further out, water freezes out and locks most of the oxygen into water ice. In this region, the carbon is in the form of gaseous CH₄.

The electron fraction roughly follows the abundance of C⁺. Charge exchange with PAHs is another important reaction type for the ionization balance. The PAHs are positively charged in the upper atmosphere, neutral at intermediate heights, and negatively charged at greater depths within the disk due to charging by electron attachment. Since the electron abundance is very low in the midplane, the PAHs can again become neutral. Charge exchange C⁺  +  PAH⁰  →  C  +  PAH⁺ and recombination C⁺ + PAH⁻ → C + PAH⁰ are the two reactions that affect the simple carbon chemistry most, leading a higher abundance of neutral carbon.

3.3. Warm atmosphere versus dust temperature

Line fluxes predicted by the representative model are presented in Fig. 5. The figure gives the CO ladder on the left and the [C i], [C ii], and [O i] atomic fine-structure lines on the right. Atomic fine-structure lines of C and CO lines up to J = 7−6 are convolved to the APEX beam and all other lines to the Herschel beam. In addition to the model with the calculated gas temperature, a model that assumes that the gas temperature is equal to the dust temperature is shown. A factor of two around the observed fluxes is shown by a gray shaded region. We consider this region to be in good agreement with observations given all the uncertainties in data and modeling (see Sect. 5.5).

The model with the calculated gas temperature fits the data well, with the exception of the high-J CO lines with J > 25 (see Sects. 4.5 and 5.5). However, the cool atmosphere with T_gas set to T_dust fails to reproduce the high-J CO ladder. The J = 16−15 transition with E_up ~ 750 K (Table 2) is already underproduced by about an order of magnitude, and the higher-J lines are underpredicted even more. For example, J = 25−24 (E_up ~ 1800 K) is underpredicted by more than two orders of magnitude. The atomic fine-structure lines have upper level energies that are much lower than the high-J lines of CO. Both the [C i] and [C ii] lines have E_up < 100 K and the line fluxes calculated with a cool and warm atmosphere are very similar. The [C ii] line is somewhat underproduced and [C i] overproduced. However, [O i] has an upper level energy of 327 K (145 μm) and 227 K (63 μm). Consequently, the [O i] lines for the cool atmosphere differ by an order of magnitude from those obtained for the warm atmosphere. We conclude that the prominent emission in high-J CO and [O i] atomic fine-structure lines is strong evidence of a higher temperature of the gas in the upper atmosphere of the disk. This agrees to the recent detection of CH⁺ (Thi et al. 2011) and OH (Sturm et al. 2010) towards this disk, which are both formed in endothermic reactions with H₂.

3.4. Origin of the line emission

Where does the observed line emission originate? The angular resolution of Herschel, corresponding to about 1000 AU for 63 μm at the distance of HD 100546, does not resolve the disk. The model however, may help us to locate the region of emission. Figure 6 gives the relative contribution to the observed fluxes. We show the integrated contribution function (ICF), which incorporates all the effects involved in the line formation: abundance of the molecule, excitation, and opacity. In addition, we provide both the radius within which 75% of the emission emerges and the radius outside which 75% of the emission emerges from.

The ICF, discussed in Tafalla et al. (2006) and Pavlyuchenkov et al. (2008), shows the relative contribution to the line emission along a ray. We obtain from the formal solution of the radiative transfer equation (Eq. (A.14)) (1)where the source function is S_ν, the dust source function is S_dust, the opacity to the observer is τ_ν and the opacity within one cell is Δτ_ν. The integration is performed over frequency. For simplicity, we take the disk to be face on and allow the rays to propagate parallel to the z axis (as e.g. in Panić et al. 2010). This avoids complications with the adopted Keplerian velocity profile but conserves the main physical conclusions. The contribution to the observed line flux from an annulus located at a radial distance r is given by 2πI(r)rdr, with the flux I(r) at radius r. Thus, the relative contribution to the observed flux is  ∝ I(r)dr. In Fig. 6, we thus show the contribution function times the radius (r × ICF), scaled to its peak value.

The low-J emission of CO J = 3−2 originates at radii¹ between r ~ 70−220 AU at z/r ~ 0.15−0.25. The density above this region is too low to contribute significantly to the emission, since the regions below are shielded by the line opacity. Thus, the far side of the disk does not contribute much to the emission, except for some radiation in line wings. Owing to the increasing density and temperature towards the inner region at a given z/r, even spatially smaller regions at radii  <100 AU contribute to the observable emission. Mid-J lines, e.g. CO J = 16−15 with E_up ~ 750 K, are not excited in the outer disk and only molecules within r ~ 35−80 AU contribute to the emission. These lines reach line opacities of around unity and the far side of the disk does not contribute to the emission in the innermost regions (radii  <30 AU). The highest-J emission, e.g. CO J = 30−29 with E_up ~ 2800 K, detected towards HD 100546 emerges from the inner r ~ 20−50 AU in a thin layer at the surface of the disk. There, CO is formed but the gas temperature is still sufficiently high to excite the lines. The opacity of these lines is optically thin and the dust extinction is not large. Thus, emission from both the near and far sides reaches the observer.

The atomic fine-structure lines of [C i] and [C ii] emerge mainly from the outer disk at radii  ~150−300 AU. Owing to the higher abundance of those species in the upper atmosphere, combined with the low critical density (~10³ cm^-3) and low upper level energy of the [C i] and [C ii] lines, these lines emerge from regions higher in the disk than the low-J lines of CO. The oxygen [O i] lines have considerably higher critical densities of  ~10⁵−10⁶ cm^-3 and also higher upper level energies. They are thus not excited in the outer disk and the line emission comes from regions within  ~50−210 AU. Since this line is optically thick, most of the emission is from the near side of the disk.

3.5. CO line profile

The predicted shapes of the CO J = 3−2, J = 10−9, J = 16−15 and J = 30−29 lines are shown in Fig. 7. For CO J = 3−2, the profile is given for the APEX beam, while the higher-J lines are calculated for the Herschel beam. The HIFI heterodyne spectrometer onboard Herschel (de Graauw et al. 2010) can spectrally resolve lines down to  <0.1 km s^-1 and covers the frequency range of CO J = 5−4 up to CO J = 16−15. For illustration purposes, we however also show CO J = 30−29, which might eventually be accessible to resolved observations with the GREAT instrument onboard SOFIA.

The modeled CO J = 3−2 line agrees well in terms of width with the spectra observed by Panić et al. (2010). A slight asymmetry in the observed spectra with a stronger blue-shifted peak is however not reproduced by the model. Panić et al. (2010) attributed the asymmetry to an asymmetric temperature structure that might be due to a warped inner disk (Quillen 2006). Modeling these asymmetries is beyond the scope of this study and would require tighter observational constraints on the density structure of the disk.

The CO J = 16−15 and J = 30−29 lines are predicted to be considerably broader than those of CO J = 3−2 and J = 10−9. While CO J = 3−2 and J = 10−9 both have a width of  ~5 km s^-1, the lines of CO J = 16−15 and J = 30−29 have widths of  ~9 km s^-1 and  ~12 km s^-1. This can be understood based on the origin of those lines from regions closer to the star. For the assumed Keplerian velocity field around a 2.5 M_⊙ star, the projected velocity along the semi-major axis is  km s^-1. For a given velocity v_proj in the observed spectra, r corresponds to the maximum radius from which the line photons can emerge. For the half width of the CO lines discussed here (2.5, 4, and 6 km s^-1) we get r = 160, 50, and 27 AU, roughly corresponding to the radii of the main emission found in Sect. 3.4. We note that the CO J = 30−29 line is broader than expected from the projected Kepler velocity at the inner rim. This is due to thermal broadening of the line, corresponding to more than 1 km s^-1 for the gas temperatures of a few 1000 K reached in the inner, upper atmosphere.

We conclude that the width of the CO lines is an interesting indicator of the radial origin of the emission and thus a probe/test of the temperature structure. With current facilities, lines up to CO J = 16−15 can be spectrally resolved, allowing us to study regions within radii of a few tens of AU as suggested by our results.

4. Results: dependence on parameters

This section presents the results of a parameter study. We vary different input parameters of the model, such as the size of the disk or the amount of gas in the disk in order to understand the dependence of the line fluxes on those parameters. Understanding these dependences will help us to interpret observational trends and also to quantify the robustness of the results, given that large uncertainties enter the modeling (Sect. 5.5).

4.1. Varying size, gas/dust ratio, and carbon abundance

Key parameters for the gas modeling are the total amount of gas and the size of the disk. In this section, we vary the outer radius of the gas disk to be either 200 or 400 AU and assume that gas/dust  =  20 or 100. A disk outer radius of r_out = 400 AU and gas/dust  = 20 (= 100) yields a gas mass of 8 × 10^-4 M_⊙ (4 × 10^-3 M_⊙). For an outer radius of r_out = 200 AU, the gas mass is 4 × 10^-4 M_⊙ (2 × 10^-3 M_⊙), thus a factor of two lower than for r_out = 400 AU.

The outer radius of the disk is difficult to determine, as discussed in Mulders et al. (2011). For HD 100546, scattered light has been detected out to  ~1000 AU (Ardila et al. 2007), but it remains unclear whether this emission originates from the disk or a diffuse remnant envelope. Thus, the disk might be significantly smaller. Modeling the low-J CO lines, Panić et al. (2010) constrained the size of the gas disk to be  ~400 AU and the gas mass to be 10^-3 M_⊙. Mulders et al. (2011) employed an exponential cutoff in surface density outside 350 AU. In other disks, Hughes et al. (2008) found that CO J = 3−2 can be more clearly explained when such a exponential cutoff is used. Here, we use a sharp cutoff at r_out = 400 AU (Panić et al. 2010), but we also consider a disk with a radius of only 200 AU (Pantin et al. 2000).

The gas/dust ratio is assumed to be either an ISM value of 100 or a lower value, as previously suggested by Chapillon et al. (2010) to explain the non-detection of [C i] 370 μm and 609 μm lines towards the Herbig Ae star CQ Tau. We adopt gas/dust = 20, which is similar to their choice. As mentioned before, these gas/dust ratios only refer to the population of small grains and are thus upper limits.

The carbon budget of the gas depends on the elemental depletion on dust grains. The warm ISM shows signs of little or no depletion of oxygen and we thus assume that all oxygen not locked up in silicates is in the gas phase. We note that oxygen may freeze-out as molecules onto the dust grains, e.g. as water ice, but this happens only in the cold shielded regions of the disk (see Sect. 3.1). Carbon, however, may be depleted by a considerable fraction by locking it up in refractory material. We chose a carbon abundance of δ_C ~ 0.1, corresponding to an abundance of 2.4 × 10^-5 relative to hydrogen². In this way, we can study the degeneracy between having a higher gas/dust ratio of 100 compared to 20 and a lower carbon gas phase abundance of δ_C = 0.1 compared to δ_C ~ 0.5. For gas/dust  = 100, we also run models with δ_C = 0.05.

In Fig. 8, we show the integrated line fluxes for the models with different radii, gas/dust ratios, and carbon abundances in the gas, as in Fig. 5. We first discuss the models with gas/dust  = 20 shown in the upper panel of the figure. The low-J and high-J lines of CO have different behaviors. As expected, the high-J lines do not depend on the outer radius of the disk since the radiation emerges from the inner disk. The optically thin high-J lines however scale with the amount of carbon in the gas. We note that this is already true for mid-J lines that are marginally optically thick. The optically thick low-J lines depend less on the amount of carbon, with a change between the model of the same radius of about a factor of two. We note that these optically thick lines may still depend on the abundance, since e.g. the line wings may remain thin. The low-J lines also depend on the outer radius of the disk. Changing the radius from 200 AU to 400 AU increases the area of the disk by about a factor of four. The low-J lines however only change by about a factor of two. This is because a considerable part of the emission originates from within 220 AU (Sect. 3.4).

As expected, the fine-structure lines of [O i] for gas/dust  = 20 do not change with the amount of carbon in the gas. They however slightly scale with the outer radius by about the same factor as the low-J lines of CO. This can be understood by these lines emerging from within about the same radial region, as discussed in Sect. 3.4.

The [C i] and [C ii] lines depend on both the outer radius and the amount of carbon in the gas phase by similar amounts. This is the result of both lines being optically thin and emerging from approximately the same region. The line fluxes indeed scale with the amount of carbon in the atmosphere, while they change by about a factor of four with radius.

The models with gas/dust  = 100, shown in the lower panel of Fig. 8, exhibit very similar dependences on the outer radius and amount of carbon in the gas phase. How do the lines change with the larger amount of gas in the atmosphere? The optically thin high-J lines of CO approximately scale with the amount of gas. However, since the thermal balance is also affected by the higher density, leading to slightly higher temperatures, the lines increase by more than a factor of five. As expected, the low-J lines scale by a smaller factor, owing to optical depth effects. The [O i] fluxes approximately scale with gas mass. The carbon fine-structure lines scale by a smaller factor, as the higher density moves the C⁺/CO transition up in the atmosphere resulting in a smaller increase in the column density.

The line ratio [C i]/[C ii] is about a factor of two lower for the higher gas/dust ratio. This is due to different branching ratios in the C⁺ destruction. The radiative association of C⁺ with H₂ leads to , while grain surface reactions of C⁺ with hydrogenated PAHs (PAH:H) leads to CH⁺. Following a chain of reactions with H₂ and electron recombinations, CH⁺ and lead to CH and CH₂. For higher densities, CH and CH₂ are turned into CO rather than being photodissociated (see Tielens & Hollenbach 1985). This only slightly shifts the CO transition, thus affects the CO fluxes only slightly, but does affect the [C i] lines.

There is some degeneracy between models having more gas (gas/dust  = 100) but less carbon (δ_C = 0.1) and models with less gas (gas/dust  = 20) and more carbon (δ_C = 0.5), since they have the same amount of carbon in the atmosphere. These models with the same outer radius agree more closely than about a factor of two in the low-J to mid-J lines of CO and [C ii] (see also Fig. 9, models with R_V = 5.5). The high-J lines of CO exhibit larger changes. As discussed before, neutral carbon changes by about a factor of three.

4.2. Varying the dust opacity

The dust opacities at UV wavelengths are another key parameter of the gas, as they control the penetration of the high-energy photons into the disk atmosphere. This affects both the chemistry and thermal balance, for example by means of the photodissociation/ionization of molecules or the photoelectric effect on the dust grains. The opacity at short wavelengths is dominated by VSGs and PAHs.

Observationally, the size of the smallest grains is difficult to determine. For HD 100546, indications of grain processing and larger minimum grain sizes were found by Ardila et al. (2007). Bouwman et al. (2003) suggested that the size distributions vary with radii, employing a population of smaller grains in the inner 40 AU and larger grains in the outer disk. Using such a grain distribution, Benisty et al. (2010) fit the SED using silicate dust grains and a surface layer of 0.05−1 μm grains at 13−50 AU and an outer disk consisting of 1 μm–1 cm size grains. In contrast Mulders et al. (2011) manage to reproduce the SED by assuming a single grain distribution of 0.1−1.5 μm size, but using a vertical structure in hydrostatic equilibrium rather than a fixed scale height. In addition to radially dependent dust properties, dust settling may also introduce a vertical stratification of dust properties (e.g. D’Alessio et al. 2006), having a similar effect on the gas as a larger scale height of the gas structure and thus a higher gas/dust ratio in the upper atmosphere (Jonkheid et al. 2006).

In this section, we study the influence of the dust opacity at UV wavelengths on the line fluxes. A detailed study of settling and radially dependent dust properties (e.g. Birnstiel et al. 2010) is beyond the scope of this study and we only globally change the dust properties at UV/Visual wavelengths. We assume the dust opacity laws given in Sect. 2.1 (Fig. 2). These are either ISM-like grains with R_V = 3.1, corresponding to a minimum grain size of  ~50 Å, grains with R_V = 5.5 corresponding to some coagulation and grain growth and considerably larger grains with a minimum size of 0.1 μm. The models with a minimum size of 0.1 μm do not contain any VSGs/PAHs and the opacity at UV/visual wavelengths is thus considerably smaller. Photodestruction of VSGs to PAHs in the upper atmosphere (Berné et al. 2009) leads to some additional vertical change in the UV opacity, even though their mass extinction coefficients in the UV are similar (e.g. Pontoppidan et al. 2007, Fig. 2). We thus consider this section as an approach to studying the influence of different UV penetration lengths on the chemistry, independent of whether the absorption is due to VSGs or PAHs or the absence of both.

In Fig. 9, we show the line fluxes obtained using the R_V = 3.1 and R_V = 5.5 dust opacities. The outer radius of the models is 400 AU, and we adopt either gas/dust  = 100 and δ_C = 0.1 or gas/dust  = 20 and δ_C = 0.5. All fine-structure lines change by less than 50% by using the two different opacity laws. The low-J lines of CO change similarly little, but the high-J lines deviate by up to about a factor of two. The changes are very similar to those for all other combinations of outer radii, gas/dust ratio, and carbon abundance.

The results obtained using the much larger 0.1–1.5 μm grains are given in Fig. 10. For gas/dust  = 100, the line fluxes are much larger than the models with R_V = 5.5 dust opacities, with increases of up to an order of magnitude. Thus, the observed line fluxes can no longer be reproduced, not even with the lower carbon abundance. To roughly fit the observations, we thus show models with gas/dust  = 10 in a second panel. The dependences of the models with large grains on the outer radius and carbon abundance is similar those of to the R_V = 5.5 models.

We conclude that the dust opacity, particularly that of the small grains governing the penetration of the UV radiation into the disk atmosphere, has major implications for the line fluxes. What is the reason for this strong dependence? In Fig. 11, we show vertical cuts through the disk in density, temperature, FUV field, and fractional abundance of CO and C⁺ for the models with different dust opacity laws. Comparing the results obtained with R_V = 3.1 and R_V = 5.5 shows that while the gas temperatures in the upper atmosphere are similar, the transition to the cold atmosphere, where gas and dust temperature are coupled, occurs at greater depth. Similarly, the C⁺/CO transition occurs at greater depths within the disk. Both results are caused by the FUV field that can penetrate deeper into the disk, closer to the midplane. While the temperature and abundance structures are qualitatively similar for the R_V = 3.1 and R_V = 5.5 models, they are considerably different for the larger grains: the gas temperature is still coupled to the dust temperature in the inner midplane (50 AU). Further out (at 250 AU), the gas temperature however no longer couples to the dust temperature. This is consistent with the much larger UV field in the midplane, which also explains the shift of the C⁺/CO transition deeper into the disk. Thus, larger grains allow UV radiation to penetrate more deeply into the disk to regions of higher density and also to heat the gas to higher temperatures. The combination of both leads to the considerable increase in line flux.

4.3. Varying PAH abundance

Given the uncertain abundance of PAHs and the difficulties in simultaneously reproducing the [C ii] fluxes and [C i] upper limits, we explore the effect of the PAH abundance on the line fluxes in this section. The PAH abundance affects the models in two ways. First, the ionization structure is changed by PAHs by means of recombination reactions such as C⁺ + PAH⁰  →  C + PAH⁺, with PAH⁰ being neutral and PAH⁺ positively charged PAHs. Second, the gas heating depends on the PAHs by means of the photoelectric effect on small PAHs.

The emission of PAH can be directly observed by its bands in the infrared (e.g. Geers et al. 2007). This emission however only traces regions where PAHs are excited by UV photons. An abundance and excitation study by Visser et al. (2007) shows that in disks around Herbig Ae stars, most of the emission from intermediate-size PAHs emerge from the uppermost layer, above the τ_V = 1 surface. Thus, PAHs do not trace down to regions where the C⁺/C transition is affected by the PAH charge exchange. The PAH abundance is thus rather uncertain.

Figure 12 presents line fluxes obtained with models that contain different PAH abundances in the outer (r > 100 AU) and inner (r < 100 AU) part of the disk. We assume a disk outer radius of 400 AU, an ISM carbon abundance of δ_C = 0.5, and gas/dust  = 20, corresponding to the representative model presented in Sect. 3.

Changing the PAH abundance at radii  >100 AU from 20% ISM to no PAHs, hardly changes the CO ladder. In addition, the [O i] fine-structure lines are hardly affected. The [C ii] line however slightly increases and [C i] decreases by a factor of about two, owing to the missing charge exchange converting C⁺ to C. Changing the PAH abundance in the outer disk to the ISM abundance increases the mid-J lines of CO, owing to the greater heating caused by the higher abundance of PAHs. We note that the high-J lines of CO are less affected than the lower-J lines, as they mostly emerge at radii  <100 AU. The fine-structure lines are only slightly affected with [C ii] being slightly weaker and [C i] slightly stronger, owing to the greater charge exchange. We also tested the influence of the charge exchange with sulfur (C⁺ + S → C + S⁺), which had been traditionally found to influence [C i] in PDRs (Hollenbach & Tielens 1997, Sect. 4.4). For the dense PDR of the disk atmosphere including a density gradient, the influence of sulfur charge exchange is however found to be smaller than that for PAHs.

Varying the PAH abundance in the inner 100 AU does not affect the atomic fine-structure lines. The high-J lines of CO however change by about a factor of two, if the PAH abundance is higher by a factor of five or set to zero. This is a result of either the higher or lower gas temperature caused by the contribution of photoelectric heating on PAHs to the thermal balance.

4.4. The stellar input spectra

The stellar UV radiation field is an important input parameter for the model, as both the heating and the chemistry depend on the input radiation field. We use observed and dereddened IUE/FUSE spectra. However, the dereddening is somewhat uncertain, particularly for shorter wavelengths. We thus also run models using a blackbody of 8000 K, 10 500 K (T_eff of HD 100546), and 12 000 K (Fig. 1). The blackbody was scaled to preserve the bolometric luminosity. Thus, the luminosity in the FUV band changed. While L_FUV/L_bol ~ 0.27 for the observed FUV spectra, it is L_FUV/L_bol ~ 0.024 (0.095, 0.15) for the blackbody of 8000 (10 500, 12 000) K.

In Fig. 13, the line fluxes of the models with the different input spectra are shown. Models with a lower FUV luminosity L_FUV have less heating and thus the mid/high-J CO and [O i] lines are weaker (Sect. 3.3). The [C i]/[C ii] ratio however changes only slightly with the UV color. We note that both C and CO are photoionized/photodissociated at similar wavelengths (λ ≲ 1100 Å) and their rates thus depend similarly on the UV color. For much colder UV spectra with very few CO dissociating photons, the gas temperature is lower owing to more efficient cooling by CO than for the atomic fine-structure lines. We conclude that for the lines discussed in this work, the UV spectrum affects the lines mostly by means of the photoelectric heating on PAHs/VSGs, which is less efficient for lower L_FUV. The [C i]/[C ii] ratio remains however similar and the results related to their flux are unaffected by an uncertainties in the dereddening of the UV spectra.

4.5. Dependence on the H₂ formation rate

A major uncertainty in the chemistry is the H₂ formation rate at high temperature. This rate may affect the line fluxes both in terms of the change in the H/H₂ transition, which also affects the C⁺/C/CO transition and the thermal balance. The heating and cooling processes affected by H₂ are cooling by H₂ lines and heating by means of UV pumping followed by collisional deexcitation, formation, and photodissociation heating. Uncertainties enter the H₂ formation rate in terms of the minimum grain size that determines the grain surface on which H₂ may form, but also through the sticking coefficient and the formation efficiency. The sticking coefficient gives the probability that H atoms stick onto the surface rather than bounce back or evaporate before reacting to H₂. This sticking coefficient and the formation efficiency depend on various parameters (e.g. Hollenbach & Tielens 1999). It may either decrease with the gas temperature when H₂ formation on physisorbed sites dominates, or stay constant when at high gas temperature the neighboring chemisorbed H atoms react and leave the surface relatively bare.

In this section, we calculate the representative model with a sticking coefficient times formation efficiency S × η = 1, that is independent of temperature. We then compare the result to our “standard” H₂ formation rate using the S × η reported by Cuppen et al. (2010), in order to capture the theoretical work of Cazaux & Tielens (2004), Cuppen & Herbst (2005) as well as the observational evidence from Habart et al. (2004). Figure 14 shows the line fluxes calculated with the models assuming different formation rates. The inset in that figure gives the temperature dependence of the H₂ formation rate.

Neither the low-J to mid-J CO lines nor the atomic fine-structure lines depend on the H₂ formation rate. The high-J lines of CO, however, increase by a factor of  ~3 when the sticking coefficient is assumed to be independent of gas temperature. This is consistent with in the inner disk (r ≲ 70 AU, z/r ≳ 0.15) hydrogen being fully molecular (see Figs. 3 and 4), which increases both the gas temperature and CO abundance in this region. We conclude that the predictions of the CO ladder are robust for the low-J and mid-J lines, although some differences due to the uncertain H₂ formation rate are expected for the high-J lines with J_up ≳ 25.

5. Discussion

5.1. Comparison with observations

How do the results of the models compare to the observed fluxes? Studying the outer radius, gas/dust ratio, and carbon abundance in Sect. 4.1 (Fig. 8), we find that the mid-J/high-J CO lines can be reproduced by either a higher gas/dust ratio and a lower δ_C or a lower gas/dust ratio and a higher δ_C. The low-J CO and [O i] lines emerge from the outer part of the disk, hence scale with the outer radius of the disk. It is thus possible to reproduce these lines for both higher gas/dust and lower gas/dust ratios, if the outer radius is adjusted accordingly. For a gas/dust ratio of 20 and δ_C = 0.5, a 400 AU disk better reproduces the [O i] fine-structure lines than the 200 AU disk, which underproduces the flux. For a gas/dust ratio of 100 and δ_C = 0.05/0.1, a smaller disk better reproduces the [O i] observations, but a 400 AU disk does not considerably overproduce the emission.

For both gas/dust ratios and outer radii, the modeled low-J lines agree reasonably well with the observations. Independent of the outer radius, the [C ii] line is underpredicted for the models that fit the CO ladder and [O i]. For gas/dust  = 20, δ_C = 0.5, and r_out = 200 AU, the [C ii] is only underproduced by about a factor of three, but that model overpredicts the [C i] line by about a factor of four. However, for gas/dust  = 100 and δ_C = 0.1, both models with r_out = 200 AU and 400 AU underproduce [C ii], but are within the upper limit to [C i]. The same is true for the model with δ_C = 0.05, which however more closely reproduces the CO ladder. Thus, none of the models are able to simultaneously reproduce both the [C ii] line and the upper limit to [C i]. We conclude that only models with a higher gas/dust ratio and low δ_C ≲ 0.1, are able to simultaneously reproduce the CO ladder, [O i], and the upper limit to [C i] but underpredict [C ii] by a factor depending on the outer radius of the disk.

Changing the dust opacities from R_V = 3.1 to R_V = 5.5 does not change these conclusions (Sect. 4.2, Fig. 9). For models with 0.1–1.5 μm grains and a much lower UV/visual extinction, the models with gas/dust  = 100 overproduce both the CO ladder and [O i]. For gas/dust  = 10 and δ_C = 0.1, both the CO ladder and the [O i] lines are fit reasonably. The [C ii] line is underproduced, but [C i] is within the upper limit for both radii. This is similar to our results for the models with a R_V = 5.5 dust opacity law, gas/dust  = 100, and δ_C = 0.05/0.1.

Decreasing the PAH abundance in the outer disk (radii  >100 AU) mostly affects the [C i]/[C ii] ratio and a model without PAHs in this region starts to approximate both lines (Sect. 4.3, Fig. 12). A higher PAH abundance in the inner disk (<100 AU) increases the flux of the mid-J and high-J CO lines and yields a slightly closer agreement with the high-J lines, despite slightly overproducing some mid-J lines. Given the large uncertainties in the gas temperature (Sect. 5.5), we however conclude that these changes are not significant.

Using an H₂ formation rate with a sticking coefficient that is independent of temperature slightly increases the high-J lines of CO, which are in closer agreement than using our “standard” sticking coefficient that decreases above a certain temperature. However, as discussed before, the uncertain gas temperature most affects the high-J lines and we thus do not consider the change in the high-J lines with H₂ formation rate as significant.

We conclude that the CO ladder together with [O i] and the upper limit to [C i] can be reproduced by models with a high gas/dust ratio and low δ_C. These models however underproduce [C ii]. Models without PAHs in the outer disk, improve the fit to [C i] and [C ii] but there is still a factor of two discrepancy for both lines.

5.2. The [C I]/[C II] ratio and the gas carbon abundance

Given the difficulties in simultaneously reproducing the [C ii] line and the [C i] upper limit, it is interesting to study the dependence of the two lines on the parameters independent of the other lines. The modeled fluxes of the [C i] 370 μm line and the ratio to the [C i] 158 μm line are shown in Fig. 15, which summarizes Figs. 8, 10 and 12.

The dependence on r_out, gas/dust, and δ_C is shown in the upper panel of Fig. 15 for models with the different dust opacities used in this work. The [C i] 370 μm line and the [C i] 370 μm / [C ii] 158 μm ratio depend similarly on r_out, gas/dust and δ_C for the different opacity laws. Exceptions are the large grain models owing to the UV photons penetrating further into the disk, thus dissociating precursors of CO, such as CH and CH₂, and affecting C/C⁺. The [C i] 370 μm line varies by about a factor of four between the models with outer radii of 200 AU and 400 AU and a similar factor between δ_C = 0.1 and 0.5. The [C i] 370 μm / [C ii] 158 μm line ratios are almost independent of both the outer radius and carbon abundance, but depend slightly on the gas/dust ratio, again since CH and CH₂ are more likely to be converted into CO than photodissociated at higher density, thus decreasing the [C i] fluxes. Thus, by increasing the gas mass by a factor of five from gas/dust of 20 to 100 only increases the [C i] by about a factor of two. This reduces the [C i] 370 μm /[C ii] 158 μm line ratio by about a factor of two. For the larger grains, this effect becomes much stronger because the lower UV extinction moves the C⁺/C transition deeper into the disk, to higher densities. We note however, that the models with large grains and gas/dust  = 100 overproduce CO and [O i] (Sect. 4.2).

The PAH abundance within 100 AU does not affect the [C i] and [C ii] fluxes. However, if PAHs are absent in the outer part of the disk (>100 AU), the [C i] fluxes decrease by about a factor of two, and the [C i] 370 μm/[C ii] 158 μm line ratio decreases by about a factor of three. As discussed before, this is due to the lower charge exchange from C⁺ to C by PAHs. Similarly, increasing the the PAH abundance from 20% ISM to ISM abundances leads to a stronger charge exchange and thus a larger line ratio.

In spite of all these parameter changes, the observed [C ii] 158 μm line and the observed line ratio [C i] 370 μm / [C ii] 158 μm cannot be reproduced by our models. Our closest match is found for the model without PAHs in the outer part of the disk (Fig. 12), but that model still differs by a factor of about two for both lines. The [C i] 370 μm / [C ii] 158 μm line ratio is surprisingly independent of the parameters, which can be understood by a similar region of origin of the two lines found in Sect. 3.4.

5.3. The origin of the [C II] emission?

Does the observed [C ii] 158 μm emission really originate from the disk? Since C⁺ is the main state of carbon in diffuse clouds and owing to the low critical density of the 158 μm line, how much could a tenuous remnant envelope or a diffuse foreground cloud contribute to the emission? We note that the discussion here assumes that the emission is centered on the source, within a  ~10′′ (1000 AU) diameter region, since extended emission has already been subtracted (Sect. 2.2).

If the extinction towards the star results from material only seen in the central pixel, how much could that contribute to the [C ii] emission? The maximum amount of emission to be expected from the foreground can be calculated from the foreground extinction of A_V = 0.36 mag (Ardila et al. 2007). Following Appendix B.1, we obtained a flux of  ~2.5 × 10^-16 W m^-2 in the Herschel beam, using a conversion of 1.9 × 10²¹ cm^-2 (Bohlin et al. 1978), a carbon fractional abundance of about 1.3 × 10^-4, and a gas temperature in the foreground above the upper level energy of the [C ii] 158 μm line (91 K). This would already exceed the flux reported in Table 1 by 60%.

How much emission could come from a compact diffuse remnant envelope? Scattered light is detected out to a distance of 1000 AU (Ardila et al. 2007) and could belong to either the disk or a remnant envelope. If a tenuous remnant envelope is exposed sufficiently to UV radiation, most of its carbon is in the form of C⁺. This nebula cannot contain CO, since this would fill in the central emission of the observed double-peak line shape (Panić et al. 2010). While the density structure of such a remnant envelope is not known from observations, assuming a ring size in the range  ~400−1000 AU at a density of  ~10⁴ cm^-3 and a thickness of  ~800 AU would yield a flux of up to 1.4 × 10^-17 W m^-2. Thus, both a foreground cloud and a remnant envelope might well contribute considerably to the [C ii] emission. We note that no [C i] emission is expected from this component, since for typical columns of C towards diffuse clouds (≲ 3 × 10¹⁵ cm^-2, Jenkins & Shaya 1979), the [C i] 370 μm line is below the detection limit of APEX.

This rough estimate shows that both foreground and remnant envelope could add substantially to the [C ii] emission, making the association of the observed flux with the disk not obvious. This question can however only be addressed with velocity resolved spectra of [C ii] 158 μm, obtained with the HIFI instrument onboard Herschel.

5.4. The carbon budget of the disk atmosphere

Our results indicate that for normal gas/dust ratios, the carbon abundance is low at δ_C ≲ 0.1. How does the gas-phase carbon abundance change as the gas evolves from diffuse clouds to dense clouds, collapsing envelopes, and eventually disks? At the earliest stage, in diffuse clouds, the carbon gas-phase abundance can be measured directly by the UV absorption of C ii, without referring to a dust model (e.g. Sofia et al. 2004, 2011). The measured amount of volatile carbon, which is not bound in dust, is constrained to be between  ~40% and 60% (δ_C = 0.4 to 0.6), with some differences between the measurements of the strong 1334 Å  feature (Sofia et al. 2011) and the weak intersystem feature at 2325 Å (Sofia et al. 2004). Assuming that the refractory material does not change as the gas evolves to cores and finally disks, these micron-sized carbon grains quickly settle down to the midplane of the disk. Thus, if a process removes gas from the upper disk, that population of carbon dust is not removed from the disk.

The volatile carbon however undergoes significant changes in the course of this evolution: from diffuse to dense clouds, the main reservoir becomes gaseous CO rather than C⁺. In pre-stellar cores, it becomes CO ice rather than CO gas (Bergin & Tafalla 2007). From pre-stellar cores to disks, it either stays in the form of CO ice, adsorbs or desorbs from grains multiple times, or is transformed into other species such as CO₂ or CH₃OH (Visser et al. 2009b, 2011). If gas is removed from the disk by some process, leading to a smaller gas/dust ratio, the volatile carbon is likely to be removed at a rate that is proportional to the quantity of gas. Thus, within the gas, the amount of volatile carbon should remain the same. Inferring a value of δ_C that is much lower than 0.4 to 0.6 as we find here means that a fraction of volatile carbon has been turned into non-volatile material in the course of the evolution from diffuse clouds to disks. Towards hot cores, where CO ice is supposed to evaporate from dust grains, Alonso-Albi et al. (2010) and Yıldız et al. (2010) found evidence of low CO hot core abundances, suggesting some transformation of CO on the surfaces of the grains into other volatile or non-volatile carbon species. For the warm disk discussed here, CO ice cannot be the main carbon reservoir. In addition, no strong CH₃OH emission is in general observed in disks, ruling out this possible carbon reservoir. However, some “mild” photochemistry could lead to more complex, less volatile organics (Öberg et al. 2009a). Some of this less volatile organic material could become similar to the so-called “CHON particles” detected in comets (Kissel et al. 1986) after further processing in the disk.

The upper layers of the disk atmosphere, probed in this work, more likely have gas/dust  >100 rather than  <100, since settling increases the gas/dust ratio in that region. In addition, as we keep the dust density structure in our work fixed, a small gas/dust ratio would mean that the relative fraction of carbon locked into dust (VSGs/PAHs) increases for small gas/dust ratios, which actually implies that the total carbon abundance is higher than the solar abundance. These arguments drive us to the conclusion that the warm gas in the HD 100546 atmosphere more likely consists of a large gas/dust ratio with a considerably smaller abundance of volatile carbon than in the diffuse ISM.

5.5. Uncertainties in T_gas

An important caveat of the models presented here is the uncertainty in the derived gas temperature. On the surface of the disk where the gas and dust temperature decouple, the equilibrium of different heating and cooling processes govern the gas temperature. As noted before, the thermal balance is a coupled problem together with the chemistry as line cooling (in [C ii], [O i], CO, ...) depends on the abundances. From PDR modeling, the thermal balance problem is known to be very sensitive to the assumptions made, such as the method used for the line radiative transfer or the H₂ physics and chemistry implemented. Röllig et al. (2007) compared the gas temperature derived by different PDR models assuming the same chemical network. They found a relatively large scatter, in particular for the model with a high density (n = 10^5.5 cm^-3) and a strong UV field (χ = 10⁵) corresponding approximately to the upper atmosphere of the outer disk. The scatter between the models is about a factor of three in T_gas (Fig. A.7) and similar for the calculated fine-structure line fluxes. In our present study, we found the high-J lines of CO to be more sensitive to the temperature (Sect. 3.3) and we thus calculated the CO ladder from the CO abundance and gas temperature obtained from the different PDR codes. The scatter in those line fluxes was found to be about 1.5 dex when outliers were removed. The chemical network, which was been taken to be the same as the one used in the benchmark study of Röllig et al. (2007), also introduces uncertainties. Vasyunin et al. (2008) performed a Monte Carlo sensitivity analysis and determined the uncertainty in the CO and C⁺ column density through a disk to be a factor of two or less.

For our work, we thus consider an agreement of a factor of two between observations and model as a good agreement. Given that many model parameters remain free even for this well-studied disk, the observed fluxes can be reproduced to a level that is much better than a factor of two. We however refrain from trying to achieve so, as we feel that given the large uncertainties, these results would not be trustworthy and provide a misleading picture.

Given the considerable uncertainty and disagreement between the models, which of our conclusions is robust? From the considerably higher gas temperature needed to reproduce the high-J CO lines, together with all PDR models agreeing that the gas temperature indeed decouples from the dust temperature at these high densities and UV fluxes, we conclude that there must be a warm gas atmosphere. We also consider the conclusions made in Sect. 4 as robust, since they mainly represent relative changes and trends that are more reliable than absolute fluxes, as long as the correct main heating and cooling processes dominate the thermal balance. Similarly, finding trends in the dependences of the line fluxes on model parameters, which can then be observationally tested using a larger sample of disks, is more promising than fitting the absolute fluxes.

5.6. The shape of the CO ladder

The high-J CO lines emerge from the upper, inner atmosphere of the disk, from regions where heating by FUV radiation leads to a gas temperature that is much higher than the dust temperature (Sect. 3.4). However, what does the shape of the CO ladder and in particular the mid-J and high-J CO lines tell us? In Fig. 16, we show the rotation diagram of the CO ladder from both observations and the representative model (Sect. 3). Goldsmith & Langer (1999) provide a detailed introduction to the rotation diagram method. The Y axis of the diagram is defined to include the solid angle of the emission, because the CO ladder is thought to emerge from different regions (Appendix B.2).

The observed CO ladder can be represented by three temperature components: A cold component with a rotation temperature of T_rot ~ 60 K (J_up = 3−7), a warm component of T_rot ~ 300 K (J_up = 14−25), and a hot component with T_rot ~ 750 K (J_up = 25−30). We did not include the upper limits to the J = 23−22 and J = 31−30 lines in the calculation of T_rot.

The low-J lines are optically thick. This leads to an increase in the flux compared to the optically thin conditions (Appendix B.2). Owing to the opacity effects, the temperature obtained from the cold component is thus not necessarily a kinetic temperature even though the lines are in thermal equilibrium. The model predicts that lines up to E_u ~ 750, J = 16−15 are optically thick (Sect. 6). The modeled CO ladder indeed shows a curved shape, up to about this upper level energy, which is characteristic of optical depth effects. However, a larger emitting area, leading to a larger dΩ can also increase the line flux in a similar way. The emitting region dΩ contributing more than 50% to the total emission (obtained from the 75% radii given in Table 2) increases from J = 16−15 to J = 3−2 by a factor of about eight, corresponding to a shift of ΔY ~ 2 in the rotation diagram.

The optically thin mid-J and high-J lines obtained from the model, can be represented well by a single temperature component of 393 K, which is somewhat higher than the observed value. The difference is however small in comparison to the large uncertainty in the gas temperature. The hot component of T_rot ~ 750 K is not reproduced well by the representative model. A higher H₂ formation rate, leading to a higher H₂ abundance in the hot, inner gas, thus greater heating and more CO, yields an additional hot component. This model more close reproduces the observations of the hot component (Sect. 4.5).

How is the CO ladder affected by the molecular excitation? To test this, we compared the CO ladder from the representative model with the same model, but assumed that the CO level population is in local thermal equilibrium (LTE). While the optically thick low-J to mid-J lines are very similar, the optically thin lines with higher J are stronger by a factor of about nine (ΔY ~ 2) in flux. The critical density of these lines is on the order of up to a few times 10⁶ cm^-3 (Table 2). In the region of the main emission (Sect. 3.4), the density is similar or higher than the critical density. The calculated level population in this region is similar to the LTE population, to within  ~10% for J = 16−15 and to within a factor of two for J = 30−29. The increase in the fluxes when assuming LTE is explained by the tenuous and hot gas at larger radii also contributing to the emission, as in the case of for example the “warm” CO finger in the outer, upper disk, discussed in Sect. 3.2.

We conclude that different effects govern the shape of the CO ladder including: (i) The temperature of the emitting gas and the density. (ii) Optical depth effects at low-J to mid-J lines. (iii) For optically thin lines, the number of emitting molecules (~N(CO)dΩ) and for optically thick lines, the emitting area (~dΩ). Depending on the transition, different effects dominate. In addition, different temperature, column density, and density structures might reproduce the same CO ladder.

6. Conclusions and outlook

We have used a new detailed physical-chemical model to calculate the atomic fine-structure emission ([C i], [C ii], and [O i]) and CO lines from CO J = 3−2 to J = 30−29 from a protoplanetary disk around a young Herbig Ae star. We have applied the model to the well-studied HD 100546 disk, but the results are more generally applicable. We have compared our results to recent observations of the object and studied the dependence of the line flux on different parameters. Our main results are:

• The high-J emission of CO can only be reproduced when the gas temperature in the inner disk is considerably higher than the dust temperature. A model with the gas temperature set to the dust temperature underpredicts the high-J lines by orders of magnitude (Sect. 3.3). The UV radiation is able to heat the gas sufficiently (Sect. 3.1).

• Both models with a gas/dust  = 100 and δ_C = 0.05 (volatile fraction of carbon δ_C) and models with gas/dust  = 20 and δ_C = 0.5 are able to reproduce the CO ladder and the [O i] lines. Models with a high gas/dust ratio overpredict the upper limit of [C i] but approximately fit [C ii]. Models with a low gas/dust ratio, however agree with the upper limit of [C i] but underpredict [C ii] (Sect. 4.1).

• None of the models reproduce [C i]/[C ii], but models with higher gas/dust ratio have a lower ratio of these lines, closer to the observed value. The line ratio is also smaller in models with larger grains that allow UV radiation to penetrate deeper into the disk (Sects. 4.2, 5.2) and models without PAHs in the outer disk, at radii  >100 AU. The ratio depends less on the stellar radiation field (Sect. 4.4).

• A substantial fraction of the [C ii] emission may emerge from a remnant envelope or a compact foreground (Sect. 5.3), hence we prefer a solution with gas/dust ratio of 100, but a small fraction of volatile carbon (Sect. 5.4). This means that a considerable part of the carbon is locked into more refractory carbon-bearing species and that the carbon partitioning between gas and dust has evolved from diffuse clouds to protoplanetary disk.

• The highest-J emission detected towards HD 100546, CO J = 30−29, emerges from radii of  ~20−50 AU. Mid-J lines, e.g. CO J = 16−15, are emitted at radii of  ~40−90 AU and low-J lines trace the outer disk (Sect. 3.4).

• The high-J lines of CO (J_up > 14) are predicted to be much broader than the low-J lines, which are observable from the ground (Sect. 3.5).

• Changing the dust opacity law to very high values of R_V allows UV radiation to penetrate considerably deeper into the disk. The resulting line fluxes for a fixed gas/dust ratio are much larger. For 0.1–1.5 μm grains, gas/dust  = 10 with δ_C = 0.1 reproduces the observations similarly well as gas/dust  = 100 and δ_C = 0.05, if a R_V = 5.5 dust opacity law is used (Sect. 4.2).

• Assuming a temperature-independent sticking coefficient for the H₂ formation only affects the high-J lines, changing the line fluxes by a factor of a few. The uncertain H₂ formation rate is thus not a critical parameter for the CO emission, except for the highest-J lines detected towards HD 100546 (Sects. 4.5, 5.6).

• The shape of the CO ladder in the rotation diagram is governed by the combination of excitation effects, optical depth, and different emitting regions of the CO lines (Sect. 5.6).

With this paper, we have introduced a new physical-chemical model and applied it to a protoplanetary disk. The model is based on our previous work (Bruderer et al. 2009a, 2010) and calculates the gas-temperature structure in a self-consistent way with the chemistry and the excitation of the molecules. The model has a wide range of applications and we plan to use it in the future to investigate both T Tauri disks as well as the outflow walls of the envelopes of YSOs, and to compare the results with Herschel observations. The model also can be relatively easily adopted to calculate three-dimensional structures and is well-suited to the analysis of upcoming high-resolution ALMA data.

---

Acknowledgments

We are grateful to Neal Evans, Ted Bergin, Gijs Mulders, Olivier Berné, and the DIGIT team for stimulating discussions and scientific support. We thank the anonymous referee and Malcolm Walmsley for carefully reading the manuscript and useful suggestions. We thank Olja Panić for providing the APEX data in electronic form. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Spinoza grant and grant 614.001.008 from the Netherlands Organisation for Scientific Research (NWO), and by the European Community’s Seventh Framework Programme FP7/2007-2013 under grant agreement 238258 (LASSIE). The work on star formation at ETH Zurich is partially funded by the Swiss National Science Foundation grant 200020-113556. S.D.D. acknowledges support by NASA grant NNX08AH28G. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).

References

Online material

Appendix A: Details of the models and benchmark tests

The purpose of this appendix is twofold. First, it presents details of the implementation of the model and second, it compares of results with other codes and benchmark problems. The model is based on those of Bruderer et al. (2009b,a, 2010) and Bruderer (2010), although it contains several improvements. The structure of the model is summarized in Fig. A.1. The modeling process starts with a gas and density distribution e.g. obtained from an analytical prescription. Different modules then calculate the dust temperature and UV radiation (Sect. A.2), abundances of atomic or molecular species (Sect. A.3), and the level population of molecules or atoms acting as coolants (Sect. A.4). The gas-temperature is then calculated from the balance between heating and cooling processes (Sect. A.5). As the abundance and level population enter the heating and cooling rates, abundances and level populations need to be re-calculated until convergence. Finally, spectral lines can be calculated and compared with observations.

To carry out the calculation, the modeled region is divided into individual cells. In this present study, we assume a axisymmetric structure. Since all geometric information is contained in one separate module, the code can in principle also be run for spherically symmetric or three-dimensional structures.

The current study applies the model to a protoplanetary disk, but the model can also be used in the context of protostellar envelopes, for example to calculate the molecular emission from UV irradiated outflow walls. For future applications, we thus present all modules in this appendix, even though they are not used in the current work (such as the calculation of the dust temperature).

A.1. Dust radiative transfer

A.1.1. Dust temperature

The dust temperature in steady state is obtained from the equilibrium between energy absorbed and emitted by dust (A.1)where the dust opacity per unit mass is κ_λ, the incident radiation field is I_λ, and the Planck function for the dust temperature is B_λ(T_dust). The radiation field consists of contributions from both the stellar radiation field and the ambient radiation by dust emission. Thus, different positions in the model are coupled by means of the radiative transfer equation. We neglect the “back-warming” of dust by hot gas and further assume a unique dust temperature for different sizes of dust.

In this present study, we implement the Monte Carlo technique proposed by Bjorkman & Wood (2001): photon packages with an energy distribution following the stellar radiation field are propagated until interaction (scattering or absorption followed by re-emission). Scattering only changes the direction of a photon package. To sample the scattering angle, we use a Henyey-Greenstein function, which defines the scattering angle by only one parameter, g =  ⟨ cos(θ) ⟩ . Absorption increases the amount of energy deposited into a cell. For a cell of dust mass M, the dust temperature T_new after absorption of N_i photon packages is given by (A.2)with the energy per photon package δL = L/N, where the stellar luminosity L is divided by the total number of photon packages N. When an absorption event increases the local temperature by , the photon-package is re-emitted at a wavelength given by (Baes et al. 2005) (A.3)This probability distribution function (PDF) accurately corrects for the incorrect spectrum of previously reemitted photons (at a too low temperature). However, it can sometimes become negative. In such cases, the approximative PDF given in Bjorkman & Wood (2001) is used. This procedure is generally very efficient as it does not require any global iteration and the dust temperature increases to the correct value. In regions of very high optical depth, however, photons can get stuck owing to a very high number of absorption and re-emission processes. We thus implement the diffusion approximation derived in Robitaille (2010), based on Min et al. (2009). Regions of very low optical depth can suffer too little absorption resulting in a noisy dust temperature (see e.g. Bruderer 2010, Fig. D.2). We implement a method proposed by Pinte et al. (2006) and measure the ambient radiation field using the estimator given in Lucy (1999) during the propagation of the photons. After all photons have propagated, the dust temperature can be improved using this value.

To benchmark the dust temperature calculation, we use the benchmark problems proposed in Ivezic et al. (1997) and compare our results with calculations obtained with the DUSTY code (Ivezic & Elitzur 1997). The benchmark problem consists of a spherical cloud with two different density gradients (either  ∝ r^-2 or constant) and four different total radial optical depth at 1 μm (τ_1  μm = 1,10,100, or 1000). The input spectra is a black body and the dust properties follow a simple analytical law. Scattering is assumed to be isotropic. The agreement between the codes is perfect, being closer than about a percent.

As a second benchmark problem, we calculate the dust temperature in a protoplanetary disk with very high total optical depth (Pinte et al. 2009). We calculate their benchmark problem for isotropic scattering with a midplane optical depth at 0.86 μm of τ = 10³ and 10⁵ (their Fig. 9). In Fig. A.3, we show the midplane temperature for the two models compared to the results obtained with MCFOST. The agreement is very good for the model with τ = 10³. For the run with τ = 10⁵, the agreement is again very good with deviations smaller than a few percent, even in the region shadowed by the dense inner rim. These deviations are comparable to those between the other codes participating in that benchmark study.

A.2. UV field

Since the calculation of the dust temperature does not explicitly compute the radiation field as a function of wavelength, the intensity at UV and optical wavelengths is calculated in a second step. This intensity is used for example to calculate photodissociation rates. The calculation of the intensity uses the same method for the propagation of the photon packages as the calculation of the dust temperature. The radiation field is then calculated with the method proposed by van Zadelhoff et al. (2003). For each wavelength bin, N photon packages are propagated. All photon packages are assigned with an initial intensity I_λ(0) = F_λS/N, where F_λ is the input flux and S the surface it passes through. This intensity then drops to I_λ(s + Δs) = I_λ(s)e^−Δτ_λ, when the package advances by Δs with an optical depth Δτ_λ. The mean intensity is calculated from (A.4)with the volume of the cell V and the sum taken over all photon packages passing through a cell. The distance and optical depth as an individual photon package i passes through a cell are given by Δs_i and Δτ_λ,i. Since we calculate the mean intensity for a few wavelengths only, the input intensity is obtained by averaging over the input (stellar) spectra in each individual spectral band. Photons from the interstellar radiation field (Draine 1978) are also accounted for.

As a benchmark test, we calculate the mean intensity in a spherical cloud, switching scattering off. We compare the calculated intensity to the analytically obtained intensity and find agreement to better than one percent. The implementation of the scattering follows the calculation of the dust temperature and was tested in the previous paragraph.

A.3. Chemistry

Abundances of molecular and atomic species are obtained from the solution of the rate equations (A.5)with the abundance n(i,t) (cm^-3) of a species i at time t. The rate of destruction and formation of a species is given by the coefficients k_ij and k_ijl, respectively. In addition, three-body reactions can be entered, but our current network does not contain any. The set of coupled, stiff, and non-linear differential Eqs. (A.5) is solved by the DVODE solver (Brown et al. 1989) for the temporal evolution of the species. We verify that the solution complies with the conservation of elements and charge. We note that the DVODE solver is widely used for applications in astrochemistry and well-tested (e.g. Nejad 2005). To obtain steady-state abundances, Eq. (A.5) is solved for dn(i,t)/dt = 0 using a globally convergent Newton-Raphson method. To close the system of equations Eq. (A.5), some rate equations are replaced by conservation equations. If the steady-state solver fails to converge, the abundances are obtained from the time-dependent solver running to a high chemical age (>10 Mio. years).

To test numerical aspects of the solver, we calculate abundances in an isothermal slab of constant density with UV irradiation. A small chemical network from the PDR comparison study by Röllig et al. (2007) is used for this test. We first calculate their problem F4 (for a density of 10^5.5 cm^-3, UV field of χ = 10⁵, and a temperature of 50 K) with the steady-state solver. In Fig. A.4, we give the abundances of H, H₂, O, O₂, C⁺, CO, and electrons and compare them to the results obtained from PDR codes. Our code to within their scatters agrees with the other codes. The scatter in the H₂ abundance is due to a different implementation of some reactions, such as H₂ dissociation including self-shielding. We also calculated the abundances with the time-dependent solver running to an old chemical age and found an agreement to better than a factor of 10^-4 with the steady-state solution.

A.3.1. Chemical network

The chemical network used in our work contains 10 elements, 109 species, and 1492 reactions. We use the same elemental composition as Jonkheid et al. (2006) and Woitke et al. (2009a), which is summarized in Table A.1. The species contained in the network are given in Table A.2. We note that is vibrationally excited H₂. PAH⁰, PAH⁺, and PAH⁻ are neutral, positively, and negatively charged PAHs respectively and PAH:H denotes hydrogenated PAHs. Species frozen-out onto dust grains are shown by JX, for example JCO is CO frozen-out.

The chemical reaction network is based on the studies by Doty et al. (2004), Stäuber et al. (2005), and Bruderer et al. (2009b), but contains some extensions. The following reaction types are included:

1.) H₂ formation on dust

For the formation of H₂ on dust grains, we implement the rate by Sternberg & Dalgarno (1989), scaled to the appropriate grain size. We use the sticking coefficients of Cuppen et al. (2010). Since the gas/dust ratio can vary over a model, rates involving small dust are scaled to the local gas/dust ratio.

2.) Freeze-out, thermal, and non-thermal desorption

Freeze-out and thermal or non-thermal desorption are implemented following Visser et al. (2009b) and references therein. The thermal desorption rates account for the abundance of the ice on the surface and change from zeroth to first order behavior when the ice contains less than a monolayer of ice. As in Visser et al. (2009b), we use the same pre-exponential factor for all species except for H₂O (Fraser et al. 2001) and CO (Bisschop et al. 2006). Binding energies follow Aikawa et al. (1997) and Sandford & Allamandola (1993).

Non-thermal evaporation by UV photons or cosmic-rays depends on a yield indicating the number of molecules desorbed per grain per incident UV photon. We assume the same yields used by Visser et al. (2009b), following Öberg et al. (2007, 2009b). To simulate cosmic-ray induced desorption, a small background UV flux corresponding to 10⁴ photons cm^-2 s^-1 (Shen et al. 2004) is added to the stellar UV flux.

3.) Hydrogenation of simple species on ices

In addition to freeze-out and evaporation, we implement the grain-surface hydrogenation of simple species. We follow the approach of Visser et al. (2011) and Hollenbach et al. (2009). The hydrogenation leads to saturation of species frozen-out, from JC  →  JCH  →  JCH₂  →  JCH₃  →  JCH₄, JN  →  JNH  →  JNH₂  →  JNH₃ and JO  →  JOH  →  JH₂O. Other grain-surface reactions to build up more complex species are not accounted for, but the focus of our work is on simple species.

5.) Gas-phase reactions

Most reactions in our network (1132 out of 1492) are pure gas-phase reactions, such as ion-neutral, neutral-neutral, or charge exchange reactions. Our network contains all reactions from the UMIST 2006 database (Woodall et al. 2007) involving the species given in Table A.2. Details of the implementation are given in Bruderer et al. (2009b).

6.) Photodissociation

Photodissociation rates are calculated from the local radiation field using the cross-sections for discrete and continuous absorption given by van Dishoeck (1988), van Dishoeck et al. (2006), and van Hemert & van Dishoeck (2008)³. For the dissociation of H₂, CO, and the ionization of C, we include the effects of self-shielding by reducing the unshielded rate with a self-shielding factor. For H₂, the factor given in Draine & Bertoldi (1996) is used. For CO, we implement the factors given in Visser et al. (2009a) and for C, the factor of Kamp & Bertoldi (2000) is used. We however limited the shielding of C by H₂ to 0.5, owing to the overlap of H₂ lines with the wavelengths at which carbon is photoionized.

The self-shielding rates depend on the column densities of H₂, C, and CO towards the UV source. We calculate these column densities together with the UV field (Sect. A.2) by averaging over the densities of all photon packages passing through a cell, weighted by their intensity.

7.) X-ray induced processes

As demonstrated by Bruderer et al. (2009b), the influence of X-rays on the chemistry can be accurately approximated by secondary ionizations induced by fast photo-electrons from H₂ and H ionization. We follow their approach and calculate the X-ray ionization rate using the cross-sections given in their Appendix A. The rates for the secondary processes are taken from Stäuber et al. (2005). We note that the cosmic-ray-induced photodissociation rates are scaled-up to the X-ray ionization rate following Stäuber et al. (2005).

8.) Cosmic-ray-induced reactions

Rates for direct cosmic-ray-induced ionization processes are taken from the UMIST 2006 network. Cosmic-ray-induced photoionization rates are also taken from that database, with the exception of CO, H₂, and H, where we follow Stäuber et al. (2005). For the CO dissociation, we use the rates of Gredel et al. (1987) implemented with the fitting equation provided in Maloney et al. (1996). We account for the ionization of CO, H₂, and H by the decay of excited He (2¹P state, 19.8 eV, see Yan 1997).

9.) PAH/small grain charge exchange/hydrogenation

We assume that PAHs and small grains are well-mixed with the gas and thus not scaled to the gas/dust ratio. For the photoionization, charge-exchange, and recombination, we use the rates given in Wolfire et al. (2003) with a collision rate parameter φ_pah = 0.5. For species heavier than hydrogen, the recombination rates are scaled to , where m_amu is the mass of the molecule.

The formation of CH⁺ and H₂ on PAHs is implemented using the rates of Jonkheid et al. (2006). We include the reactions Reactions (A.6)–(A.8) run at 10^-9 cm³ s^-1, while (A.9) runs at the same speed as the recombination of C⁺ with PAH⁰.

10.) Reactions with (excited state H₂)

Vibrationally excited H₂ is included into the chemistry with a two-level approximation. denotes H₂ in a vibrationally excited pseudo-level with an energy corresponding to 30 163 K (London 1978). The energy stored in the vibrational state can be used to overcome an activation barrier (e.g. Sternberg & Dalgarno 1995). Following Tielens & Hollenbach (1985), the rate coefficients of a species with can be approximated using the rate for H₂ by reducing the exponential factor by 30 163 K. For a H₂ rate  ∝ exp(− γ/T_kin), the exponential factor for is thus γ ^∗  = max(0,γ−30   163K). This approximation however overestimates the rate coefficients, if γ is too large and we switch off reactions with γ > 20   000 K.

The population of is determined by collisions with H and H₂ (Tielens & Hollenbach 1985), spontaneous decay (2 × 10^-7 s^-1; London 1978) and the pumping by FUV photons, which is assumed to be ten times the photodissociation rate.

A.4. Molecular/atomic excitation

The population n_i (in cm^-3) of a level i of a molecular or atomic species is determined by the rate equation (A.10)with the rate coefficients for the transition between level i and j, (A.11)with the level energy E_i, the collisional rate coefficients C_ij and the Einstein A and B-coefficients A_ij and B_ij, respectively. The rate of chemical formation of a molecule forming into the state i and destroyed from that state is given by F_i and D_i, respectively (Stäuber & Bruderer 2009). In the following, we assume steady-state conditions and set dn_i/dt = 0. To calculate the ambient radiation field (A.12)the normalized line profile function φ_ij(ν) and the intensity depending on the frequency and direction is needed in every cell of the model. The intensity along a ray and for one frequency is calculated from the radiative transfer equation (A.13)with the absorption and emission coefficient α_ν and j_ν, respectively and the distance ds that the ray propagates. Both coefficients contain contributions from both line and dust absorption and emission. Since the level population enters the absorption and emission coefficients, Eqs. (A.10), (A.12), and (A.13) form a coupled problem, which is computationally very demanding to solve (see for example Hogerheijde & van der Tak 2000; Brinch & Hogerheijde 2010).

In this work, we employ an approximation similar to the escape probability method (e.g. van der Tak et al. 2007). In this approximation, the physical conditions and thus the source function S_ν = j_ν/α_ν is assumed to be constant in the modeled region. The ambient radiation field  ⟨ J_ij ⟩  can then be expressed by an analytical expression, which greatly facilitates the calculation. For our models with physical conditions strongly varying with position, the approximation of a single zone with constant physical conditions is not applicable. The solution of Eq. (A.13) then reads (A.14)using the definition of the optical depth, dτ_ν = α_νds and the background intensity I_ν,bg. Similar to Poelman & Spaans (2005), we solve for the excitation at every position, although we approximate the local radiation field by keeping the source function constant along a ray. The time-consuming integration in Eq. (A.14) reduces in this way to a simple calculation of the opacity along a ray from the current position to the edge of the modeled region. This approximation allows us to solve for the excitation of complex problems (many lines, high optical depth) in a short time. For example the excitation of water in the models presented in this work can be solved in a few minutes using a standard PC.

To include pumping by dust and the velocity structure, we follow the approach of Takahashi et al. (1983). For the dust and line opacity τ_ν,L and τ_ν,D, along a ray, we calculate The level population is then obtained from Eq. (A.10) with (A.17)and . The local dust temperature is given by T_D and the background radiation (CMB) is given by I_ν,bg. For the calculation of ϵ_ij and η_ij, either 6 or 26 directions and of order 40 frequency bins per line are used.

Once the level population is determined, the molecular cooling rate can be calculated and synthetic maps can be derived by solving the radiative transfer equation (Eq. (A.14)). Our raytracer produces images in fits format. To properly resolve the inner structure (hot core or inner disk), we shoot several rays per pixel through that region and average over the intensity. Along the ray, the integration steps within one computational cell are refined, according to the velocity structure in order to properly resolve narrow lines.

Benchmark test

To benchmark the excitation calculation, we calculate the problem proposed by van Zadelhoff et al. (2002), which consists of the calculation of the HCO⁺ excitation in a spherical, infalling class 0 cloud core. The physical conditions are given in Fig. A.5. The abundance of HCO⁺ is fixed to either 10^-8 or 10^-9 relative to the H₂ density. The difficulty of this benchmark problem lies in the high optical depth reached in the lines and that collision partner densities are below the critical density. The level population is thus far from the local thermal equilibrium and radiative excitation is important.

In Fig. A.6, we compare the solution obtained with this method with the solution calculated by RATRAN (Hogerheijde & van der Tak 2000). That code is widely used and has been well-tested against other codes. For an abundance of 10^-9, the normalized level population agrees to within 30% with the solution provided by RATRAN. An exception is the J = 3 level in the inner part of the cloud. The deviations are larger for an abundance of 10^-8, but still mostly within 30%. The intensity is obtained with the raytracer implemented in our code and convolved to the telescope beam. In Fig. A.6, the J = 1 → 0 and J = 4 → 3 lines are shown for a beam of 29′′ and 14′′, respectively. We also indicate the intensities obtained assuming local thermal equilibrium (LTE) and “thin” conditions, setting the ambient radiation field to 0. As the level populations, the line fluxes mostly agree to within about 30% of the intensity obtained with RATRAN. We note that the intensities obtained with RATRAN are calculated with their raytracer (SKY) and convolved using the MIRIAD package (Qi 2005). Thus, this benchmark problem also tests our raytracer and convolution routine.

The agreement found in with this benchmark problem is similar to the findings of Woitke et al. (2009b), who calculated the water emission from a Herbig Ae disk calculated with a similar method and compared them to RATRAN. We consider the agreement as reasonable and sufficient for our application, since the uncertainties entering through the chemical network calculation are much larger. It also shows that this method is a much better approximation than e.g. assuming LTE.

A.5. Thermal balance

The gas temperature is obtained from the equilibrium between the heating and cooling rates (A.18)where ϵ is the internal energy of the gas and Γ_i and Λ_i are different heating and cooling rates (for example, line cooling or photoelectric heating on dust grains). This non-linear equation is solved by a secant method starting from the dust temperature. If the convergence is insufficient, the process is restarted with a different value. Convergence is reached if the heating and cooling rate agree to within a predefined threshold in all cells (typically 1%). We implement the following heating and cooling rates (see Bruderer et al. 2009a):

1.) Photoelectric heating

Photoelectric heating on large graphite and silicate grains is implemented following Kamp & van Zadelhoff (2001). The required grain absorption cross-section, averaged over the FUV band, is taken to be consistent with the calculation of the UV field and the dust temperature. As for the chemistry, the rates involving dust grains are scaled to the local gas/dust ratio. Photoelectric heating on small dust grains is implemented following Bakes & Tielens (1994). We also implement their recombination cooling rate.

2.) Gas-grain heating or cooling

Depending on the difference between gas and dust temperature, gas-grain collisions can either heat or cool the gas. In regions of very high density, gas-grain collisions couple the gas temperature to the dust temperature (e.g. Doty & Neufeld 1997). We implement the heating and cooling rates following Young et al. (2004).

3.) H₂ heating

Molecular hydrogen can contribute to the heating and cooling in different ways: (i) Line cooling. (ii) Heating through collisional deexcitation of FUV pumped H₂. (iii) Formation heating. (iv) Photodissociation heating. For (i) and (ii), we use the analytical fit of Röllig et al. (2006) to the exact multilevel treatment by Sternberg & Dalgarno (1995). For (iii), we follow Kamp & van Zadelhoff (2001) and assume that one third of the binding energy of  ~4.5 eV is released to the gas for heating. Photodissociation (iv) carries away about 0.4 eV in the form of heat to the gas (e.g. Jonkheid et al. 2004). We note that the H₂ dissociation rate is taken to be consistent with the chemistry.

4.) Cosmic-ray heating

Following Cravens & Dalgarno (1978) and Glassgold & Langer (1973), a cosmic-ray ionization of H₂ (H) releases 8 eV (3.5 eV) to the gas.

5.) X-ray heating

We implement the X-ray heating rates following Dalgarno et al. (1999) for a H₂, H, and He mixture. The energy deposition rate per particle H_x is calculated using the cross-sections provided in Bruderer et al. (2009b).

6.) Line cooling (molecular and atomic fine-structure lines)

Cooling rates of molecular lines and atomic fine-structure lines are derived from the level population calculation (Sect. A.4). We account for O, C, C⁺, CO, ¹³CO, H₂O, and OH. The abundance of ¹³CO is obtained from scaling the CO abundance by ¹²C/¹³C = 70 (Wilson & Rood 1994). Molecular data are taken from the LAMDA database (Schöier et al. 2005) and include the collisional excitation with H, H₂, and electrons for O, C, and C⁺, respectively. For the molecular species, only excitation with H₂ is available.

7.) Hydrogen Lyα and [O I] 6300 Å  cooling

Hydrogen Lyα cooling and cooling by neutral oxygen [O i] by means of the metastable ¹D-³P line at 6300 Å  is included following Sternberg & Dalgarno (1989).

8.) Carbon ionization

Carbon ionization delivers about 1 eV to the gas (Jonkheid et al. 2004). The ionization rate is assumed to be consistent with the chemical network including the effects of self-shielding.

A.5.1. Benchmark test

To test the implementation of the thermal balance calculation, we run the benchmark problem of the PDR comparison study of Röllig et al. (2007), which consists of four slab models with a density of 10³ and 10^5.5 cm^-3 and an incident UV radiation field of χ = 10 or χ = 10⁵ (in units of Draine ISRF, Draine 1978). To eliminate the differences caused by the adopted chemical network, the benchmark study is run with the same simple chemical network consisting of only the 31 species used in the PDR comparison study. The parameters given in Table 5 of Röllig et al. (2007) are implemented.

The gas temperature calculated with our code is shown in Fig. A.7, together with the results from other PDR models. The gas temperature derived with our code is in reasonable agreement with that found by other PDR codes. Differences between the different codes are however considerable, particularly for the high-density (10^5.5 cm^-3) model with high UV intensity (χ = 10⁵). This model is, however, closest to the conditions found in disk atmospheres and outflow walls.

Appendix B: Analytical estimates of fluxes

B.1. Analytical estimate of the [C II] flux

To analytically estimate the [C ii] flux that might originate from a remnant envelope or the foreground, we use (B.1)for the flux F, the integrated intensity I, and the solid angle dΩ. The solid angle of the emitting area is calculated from the projected area A and the distance d to the source. The integrated intensity I is obtained from the line frequency ν, the Einstein-A coefficient A_ul, the column density of C⁺N_C⁺, and the normalized level in the upper level x_u. Assuming the level population in the LTE, which is the case above the critical density of a few times 10³ cm^-3, we reach a maximum x_u = 2/3, given by the statistical weights of the C⁺ levels. These equations assume optically thin emission. For a line width of 1 km s^-1, τ = 1 at the line center is reached for a column density of a few times 10¹⁸ cm^-2.

B.2. Rotation diagram of CO

For the rotation diagram of CO, we consider the optically thin flux of CO, which is assumed to be thermalized, (B.2)using the same notation as in Appendix B.1 and the CO column density N(CO), the upper level energy E_u, and the partition function Q(T). Rearranging yields (B.3)

with the column density per level . Thus, (B.4)We note that the rotation diagram is often given by assuming that F_ul = dΩ_beamI_ul, for the solid angle of the telescope beam , thus the molecules are equally distributed over the beam. This is however inappropriate for the CO ladder probing a wide range of temperatures and subsequently different emitting regions (Sect. 3.4).

How do opacity effects alter the rotation diagram? Assuming a square line profile of (velocity) width Δv for simplicity and neglecting the background intensity, we obtain for the integrated intensity using Eq. (A.14)), (B.5)For τ_ul ≪ 1, , recovering the expression used in Eq. (B.2). For τ_ul > 1, , which is always larger than the optically thin expression. Thus, optical depth effects shift points above the optically thin lines in the rotation diagram.

All Tables

All Figures