Online material

Appendix A: New features in ProDiMo

A.1. Fixed disk structure

In contrast to earlier publications, we use a parametrised description for the shape and distribution of gas and dust in the disk in this paper (A.1)between an inner and outer disk radius, R_in and R_out, respectively, with sharp edges. ρ(r,z) is the local gas mass density. H(r) is the vertical scale height of the disk, assuming to vary with radius as (A.2)H₀ is the reference scale height at reference radius r₀. ϵ is the column density power-law index and β the flaring power. The constant ρ₀ is adjusted such that the integrated disk mass 2π ^! ρ(x,z) dzr dr equals M_disk.

A.2. Dust size distribution and dust settling

The dust grains are assumed to have a power-law size distribution in the unsettled case as (A.3)between minimum and maximum grain radius, a_min and a_max, respectively. The free constant in Eq. (A.3) is adjusted to result in the prescribed unsettled dust/gas mass ratio ρ_d/ρ.

A very simple recipe has been implemented to account for the major effects of vertical dust settling. We assume that the dust grains are distributed vertically with a smaller scale height (A.4)where H(r) is the gas scale height, and s and a_s are two free parameters. Since ProDiMo can self-consistently calculate the vertical stratification of gas from the resulting gas temperatures and mean molecular weights, in which case H(r) does not exist, Eq. (A.4) can not be used directly in the general case. Instead, we make use of the equation that defines H(r), namely c_T = H(r) Ω(r) and write (A.5)where c_T is the local gas isothermal sound speed, is the reduced variant for dust size a, Ω = v_Kepler/r is the angular velocity and v_Kepler the Keplerian velocity. Equation (A.5) is then used to calculate the vertical distribution of dust particles of different sizes with respect to the already determined gas stratification by treating them like an independent fluids with lower c_T as compared to the gas (A.6)where belongs to the gas, and belongs to dust particles in  [a,a + da] . g_z is the local z-component of gravity. The solution of the differential Eq. (A.6) is (A.7)where f′(a,z) da = ρ′(a)/m_d(a) is the settled dust size distribution function, and c_T(z) and f(a,z) ∝ ρ(z) are considered as known. The constant const(a) is then determined to make sure that the vertical column density of dust particles for every size is conserved (A.8)According to these assumptions, all dust quantities that used to be constant, such as the local dust/gas mass ratio ρ_d/ρ, the dust moments  ⟨ a ⟩ ,  ⟨ a² ⟩ ,  ⟨ a³ ⟩ , and the dust opacities per mass, become spatially dependent quantities. The dust absorption and scattering opacities are calculated by applying effective mixing (Bruggeman 1935) and Mie theory, based on f′(a,z).

A.3. PAH ionisation equilibrium and PAH-heating

We consider a typical size of PAH molecules with N_C = 54 carbon atoms and N_H = 18 hydrogen atoms (circumcoronene), motivated by studies that PAHs with much less carbon atoms would not be stable around young stars on timescales of a few Myr. Much larger PAHs are not consistent with the spatial extent of observed PAH emission in various bands (see e.g. Visser et al. 2007). We include PAH⁻, PAH, PAH⁺, PAH²⁺ and PAH³⁺ as additional specimen in the chemical reaction network. Circumcoronene is probably among the smallest PAHs that can survive in a disk around Herbig Ae stars (Visser et al. 2007). The following processes are considered in detail: PAH-photoionisation (Tielens 2005), electron recombination and some charge exchange reactions (Wolfire et al. 2008; Flower & Pineau des Forêts 2003). These reactions do not change the basic PAH lattice, but only affect the charging of the PAH molecules. Hence, the total amount of PAHs is conserved and treated like an element with given element abundance (A.9)where being the standard ISM particle abundance with respect to hydrogen nuclei (Tielens 2008), and f_PAH is the fraction thereof assumed to be present in the disk.

Concerning the photo-ionisation reactions, we consider the PAH absorption cross sections from Li & Draine (2001, see their Eqs. (6) −  with recent updates for the resonance parameters from Draine & Li (2007). We use their neutral PAH cross section for charge k = 0, and the charged cross section otherwise. The photo-ionisation rates  [s^-1]  for PAH molecules with charge k are calculated as (A.10)where ν [Hz] is the frequency, J_ν(r,z)  [erg/cm²/s/Hz/sr]  the mean intensity as computed by the dust continuum radiative transfer, and s_ν(r,z) the following self-shielding factor and are the radial inward and vertical upward hydrogen nuclei column densities [cm^-2] in the disk, as measured from point (r,z). For simplicity, we put . The photo-electron yield is taken from (Jochims et al. 1996) (A.13)According to Weingartner & Draine (2001), the ionisation potentials IP^k depend on PAH-size and charge k. Using their Eqs. (1) and (2) for N_C = 54, the results are IP⁻¹ = 3.10 eV, IP⁰ = 6.24 eV, IP^+ 1 = 9.38 eV and IP^+ 2 = 12.5 eV. The formula is an approximation that reproduces the ionisation potential of benzene as the smallest PAH and graphite as infinitely large PAH. The threshold wavelengths are given by .

With these photo-ionisation and recombination rates, the local particle densities of PAH⁻, PAH, PAH⁺, PAH²⁺ and PAH³⁺ are calculated consistently with the gas-phase and ice chemistry, with a considerable influence on the resultant local electron density. Once these particle densities have been determined, the total PAH heating rate  [erg/cm³/s]  can be calculated as (A.14)The PAH recombination cooling rate is calculated as (A.15)where is the PAH recombination rate coefficient, n_e is the electron density and the gas temperature. As a general result, we obtain a top-down layered PAH charge structure, with PAH²⁺ in the uppermost layers with reduced PAH-heating, but PAH⁻ and intensified PAH-heating close to the midplane. For massive disks (unlike ET Cha), where the electron density virtually vanishes in the deepest layers, the PAH charge becomes neutral again in the midplane.

A.4. Fluorescent UV-pumping

We have replaced the small model atoms for O, C and C⁺ (as listed in Table 4 of Woitke et al. 2009) by larger ones, see Table A.1 (this paper), to account for fluorescent UV and optical pumping of the lower levels responsible for the far IR fine-structure lines. The new model atoms have collisional data only among the lowest states, but much more radiative data. The additional energy levels and transition probabilities are taken from the National Institute of Standards and Technology atomic spectroscopic database (Ralchenko 2009). We selected all evaluated transitions longward of 912 Å. Rovibronic transitions are permitted and their probabilities are of the order of (1−10⁶) s^-1, much higher than typical collision rates (~10^-9 s^-1).

Collisional data connecting the two first electronic-excited states at 22830 K and at 48370 K of neutral oxygen and the three spin-orbit split ground states ³P exist with electrons and hydrogen atoms as collision partners (Störzer & Hollenbach 2000; Krems et al. 2006). For C and C⁺, we use the rates collected in the Leiden Lambda database (Schöier et al. 2005), whose collision rates are mostly limited to the ground state electronic levels.

A.5. [OI] 6300 Å -pumping by OH photo-dissociation

Photo-dissociation of OH by absorption of UV photons into the 1 ²Σ⁻ state, λ ≃ (973−1350) Å, produces an electronically excited oxygen atom O(¹D) as (A.16)while absorption into the 1 ¹Δ and 3 ²Π states, λ ≃ (1100−1907) Å, leads to the ground state ³P. According to van Dishoeck & Dalgarno (1984), about 55% of the OH photodissociations by an interstellar UV field⁶ result in an O-atom in the state, which happens to be the upper level of the 6300 Å line . This chemical pumping was suggested by Acke et al. (2005) to play a major role in the formation of the 6300 Å line around Herbig Ae/Be stars. We take this effect into account by introducing a quasi-collisional pumping rate  [s^-1]  (A.17)where is the OH photodissociation rate  [1/s] . We apply Eq. (A.17) to l =  { 1,2,3 } , where 1 = ³P₂, 2 = ³P₁, 3 = ³P₀, and u = 4 = ¹D₂. We add these excitation rates to the other collisional excitation rates C_l u (see Sect. 6.1 in Woitke et al. 2009). In writing Eq. (A.17), we rely on statistical equilibrium and assume that the majority of chemical reaction channels forming OH start from the 3 lowest ³P levels with n_O ≈ n_O(1) + n_O(2) + n_O(3), making sure that , where are the fractional populations.

A.6. H₂-pumping by its formation on dust surfaces

The formation of H₂ on grain surfaces liberates the binding energy, causing the newly formed H₂ molecules to leave the surface in a vibrationally highly excited state. Duley & Williams (1986) estimated that the H₂ molecules will be released in vibrational states as high as v ≲ 7, but in the rotational ground state J = 0 for p-H₂ and J = 1 for o-H₂. (A.18)We treat this formation-pumping by (A.19)where is the H₂ formation rate on dust surfaces  [1/s] , u is the index of the highest state with v ≤ 7 and J ≤ 1 that still has collisional de-excitation rates in the database, and Eq. (A.19) is applied to all states l < u, i.e. to all states that have lower excitation energy than u. n_tot is the total ortho or para H₂ particle density, respectively, and α = n_o/p − H2/n_H2 is the prescribed ortho/para-H₂ fraction.

Our formulation of the H₂-pumping, therefore, leads to a constant de-population of all states l < u with a certain timescale (simulating the disappearance of H₂ in all states due to other chemical channels that are feeding the H₂ formation) and an equally strong population  [cm^-3 s^-1]  of the state u due to the formation of H₂ on grains.

The condition of state u being collisionally connected avoids artifacts at very low densities, where collisions are rare and the pumping would lead to an almost complete de-population of all low-lying states. According to our current collection of collisional H₂ de-excitation rates, this condition implies u = (v = 3,J = 0/1). We checked that our modelling of the o-H₂v = 1 → 0 S(1) line at 2.122 μm does not depend much on the choice of this upper level u as long as v > 1.

A.7. Line transfer

During solving the chemistry and energy balance of the gas, ProDiMo calculates the various level populations of atoms, ions and molecules by means of an escape probability method (see Sect. 6.1.1 in Woitke et al. 2009). These are stored and used later to perform a formal solution of line transfer for selected spectral lines.

The continuum  +  line radiative transfer equation is given by (A.20)where I_ν is the spectral intensity and s the distance on a ray. The source function and extinction coefficient are given by Assuming isotropic dust scattering, the continuum and line transfer coefficients are given by where and are the local dust absorption and scattering coefficients, B_ν is the Planck function, J_ν is the local mean intensity (taken from the results of the continuum radiative transfer), ν is the line centre frequency, n_u and n_l are the level populations  [cm^-3]  in the upper and lower state, respectively, and A_ul and B_ul,B_lu are the Einstein coefficients.

The profile function φ_ν  [Hz^-1]  is assumed to be given by a Gaussian with thermal  +  turbulent broadening where v′ is the observers velocity with respect to the star along backward ray direction n, v is the 3D-velocity of the emitting gas with respect to the star (assumed to be given by Keplerian orbits), and x the local velocity shift. The velocity width is and the thermal velocity where m is the mass of the line emitting species.

We use the same setup of parallel rays as described in Thi et al. (2011, see their Sect. 2.3), organised in about 150 log-equidistant concentric rings in the image plane, each subdivided into 72 angular segments. Equation (A.20) is solved numerically on each of these rays in the observer’s frame on 151 velocity grid points to sample v′, where the transport coefficients are pre-calculated on the grid points and later interpolated along the rays, whereas the profile function is always calculated from scratch. The numerical scheme features a variable step size which is controlled by comparing the results after two consecutive steps with the results obtained after one step of double size.

A.8. Chemical heating

By definition, exothermic chemical reactions convert chemical potential energies into heat, whereas endothermic reactions consume internal kinetic energy and actually cool the gas. We calculate this chemical heating/cooling rate  [erg/cm³/s]  as (A.29)where r is an reaction index, R(r) is the reaction rate  [1/cm³/s]  and ΔH_r  [erg]  is the reaction enthalpy, which is positive for exothermic reactions and negative for endothermic reactions (A.30) [erg] is the heat of formation of the chemical species involved in the reactions (pr means products, ed means educts). By simplification, we neglect the temperature-dependence of ΔH_f, and take the values for the heat of formation for all species from (Millar et al. 1997), who list [kJ/mol] at 0 K in their Table 2.

The details of exothermic reactions are quite difficult and often not precisely known. The excess chemical binding energy is seldomly released in form of kinetic energy directly (although there are some exceptions like, for example, dielectronic recombination which is radiationless). But in the vast majority of cases, the reaction will create products in some kind of electronic, vibrational and rotational excited states.

The energy temporarily stored in these excitational states can then be either radiated away (in which case there is only little net heating), or subsequent collisional processes can thermalise the excess energy. The ratio between these two competing processes depends on the critical density for collisional de-excitation, which strongly depends on type and amount of excitation, for example  ~10¹⁵ cm^-3 for permitted electronic transitions down to  ~10³ cm^-3 for rotational transitions. Further complications arise in possibly high optical depth where created photons cannot escape directly, scatter around and drive secondary processes, with a higher probability to get (partly) thermalised. To avoid all these complications, we have simply introduced an efficiency in Eq. (A.29). We also exclude certain types of reactions from Eq. (A.29): reactions which are known to produce or absorb photons, cosmic ray and cosmic ray-induced reactions, X-ray primary and secondary reactions, reactions on grain surfaces, and reactions which are energetically treated in larger detail elsewhere. Some of the most important reactions are found to be where we put for the last two “collider” reactions, assuming that these highly endothermic reactions, which only occur in hot gases, are in fact driven by the kinetic energies of the colliding atoms and molecules.

	Fig. A.1 Sketch of energetics in astrochemical reaction cycles.
Open with DEXTER

In kinetic chemical equilibrium, as assumed in this paper, there is no net formation/destruction of any chemical species, but the chemistry is organised in complicated reaction cycles as sketched in Fig. A.1. In the absence of cosmic rays, UV photons and X-ray irradiation, only the energetically most favourable chemical configurations, like H₂,H₂O,CO,CH₄,NH₃ etc., are abundant. However, due to the impact of these high-energy photons and particles, much less stable molecules, atoms and ions are continuously created which, via some complicated chemical paths, come cascading down in energy again, eventually re-forming the abundant stable molecules. Along these paths, some of the chemical potential energy can be lost in form of secondary photons. The net effect of this cycle is typically a heating, because we have – by definition – excluded the primary UV, X-ray and CR reactions from Eq. (A.29). The chemical heating is hence yet another way to partly thermalise the energy of incoming high-energy photons and particles through secondary exothermic reactions.

To our surprise, the chemical heating, even with low efficiency , results to be an important heating process in protoplanetary disks, in particular at the bottom of the warm molecular layer where many of the near-mid IR spectral lines are formed, and densities are of order 10⁹−10¹¹ cm^-3, which, as far we are aware of, has not been noticed so far in other disk modelling papers.

Appendix B: The outflow model

We have used the HVC [OI] 6300 Å and [SII] 6731 Å line luminosities to estimate the outflow mass loss rate in two independent ways as suggested by Hartigan et al. (1995, see their Eqs. (A8) and (A10)). Assuming the same values for electron density (only necessary for the OI line) and the projected velocity of the jet as used by in Hartigan et al., we obtain an outflow mass loss rate of Ṁ_outflow = (0.9−2) × 10^-9   M_⊙/yr, similar to the mass accretion rate as estimated by (Lawson et al. 2004).

Hollenbach (1985) showed that it is possible to relate the mass outflow rate to the [OI] 63.2 μm line luminosity by assuming that each atom in the outflow passes through a single shock wave and that the gas with T_g < 5000 K cools by radiating only in the [OI] 63.2 μm line. This procedure is likely to overestimate the resulting line flux since other cooling processes can be active as well, and since the jet material typically passes through several shock waves close to the star. Nevertheless, we use the formula derived by Hollenbach (1985) to obtain a rough estimate of the [OI] 63.2 μm line luminosity as would be expected for an Ṁ_outflow = (0.9−2)  ×  10^-9   M_⊙/yr outflow. The result is a [OI] 63 μm line luminosity equal to (0.6−1.5) times the observed value. This estimate suggests that a substantial part of the [OI] 63.2 μm line as seen by Herschel might originate in the outflow rather than in the disk. We note, however, that the peak of the HVC is shifted only by about −42 km s^-1 with respect to the stellar velocity, and the formulae used above may not be appropriate to such unusually low outflow/shock velocities.

Another way to estimate the [OI] 63.2 μm emission line flux from an outflow is to use the shock models of (Hartigan et al. 2004), which predict the line ratio [OI] 6300 Å/[OI] 63.2 μm for a variety of outflow shock parameters. The authors find 6300/63.2 line ratios of about (0.8−2). For ET Cha, the measured HVC [OI] 6300 Å line flux is (37−87) × 10^-18 W/m², hence the predicted [OI] 63.2 μm outflow line flux should be (30−175) × 10^-18 W/m², which is just consistent with the line flux as observed with Herschel, (30.5 ± 3.2) × 10^-18 W/m². If we include the LVC and work with the total [OI] 6300 Å emission line flux, the results become inconsistent, i.e. the outflow model alone would already over-predict the measured flux.

Appendix C: Details of the best-fitting disk model

	Fig. C.1 Dust opacities assumed in the best-fitting disk model, calculated with effective mixing and Mie theory according to the parameters listed in Table 5. Note that the absorption coefficient (red dotted) is higher than the scattering (black dashed) also at optical and near IR wavelengths, which is untypical for pure silicates.
Open with DEXTER

Figure C.1 shows the dust opacities assumed in the best-fitting model, which results to be absorption-dominated and to scale roughly like λ, except for the 10 μm and 20 μm silicate features.

Fig. C.2 Physical details of the best-fitting disk model. Upper row: total hydrogen nuclei density n ⟨ H ⟩   [cm-3] , with overplotted contours for visual extinction AV = 1, and gas temperature Tgas = 1000 K, and strength of UV radiation field with respect to interstellar standard χ. Second row: dust and gas temperature structures, Tdust and Tgas (note the different scaling). Lower row: most important heating and cooling processes.

Open with DEXTER

Figure C.2 visualises the densities, UV radiation field strength, and resulting gas and dust temperatures of the best-fitting model. We obtain the typical T_dust-pattern for passive disks (see e.g. Pinte et al. 2009) with a shadow casted by the innermost disk regions and dust temperatures of the order of 1450 K at the inner rim at 0.022 AU, and  ~35 K at the outer disk edge located at 8.2 AU. The gas temperature structure T_gas(r,z) follows the dust temperature structure in the midplane (where thermal accommodation dominates) but shows substantial deviations at higher altitudes where the disk starts to become optically thin to UV radiation.

In particular, there is a “hot finger” at relative heights z/r ≈ 0.3−0.6 in this model, stretching out from the inner rim to about 4 AU in radius. Here, gas temperatures of the order of 5000 K are achieved in this model due to PAH heating versus FeII and MgII line cooling. Just below the hot atomic layer, there is a warm molecular layer with temperatures 300−1500 K, heated by PAHs and chemical heating, balanced by CO ro-vibrational cooling. The outermost layers beyond 4 AU are featured by an equilibrium between PAH and chemical heating, versus OI, OH, H₂O and CO rotational line cooling.

	Fig. C.3 Chemical details of the best-fitting disk model. Upper row: electron concentration nel/n ⟨ H ⟩  and mean PAH-charge, with overplotted contours for ionisation parameter log (χ/n ⟨ H ⟩ ). Second row: neutral carbon and CO-concentration, the CO-plot includes contours for dust temperature (blue) and gas temperature structure (white). Lower row: OH and H2O concentration.
Open with DEXTER

Figure C.3 shows some details about the chemical structure of the disk. There is a quite sudden transition of the charge of the PAHs from 3 to –1 where the ionisation parameter χ/n _⟨ H ⟩  is about 100. The PAHs are mostly negatively charged then, except for the outer midplane where the gas almost completely neutralises. The neutral carbon, CO, OH and H₂O concentrations are typical of a photon dominated region (PDR-structure).

Fig. C.4 Density structure of inner rim and impact on SED. The upper left plot shows the prescribed density structure of our best-fitting model. Dashed contour lines in the upper panel refer to the calculated dust temperatures. The upper right plot represents a model with identical parameters where, however, the vertical disk structure is calculated consistently with the resultant gas temperatures and mean molecular weights, which results in a flatter midplane close to the star and too little near-mid IR excess (lower right plot). Note the “soft inner edge” (see Woitke et al. 2009). The lower left plot compares the prescribed scale height of the best-fitting model (full line) with the scale heights resulting from the model with calculated vertical disk stratification. The red dotted line shows these results in a deeper layer at z/r = 0.05, and the blue dotted line at z/r = 0.5.

Open with DEXTER

Appendix D: Variation of input physics

D.1. Difference between dust and gas temperature

A control model where the gas temperature is assumed to be equal to the dust temperature resulted in the line fluxes shown in Table D.1. We conclude that modelling the gas energy balance, mostly leading to T_gas > T_dust in the line forming regions, is absolutely essential to understand the gas emission lines from protoplanetary disks.

D.2. Influence of X-rays

We have run a comparison disk model with X-ray heating and chemistry included, as has recently been implemented by Aresu et al. (2011), with X-ray luminosity L_X = 6 × 10²⁸ erg/s (XMM observations by López-Santiago et al. 2010). We assumed an X-ray emission temperature of T_X = 10⁷ K and a minimum energy of X-ray photons of 0.1 keV. This model does not result in any observable changes in the calculated line fluxes. The modification by X-rays are only –0.5%, 1.1%, 1.4% and 0.05% for CO J = 3 → 2, [OI] 63 μm, [OI] 6300 Å (LVC) and o-H₂ 2.122 μm, respectively. Since the X-rays are attenuated by gas in the model, but the UV photons by dust, and our best-fitting model is extremely gas-rich (assumed gas/dust mass ratio  ~23 000, see Table 5), the X-rays do not penetrate deep enough to change the energy balance in the line emitting regions.

D.3. Influence of viscous heating

We have run a comparison disk model with viscous (gas) heating included, via the formula of Frank et al. (1992)(D.1)where ρ is the gas mass density, ν_kin = α c_T H_p the viscosity, α the viscosity parameter, the isothermal sound speed, H_p = c_T/Ω_kep the pressure scale height, and the Keplerian angular velocity. Putting the viscosity parameter to α = 0.01, we find no significant changes in the far-IR and (sub-)mm lines, but modest changes in the calculated line fluxes in the optical and near-IR. The enhancement by viscous heating is 1.0%, 2.4%, 25% and 16% for CO J = 3 → 2, [OI] 63 μm, [OI] 6300 Å (LVC) and o-H₂ 2.122 μm, respectively. Since the viscous heating scales as Γ_vis ∝ ρ, but most cooling processes scale as Λ ∝ ρ², the effect of the viscous heating is actually strongest in the low density uppermost disk regions, which is counter-intuitive. It renders the gas temperature in these layers unbound (artificially limited by 20 000 K), as no implemented cooling process is able to balance the viscous heating according to Eq. (D.1) in a thin gas. We therefore refrain from discussing the effects of viscous heating any further in this paper.

D.4. Self-consistent disk structure

Figure C.4 visualises the density structure as resultant from a self-consistent disk model where the vertical disk stratification is a result of radiative transfer, chemistry, and gas energy balance, assuming vertical hydrostatic equilibrium. We find a good match concerning the flaring angle (the slope in the lower left plot), but a generally flatter disk structure, more condensed toward the midplane. The z-dependent scale heights from the self-consistently calculated disk model are calculated as (D.2)where p denotes the gas pressure. In the lower left plot of Fig. C.4 we have plotted two scale heights from the self-consistent disk model, one measured close to the midplane (at relative height z/r = 0.05, red dotted line) and one measured high above the midplane (at z/r = 0.5, blue dotted line). The difference between these two scale heights is a natural consequence of our calculated gas temperature structure, with cold conditions in the midplane and a warm/hot disk surface.

We observe a fair match of the scale heights at z/r = 0.5 with the prescribed scale-heights from our best-fitting model, but in the midplane the scale heights of the best-fitting model are about a factor of 2−3 too large. Figure C.4 also demonstrate that, consequently, we loose our SED-fit when using the self-consistent model. The self-consistent model intercepts less star light, resulting in a significant cooling and flux deficit in the near and mid IR as compared to the observations.

D.5. Influence of non-radiative dust heating

	Fig. D.1 Comparison between predicted SEDs of two models with and without non-radiative dust heating. The black model shows once more the SED of our best fitting model, without non-radiative dust heating. In the green model, non-radiative dust heating through thermal accommodation is included. All model parameters are identical otherwise.
Open with DEXTER

In our best-fitting model, we have ignored non-radiative heating/cooling of the dust grains when determining the dust temperature structure T_dust(r,z). However, since the gas is typically warmer than the dust, inelastic gas-grain collisions (thermal accommodation) lead to a collisional, i.e. non-radiative, heating of the dust. If we include this effect (see Eqs. (14) and (108) in Woitke et al. 2009), we do not observe much of an effect on the calculated line fluxes, but the dust temperatures result to be slightly higher, with noticeable effects on the SED, see Fig. D.1.

According to this model with the dust in non-radiative equilibrium, the temperature contrast between gas and dust is mainly driven by exothermic chemical reactions which are active even in quite deep and dense layers (see Fig. C.2), causing an overall heat transfer from gas to dust in the disk of  ~4.5 × 10^-3   L_⊙, i.e. about 5% of the stellar luminosity. It is this additional energy input that leaves the disk in form of additional mid IR continuum photons, causing the depicted variations in Fig. D.1.

The effect of non-radiative dust heating on the SED is similar to increasing the scale heights. We have performed an additional run of the evolutionary strategy with enabled non-radiative dust heating. This run did not entirely converge. The final parameter set had a gas mass of 1.2 × 10^-3   M_⊙, a scale height of only 0.007 AU, i.e. a reduction of  ~35% with respect to our best-fitting model. Unfortunately this is a slow and unstable option, because there is an additional outer iteration necessary between gas and dust temperature determination, to achieve consistent results, which requires about 3 times more computational time.

D.6. Influence of treatment of H₂-formation

The formation of H₂ on dust grain surfaces is one of the most important first steps to initiate a rich molecular chemistry. It has profound effects also on other abundances, for instance the C⁺/C/CO transition (because of the mutual H₂/C shielding) and on the formation of OH and H₂O. Yet its rate is still rather uncertain. Our default choice is to calculate this key chemical process according to (Cazaux & Tielens 2004, called “model B” in Table D.2). We have run two comparison models with two other formulations, one with a typically lower H₂-formation rate according to (Sternberg & Dalgarno 1995). And one with a typically higher rate according to (Cazaux & Tielens 2010).

Concerning the formulation of Sternberg & Dalgarno (1995), which is valid for standard ISM size distribution and dust/gas ratio only, we add a scaling factor to account for deviations of the total dust surface per hydrogen nucleus in disks as (D.3)

The H₂-formation rate coefficient needs to be multiplied by the neutral hydrogen atom density n_H to get the H₂-formation rate in  [cm^-3 s^-1] . The normalisation factor in Eq. (D.3) results from  ⟨ a² ⟩ _ISM  [n_dust/n _⟨ H ⟩ ] _ISM = 5.899 × 10^-22   cm² under interstellar conditions, i.e. for a_min = 0.005 μm, a_min = 0.25 μm, p = 3.5, ρ_gr = 3 g/cm³ and ρ_dust/ρ_gas = 0.01.

Table D.2 shows a strong dependence of the predicted o-H₂ 2.122 μm line on the assumed H₂-formation rate on grains. The Sternberg & Dalgarno (1995)-formalism results in a factor of 0.77 weaker and the Cazaux & Tielens (2010)-formalism in a factor of 4.8 stronger line flux. This is a daunting example of hidden uncertainties in astrochemical modelling. We note that the Sternberg & Dalgarno (1995)-formalism gives a lower FWHM ≈ 30 km s^-1 that brings this model closer to the observations. All other calculated gas emission lines (including those of H₂O) are less affected.

© ESO, 2011