© ESO, 2011

1. Introduction

In the recent decade, modelling of the gas chemistry and physics in protoplanetary discs has become a field of its own, following the lead of the dust SED modelling with multi-dimensional continuum radiative transfer (e.g. Chiang & Goldreich 1997; D’Alessio et al. 1998; Dullemond et al. 2002). The impact of the gas chemistry and energy balance is greatest in the disc surface layers where most of the infrared and sub-mm line emission arises (e.g. Kamp & Dullemond 2004; Jonkheid et al. 2004; Gorti & Hollenbach 2004; Nomura & Millar 2005). The complexity of the gas chemistry and individual heating/cooling processes, as well as their mutual dependence, makes gas structure modelling more complicated, physically more demanding, and computationally more time consuming. As a result, most modelling has so far concentrated on either single selected objects (e.g. Qi et al. 2003; Semenov et al. 2005; Qi et al. 2008; Gorti & Hollenbach 2008; Henning et al. 2010), prototypes (e.g. Aikawa et al. 2002; Meijerink et al. 2008; Goicoechea et al. 2009), or small model series studying the dependencies of specific parameters on the gas emission from discs (e.g. Aikawa & Nomura 2006; Jonkheid et al. 2007; Nomura et al. 2007; Glassgold et al. 2009; Kamp et al. 2010).

Current and upcoming observing facilities such as Herschel, ALMA, JWST/MIRI, SPICA/SAFARI are excellent tools for studying the gas in large samples of protoplanetary discs. There is therefore a clear need for a simultaneous and consistent modelling of dust and gas in protoplanetary discs, with the aim to identify the diagnostic power of certain continuum and line observations for the physical, chemical and temperature structure in the discs. These identifications must be based on a wide variety of disc models covering a wide spread in stellar, disc, dust, and gas parameters.

We jointly produced a grid of 300   000 disc models (Disc Evolution with Neat Theory – DENT) to study the predictive power of individual gas lines, line ratios, and continuum tracers (Woitke et al. 2010). The line emission studies concentrate on a series of fine structure (O i and C ii) and molecular lines (high J CO and H₂O) that are observable with Herschel, as well as on a number of CO sub-mm lines accessible from the ground. A first analysis of the grid has shown that, for example, the [O i] 63 μm emission strongly depends on the amount of UV excess coming from the star (Woitke et al. 2010). Subsequent application to the first findings from the Herschel open time key program GASPS (Pinte et al. 2010) has revealed that the measured [O i] fluxes for a sample of 30 discs indicate that discs around T Tauri stars indeed require extra heating of the disc surface – either through a UV excess and/or X-ray irradiation – while the emission from discs around Herbig Ae stars can generally be matched by photospheric heating alone. This paper now presents a more systematic and complete study of the gas properties and diagnostics from the DENT grid. The dust diagnostics will be discussed in a separate paper.

In the following, we briefly summarize the DENT grid approach (Sect. 2) as presented in Woitke et al. (2010) before we discuss the statistics of gas chemistry and thermal structure of the 300   000 models in Sect. 3 and the line diagnostics in Sect. 4. We present our current best strategies for deriving disc gas parameters in the final conclusions (Sect. 6).

2. The DENT grid

The DENT grid has been described for the first time in Woitke et al. (2010) and we summarize here only a few key aspects. More details can be found in the Appendix.

The grid is computed using the 3D Monte Carlo radiative transfer code MCFOST (Pinte et al. 2006, 2009) and the gas chemistry and physics modules of ProDiMo (Woitke et al. 2009a; Kamp et al. 2010). The column densities, flaring, dust size distribution and settling are parametrized with power laws, assuming a prescribed global gas-to-dust ratio. Mie theory is used to compute the dust optical properties assuming an astronomical silicate composition. The central star is parametrized by its mass and age, as well as strength of the UV excess. Dust and gas temperatures are calculated separately from Monte Carlo continuum radiative transfer and a detailed heating/cooling network, respectively (see Sects. 2.1 and 2.2).

We have chosen 11 free model parameters for the grid, leading to a total of 322   030 disc models. The range of parameter values explored in DENT is discussed further in the Appendix. The total computing time of the grid is 200 000 CPU hours, an average of 40 CPU minutes per model. The total grid took 3 weeks on 400 processors and generated 1.8 TB of data.

The grid is very coarse in nature, e.g. stepping through gas mass in factors of ten. Hence, we cannot expect any method that we develop to measure disc properties such as the disc gas mass from gas line emission or line ratios to be more accurate than that step size, in this case an order of magnitude. The same holds for other quantities such as the outer disc radius, the flaring index, the dust-to-gas mass ratio etc.

For the following analysis, it is also important to keep in mind that the grid of disc models is not constructed to reflect the statistics with which discs occur in nature, e.g. the fact that some objects are rarer than others. There could be some disc types that we missed in our approach, or others that simply do not occur in nature. We take the latter into account by excluding certain parameter combinations in the analysis of the grid. Models with a flaring index of β = 0.8 are for example more appropriate for late stages of disc evolution such as debris discs and so they are excluded from the study of line diagnostics in younger discs (noted in the figure captions whenever appropriate). A global disc analysis using the Toomre criterium shows that 94% of the disc models in the DENT grid are gravitationally stable (see Appendix A.3).

We chose parametrized density structures to avoid imposing the complex gas structure that is based on presumptions mainly derived from disc modelling approaches such as those of Gorti & Hollenbach (2004), Nomura & Millar (2005) or Woitke et al. (2009a). While it is true that the more consistent vertical disc structure models need less assumptions than the simple parametrized models used in DENT, even a full ProDiMo including the vertical structure calculation (Woitke et al. 2009a) is not guaranteed to capture all relevant physics. Using parametrized disc structures enables us to explore a wide unbiased parameter space. From these simplified disc structures, we can still assess the relative influence of several of the key parameters on line and continuum fluxes. Even though the modeled fluxes are not absolute, we can still test fundamental aspects of our disc understanding, such as uncertainties in various observational techniques to derive disc dust and gas masses.

2.1. MCFOST

MCFOST is a 3D Monte Carlo continuum and line radiative transfer code. The details of the numerical schemes are presented in Pinte et al. (2006) and Pinte et al. (2009). In short, the code stochastically propagates photon packets in 3D through the disc. The transport of packets is governed by successive scattering, absorption and re-emission events, which are determined by the local dust properties (opacity, albedo, scattering matrix) and temperature. MCFOST uses the immediate re-emission concept (Bjorkman & Wood 2001) with a continuous deposition of energy to estimate the mean intensity (Lucy 1999). In the optically thick regions of the disc, the code uses a diffusion approximation method to converge the temperature structure. Emerging SEDs are generated by a ray-tracing run from the temperature and radiation field estimated by the Monte Carlo run.

For the DENT grid, 10⁶ packets were used to compute the temperature structure. The radiation field was converged using 10⁴ packets per wavelength bin at wavelength greater than 0.5 μm. As a reliable estimate of the UV field is critical for the chemistry module in ProDiMo, we used 10⁵ packets per wavelength bin for wavelengths shorter than 0.5 μm. This ensures that the UV radiation field is well converged, even in the deep regions of the disc.

MCFOST also includes a NLTE line transfer module which is not fully used here, but whose ray-tracing module is used instead to compute the emerging spectral lines from the level populations estimated by ProDiMo.

2.2. ProDiMo

ProDiMo calculates radiation thermo-chemical models of protoplanetary discs. It is described in detail in Woitke et al. (2009a) and Kamp et al. (2010) and we summarize here the most salient features. Most noticable, the code is not used in its full capacity, because the self-consistent computation of the vertical hydrostatic equilibrium structure is omitted.

The code takes as input the dust density and temperature, n_dust(r,z), T_dust(r,z) and dust properties as provided by MCFOST. In addition, the mean radiation field J_ν(r,z) from MCFOST is passed on to ProDiMo. There are two grid parameters that apply only to the gas part of the model: the gas mass M_gas, and hence the dust-to-gas mass ratio δ, and the fractional UV luminosity defined as f_UV = L_UV/L _∗  in the wavelength range between 91 and 250 nm.

The chemical network consists of 9 elements, 71 species connected through 950 reactions (neutral-neutral, ion-molecule, photoreactions, cosmic ray reactions and adsorption & desorption of CO, CO₂, H₂O, NH₃, CH₄). The background radiation field J_ν enters the computation of the gas chemical and energy balance at several points, namely in the detailed frequency dependent integration of the UV ionization and dissociation cross sections, in the radiative pumping rates involved in the non-LTE modelling of atoms and molecules (O i, C i, C ii, Mg ii, Fe ii, Si ii, S ii, CO rotational & ro-vibrational, o-H₂ & p-H₂ ro-vibrational, o-H₂O & p-H₂O rotational), and in the photoelectric heating rates. The atomic and molecular data used in the statistical equilibrium calculation is described in detail in Sect. 6.1.5 of Woitke et al. (2009a).

3. Gas physics and chemistry

As a first step in the statistical analysis of the DENT grid, we discuss here how the masses and average temperatures of a number of key species (O, CO, C⁺ and H₂O) depend on the grid parameters (see Appendix A.4). With this method, we can obtain a first understanding of the disc chemistry (atomic versus molecular), the relevance of ices (how much CO or H₂O is bound in ice form) and the thermal disc structure (gas dust de-coupling and cool surface layers).

3.1. Oxygen and CO masses

Oxygen is the third most abundant element in the interstellar medium. In dense molecular clouds or protoplanetary discs, oxygen can be either in the gas or ice phase (both are included in our chemical model); in the gas phase, it can be atomic (O) or molecular (CO, O₂, H₂O). By measuring all possible phases, we can constrain the total oxygen mass in a disc.

Figure 1 shows the statistical distribution of atomic oxygen masses. Low mass disc models have almost no CO, H₂O or ices (X(CO) ≲ 10^-3, X(H₂O) ≲ 10^-4) and hence their atomic oxygen mass fraction is close to maximum. Assuming all oxygen is atomic, we obtain a mass fraction of X(O) = 2.3  ×  10^-3. Higher mass disc models show a wider distribution of atomic oxygen mass fractions in Fig. 1. This is due to a higher molecular fraction (dust shielding and gas self-shielding), e.g. CO, CH₄, and H₂O. Studying the CO mass fraction distribution (Fig. 2), the number of discs with masses greater than 10^-5   M_⊙ that have a particular CO mass fraction is rising all the way towards the maximum possible CO mass fraction. The lowest mass discs (10^-6 ≤ M_gas ≤ 10^-7   M_⊙) show instead a broader and flatter distribution; ice formation and photodissociation for those models are extremely model dependent. Also, none of those models ever shows all carbon being bound in CO. The dust-to-gas mass ratio lowers for example the peak of the CO mass fraction by a factor 5 going from 0.001 to 1 (Fig. 3). This is due to the average dust temperature being lower in the high dust-to-gas models and hence more CO freezing out.

3.2. Water masses

The mass fraction of gas phase water depends strongly on the surface density power law index, i.e. on the distribution of mass within the disc. Massive disc models that concentrate the mass in the inner dense disc (ϵ = 1.5), convert a substantial amount of their oxygen into water (Fig. 4). For shallower density gradients, the outer disc contains more mass than the inner disc; here temperatures are low so that water is predominantly frozen on the grain surfaces. What remains is a low abundance of gas phase water above the ice layer. Models that have an inner gap, i.e. R_in = 10 or 100 times the sublimation radius tend to have lower mass fractions of gas phase water. Those models lack the big inner gas phase water reservoir (Woitke et al. 2009b).

3.3. Oxygen gas temperatures

We calculate average gas temperatures for various species from the models following (1)where n_sp is the species density and m_sp its mass. M_sp is the total mass of a species in the disc model. The oxygen gas temperature distribution peaks at 20 K and the median gas temperature is 40 K. The distribution shows that  ~ 60% of the models stay between 15 and 70 K. The temperature distributions for specific subsamples of parameters, such as the the stellar luminosity and disc gas mass, are shown in Figs. 5 and 6. If we split the sample by luminosity, the higher luminosity models show on average higher gas temperatures. This can be directly attributed to the higher energy deposition rate which leads to an increased gas heating. Looking at low mass disc models (Fig. 6), we note that their average gas temperature is almost a factor three higher than that of high mass disc models. This is explained by optical depth effects: in models with moderate optical depth, a higher fraction of the total mass can be heated by irradiation. Hence, models with lower disc mass tend to have on average higher gas temperatures. These general trends are based on a huge number of disc models and the actual distribution of gas temperatures is fairly wide with a long shallow tail towards very high gas temperatures (Fig. 6).

The O i gas temperature statistics (and emerging line fluxes) have already been used in Pinte et al. (2010) to identify the need for an extra heating mechanism from the [O i] line fluxes of T Tauri stars. In that particular case, the stellar luminosity of the low mass stars does not provide enough energy to explain the observed flux levels. Our models show that a UV excess (f_UV = 0.1), which can originate either from a chromosphere or from ongoing accretion, can provide this extra energy and reproduce the observed fluxes. X-ray heating is an alternative process for T Tauri discs (Nomura et al. 2007; Gorti & Hollenbach 2008; Aresu et al. 2011).

Figure 7 shows the impact of dust settling on the coupling of gas and dust temperature for the regions in which atomic oxygen dominates. Since atomic oxygen is widely distributed in the disc, we can take this to first order as an estimate of the amount of thermal gas-dust de-coupling. Well mixed models, s = 0, show a fairly narrow peak around T_g = T_d and a tendency to have gas temperatures higher than dust temperatures (often in the form of a hot surface gas layer). The settled models, s = 0.5, show on average gas temperatures that are lower than those of the dust. This is not due to the presence of cool surface layers in settled models (see Sect. 3.6); the mass averaged O i gas temperature is hardly affected and so the reason lies in the average dust temperature. In settled models, only the small grains stay behind in the disc surface, while the larger ones settle towards the midplane (see Table A.1). Hence, the effective grain size in the surface of settled models is smaller compared to non-settled models, thereby decreasing the dust emissivity at long wavelength. This leads to slightly warmer dust.

3.4. CO gas temperatures

Figures 8 and 9 show the distribution of CO mass averaged gas temperatures in the disc models as a function of the outer radius and the surface density power law index. Given our choice of surface density power laws, it is evident that the CO mass is dominated by the material close to the outer radius. Since the gas temperature is also a clear function of distance, models with larger outer radius show on average lower CO gas temperatures. Figure 8 also illustrates the limitations of the concept of a single canonical CO temperature T_CO to convert submm line fluxes into disc gas masses. We will get back to this point when looking for possibilities to estimate disc gas masses from line fluxes in Sects. 5.2 and 5.3. The surface density power law exponent also impacts the distribution of mass in the disc. Hence, models with a low ϵ put more mass in the outer disc compared to the inner disc; hence, they show on average lower CO gas temperatures. In addition, the trends with luminosity and disc mass found for oxygen, hold also for CO.

For a high percentage of models, the average gas temperatures for atomic oxygen and CO are very similar. Except for the lowest mass discs, both temperatures are within a factor two of each other. This also holds across all gas-to-dust mass ratios as can be seen from Fig. 10.

3.5. Ionized carbon in discs

Ionized carbon forms a thin skin on the disc surface with a typical column density of 10¹⁷⁻¹⁸ cm^-2. Very low mass disc models, M_gas ≤ 10^-5   M_⊙, are so optically thin that carbon is fully ionized throughout the disc resulting in a constant M(C ⁺ )/M_gas fraction of  ~1.2  ×  10^-3. This is only weakly dependent on detailed disc parameters.

As soon as the ionizing UV radiation is blocked by sufficient dust opacity, carbon turns atomic and subsequently molecular. The higher mass models show thus a decreasing fractional C⁺ mass, with the total M(C⁺) levelling off at  ~10^-8   M_⊙. The variation between models can be an order of magnitude mostly depending on the choice of disc extension and stellar luminosity (Fig. 11): the C⁺ mass is higher for larger outer disc radii and higher luminosity stars. The gas residing in the outer disc is generally more optically thin – also to interstellar UV radiation – than the gas in the inner disc regions.

The dependency of the C⁺ gas temperatures on stellar luminosity is less pronounced than for atomic oxygen. The distribution has a strong peak at low temperatures around 20 K except for luminosities L_∗ > 3.6   L_⊙ (Fig. 12). There is a clear dependency on the outer radius similar to CO, with large models (R_out = 300,500 AU) having an average C⁺ gas temperature of  K. The small models (R_out = 100 AU) show a much broader flat distribution between 25 and 100 K. The reason is that the C⁺ mass averaged temperature is dominated by the outermost regions – just as for CO – since it is evenly distributed over the entire disc surface. Hence, the inner radius plays no significant role here.

3.6. Cool surface layers

A subset of DENT grid models shows very cool surface layers on top of warmer midplanes, i.e. models with a positive gas temperature gradient (as seen from the surface). These are very flat, massive discs in which the grains have started to settle. Typical models in this category show M_gas ≥ 0.01   M_⊙, no UV excess, settled dust (s = 0.5) and are non-flaring (β = 0.8,1.0). These models appear here for the first time because of the wide unbiased parameter study.

A related effect has been seen earlier in models by Aikawa & Nomura (2006) who studied grain growth in disc models that do solve iteratively for the vertical hydrostatic equilibrium. For well mixed models with a large maximum grain size a_max, those authors find a similar effect. Photoelectric heating is diminished in those models resulting in lower gas temperatures in the disc atmosphere. On the other hand, optical depth is reduced as well allowing stellar UV photons to penetrate deeper and heat the gas and dust at intermediate heights. The authors also find that discs with larger a_max flare less than those with small molecular cloud like grain size distributions. This is an immediate result of the lower gas temperatures in the surface layers. For their largest maximum grain size a_max = 10 cm, they find that the gas in the disc atmosphere becomes cooler than the dust. In a subsequent paper Nomura et al. (2007) included the effect of grain settling and X-ray heating and show that – compared to the standard well-mixed model with molecular cloud like grains – this indeed lowers the gas temperature in the surface and increases the gas and dust temperatures close to the midplane. The gas temperature gradient becomes much flatter, but there is no pronounced positive temperature gradient. Unfortunately we cannot disentangle from their published model the effect of additional gas heating by X-rays and decrease in gas photoelectric heating due to dust growth and settling.

Jonkheid et al. (2007) have studied a series of four Herbig discs with an increasing gas-to-dust mass ratio due to grain growth and settling. The basic dust disc models are taken from Dullemond & Dominik (2004), who use the dust midplane temperature to set the scale height of the disc. Jonkheid et al. (2007) do not iterate the vertical structure based on the gas temperature and so their approach is somewhat similar to our choice of parametrized disc models. They find that a large amount of the disc gas mass is indeed cooler than the dust and that the highest gas-to-dust mass ratios lead indeed to a positive gradient in the vertical gas temperature.

In nature, we do not expect such flat, massive discs, with strong settling to be very common. However, the presence of such cool surface layers would have implications for the disc surface chemical composition and line emission. At cold temperatures, neutral-neutral chemistry becomes inefficient and the formation of warm water in the disc surface could be affected. Self-absorption in the dominant cooling line [O i] 63 μm (see Sect. 4.3) and temperature gradients as measured from various ¹²CO and ¹³CO lines could reveal such cool surfaces.

4. Line diagnostics

In the following sections, we explore the diagnostic power of individual emission lines from the DENT grid. We check especially for clean dependencies on disc input parameters such as outer radius, disc mass, dust-to-gas mass ratio, surface density power law exponent etc.

4.1. The C⁺ fine structure line at 158 μm

Since ionized carbon originates from a thin skin around the discs, its mass (column densities) tends to be constant until the disc gas mass itself becomes so low that the disc is transparent to UV photons. The total [C ii] 158 μm flux correlates well with the mass in C⁺ indicating that the line is mostly optically thin.

Figure 13 illustrates this effect nicely through the statistical distribution of [C ii] 158 μm fluxes as a function of disc mass. Down to disc masses of 10^-4   M_⊙, the discs stay optically thick and the C⁺ mass fraction and hence its emission stays constant around  ~3  ×  10^-18 W/m². The large width of the flux distribution for these models reflects the spread in stellar luminosity, UV excess, outer disc radius, flaring angle; the higher fluxes results from models with high luminosities (stellar and/or UV – see Fig. 14), large outer radii, and flaring geometry. The [C ii] line flux depends only weakly on the disc surface density gradient, the gas-to-dust mass ratio, settling, the minimum grain radius a_min, and the inner disc radius. For lower gas masses, the peak of the [C ii] distribution shifts towards lower fluxes and the distribution becomes much narrower.

The DENT grid was generated using only two values of UV excess, 10% and 0.1% of the total stellar luminosity. While these values bracket the UV excess observed for young T Tauri stars, it can lead to an overestimate of the UV radiation for Herbig Ae/Be stars that do not possess an excess in the first place. Figure 14 shows that the [C ii] emission increases by an order of magnitude for the high UV models. The Herbig Ae/Be models with pure photospheric UV fluxes would fall even below the red points in Fig. 14.

In practice, the [C ii] 158 μm line is easily contaminated by the low density surrounding gas, even with Herschel/PACS. Although the line is potentially a good tracer as shown above, it requires a very careful observational procedure and data reduction to ensure that the flux measured is indeed purely from the disc.

4.2. Oxygen fine structure lines at 63 and 145 μm

In the DENT grid, 80% of the disc models show a 145   μm line that is a factor 10 − 100 weaker than the 63   μm line. Given that the Herschel sensitivity limits for both lines are very similar, this clearly favors the detection of the 63   μm line.

Woitke et al. (2010) have shown that the 63   μm fine structure line could be promising as an order of magnitude disc gas mass estimator. Their Fig. 4 illustrates that if we can estimate the average excitation temperature of neutral oxygen, the measurement of the 63   μm line can be used to derive gas masses in the context of a statistical analysis of a large sample of discs such as the GASPS sample.

4.3. The [O i] 63/145 line ratio

In PDRs, low [O i] line ratios ( < 10) are expected if both lines are optically thick (Tielens & Hollenbach 1985) and the gas is cooler than  ~200 K. Liseau et al. (2006) argue that confusion with cool envelope gas might explain ratios of a few in observations of YSO’s. Lorenzetti et al. (2002) suggested that high optical depths in the [O i] 63 and 145 μm lines can explain the [O i] 63/145 ratios of 2 − 10 observed by ISO/LWS in massive discs around early-type HerbigAeBe stars.

The grid results presented here confirm that the [O i] 63/145 ratios for massive (>10^-2   M_⊙) discs is of the order of 1 − 10. In our models, we observe that [O i] 63/145 < 10 can be produced by the discs themselves, when [O i] 63 μm is about to go into absorption, whereas [O i] 145 μm is still in emission. As a fundamental line, [O i] 63 μm can go into central absorption if the emitted light has to pass through another cold neutral oxygen layer above, whereas the [O i] 145 μm line practically never goes into absorption. These effects happen for massive discs only, in particular for settled dust distributions, and are strongly dependent on inclination, see Fig. 15 and also Sect. 3.6. The hypothesis that such models exist is a unique testcase for the SOFIA/GREAT instrument that can observe the [O i] 63 μm line with high spectral resolution and thus check for self-absorption within the line profile.

The median [O i] line ratio for the entire DENT grid is 25. Contrary to PDRs, the [O i] 63/145 line ratio of disc models is not sensitive to the average oxygen gas temperature over a wide range of temperatures between 50 and 500 K. The line ratio correlates instead with e.g. the dust-to-gas mass ratio (Fig. 16).

4.4. The CO low J rotational lines

The low J¹²CO rotational lines are optically thick and arise mainly from the outer disc surfaces. Hence they form an excellent tracer of the outer disc radius. The larger the disc, the larger the emitting surface. The oxygen fine structure line at 63 μm is not affected by the disc outer radius for R_out ≥ 100 AU (Kamp et al. 2010). Figure 17 shows the mean [O i] line flux plotted against the CO 2 − 1 flux with the outer disc radius coded in colour. In this diagram, the disc models clearly separate into small ones with R_out = 100 AU and ones that are larger than 300 AU. This indicates that a combination of the oxygen and CO line can help to remove the outer disc radius degeneracy in measuring disc masses.

Detection limits with current radio telescopes such as APEX or SMA are a few times 10^-20 W/m². With ALMA coming along soon, we will be able to detect much fainter discs down to  ~10^-22 W/m². Figure 18 shows the CO 2 − 1 line flux versus the dust continuum at 1.2 mm with the respective ALMA band 6 detection limits. These grid results suggest that ALMA is sensitive enough to detect disc masses of the order of 10^-6   M_⊙.

4.5. The [O i] 63/CO 2–1 line ratio

Figure 19 illustrates that the [O i] 63/CO 2 − 1 line ratio correlates with the average oxygen gas temperature in the range 15 − 70 K. If the [O i] 63 μm line is optically thin, it depends on disc gas mass and average oxygen gas temperature (Woitke et al. 2010). The CO 2 − 1 line is always optically thick (down to M_gas ~ 10^-4   M_⊙) and is therefore a function of the CO gas temperature and disc surface area, i.e. R_out. Using the two lines [O i] 63 μm and CO 2 − 1, we also get a handle of the amount of extra UV radiation irradiating the disc (see Fig. 20).

From Figs. 17, 19 and 20, it becomes clear that this [O i]/CO line ratio captures in a complex way a number of disc properties such as the average gas temperature, the disc outer radius, the flaring index, and the amount of extra UV irradiation (f_UV). Thus this line ratio is a powerful means to break the modelling degeneracy and find a suitable method to estimate gas masses even in the case of very few observational constraints. This is especially relevant for large surveys such as GASPS (Dent et al. 2011, in prep.) that contains many objects where only the [O i] 63 μm fine structure line is detected and ancillary data is mostly restricted to SED’s and CO submm lines.

4.6. Water lines

Due to the molecular structure, water emission lines are not only difficult to calculate (non-LTE, population inversion, optical depth effects), but also difficult to interpret in terms of global disc parameters such as gas mass or surface density profiles (Cernicharo et al. 2009). Woitke et al. (2009b) traced the origin of water lines in the discs around Herbig Ae stars back to three main water reservoirs: (1) a main water reservoir in the central midplane out to  ~10 AU (containing most of the water mass); (2) a cold water belt at r ≈ 10−100 AU and z/r ≈ 0.05 (originating from desorption of water ice mantles); and (3) a hot water surface layer at distances of 2−25 AU (originating from an active neutral-neutral chemistry at T > 200 K). Water lines that originate deeper in the disc suffer from the problem that gas and dust temperature couple at some fiducial depth, making it potentially impossible to see the line emission above the strong dust continuum (Woitke et al. 2009b; Pontoppidan et al. 2010).

Figure 21 shows that for massive discs (M_gas ≥ 0.01   M_⊙), the water lines are generally stronger in the settled models. The effect of dust settling on the water chemistry and line emission is very complex as it affects the available dust surface area, the temperature gradients and the continuum optical depth. Disentangling those effects simply from the grid observables is impossible. Hence, we discuss here two exemplary models to illustrate the general trend displayed in Fig. 21. Picking a settled and non-settled flaring disc model with a_min = 1.0   μm (M _∗  = 2   M_⊙, T_eff = 4960 K, 3.7   L_⊙, R_in,out = 17,500 AU, M_gas = 0.1   M_⊙, δ = 0.01, ϵ = 0.5, β = 1.0), we find that the dust in the surface layers of the settled model is warmer on average, while the gas temperature barely changes (using Eq. (1) for oxygen, it changes from 13 to 16 K). This was already discussed in Sect. 3.3.

Even though the temperatures are low, there is very little grain surface area available for e.g. H₂ or ice formation. The H₂ formation on grain surfaces is thus inefficient and atomic hydrogen abundances are high in these settled models. As a consequence of the lack of grain surfaces in the settled model, molecular hydrogen forms via the H⁻ route and the atomic hydrogen abundance stays a factor 100 higher than in the non-settled models. This leads to an efficient cold water formation route driven by radiative association The relevance of radiative association for cool atomic H regions was already pointed out by Glassgold et al. (2009), but found to be irrelevant for the inner disc structures that they studied. These reactions have low rate coefficients¹, but with atomic hydrogen being a factor 100 more abundant in settled models, they dominate the formation route of water. The surface layers (z/r > 0.2) of the settled model show water abundances that are several orders of magnitude higher than those of the unsettled model.

5. Diagnostic tools

In the following, we discuss a few methods to measure the disc gas masses from gas line observations. It is evident by the way how we set up the parameter space that the DENT grid does not reflect a proper distribution of disc properties, e.g. number of T Tauri versus Herbig discs, distribution of disc masses reflecting disc lifetimes etc. Hence, our aim is not to decide which method is better, but to put the various methods into the context of disc modelling and discuss some of the limitations and uncertainties that become evident if one confronts the simple methods with a wide variety of disc appearances (sizes, flaring indices, gas-to-dust mass ratios, etc.).

The most common technique to estimate disc gas masses has been the use of the mm dust continuum flux. At those wavelengths, the disc is assumed to be optically thin. Hence, the mm-flux should reflect the entire dust mass contained in large grains. The main uncertainties in converting the mm-flux into a total dust and subsequently gas mass are the dust opacities κ_ν ~ κ₀(ν/ν₀) ^− β and the dust-to-gas mass ratio δ. An average opacity index β ~ 0.46 has been recently found for the discs in the ρ Ophiuchus star forming region (Ricci et al. 2010); these low values of β (compared to 1.7 for the dust in the interstellar medium) are generally interpreted in terms of grain growth, see also Draine (2006). However, recent interferometric data suggests that β varies with location in the disc (e.g. Isella et al. 2010). The dust-to-gas mass ratio is generally assumed to be the canonical value of δ = 0.01 found for the interstellar medium neglecting potential differences in dust and gas evolutionary timescales. Even if the dust mass could be estimated very precisely, this method would still leave us with large uncertainties on M_gas because of the unknown conversion factor. Hence, the intrinsic uncertainties of this method are much greater than an order of magnitude.

Very few direct estimates of the gas-to-dust mass ratio exist for discs. An example is the study by (Glauser et al. 2008) that estimates dust and gas mass on the same disc line-of-sight in a disc. The derived ratio of  is compatible with the canonical value of 100.

The following three sections discuss potential methods only based on gas emission line fluxes, their limitations, and uncertainties.

5.1. [O i] and [C ii] fine structure line ratios

Figure 22 shows a two-colour line flux diagram suggesting a method how to determine the disc gas mass purely from three line flux observations. As was argued by Kamp et al. (2010), the [O i] 63/145 line ratio shows a clear correlation with gas mass. This plot includes the results for all types of central stars, low and high excess UV, flared and non-flared discs, continuous discs as well as discs with big inner holes, and varying dust parameters and outer disc radii. Despite this variety of disc types, all models with a certain disc gas mass fall into certain regions in this two-colour diagram. Apparently, the main dependencies on f_UV and β tend to cancel out here, because these parameters shift all line fluxes up and down equally.

However, the diagram looks different for other sub-selections of models as the one stated in the caption of Fig. 22. A steeper surface density power law (higher ϵ) weakens in particular the [C ii] 158 μm flux and hence shifts the points rightward in the figure.

Selecting only settled models has shifted the location of massive models (M_gas ≥ 0.01   M_⊙) downward in Fig. 22 compared to a similar plot for unsettled models. The reason for this is mainly [O i] 63 μm self-absorption in the cool surface layers as discussed in Sect. 3.6.

In addition, there are other parameters such as the outer disc radius R_out and the dust size parameters, which introduce additional non-trivial dependencies and also a large scatter in the various mass bins. Thus a mass determination exclusively based on these three line ratios would implicate large errors especially for the higher mass disc models.

5.2. Low J CO lines

For many objects, the [O i] 145 and [C ii] 158 μm lines will not be available as they are often at least a factor 10 fainter than the [O i] 63 μm line. Therefore, we study in the following also the CO low J lines that are accessible from the ground. Especially with ALMA coming up, these lines will be detectable for a large number of discs (see Sect. 4.4).

The low J¹²CO lines are often optically thick and can hence provide at most a lower limit to the total gas mass. Dent et al. (2005) present a formula to convert the CO 3−2 line flux I_{CO   3−2} into a gas mass assuming LTE and a constant ¹²CO/H₂ conversion factor of 5  ×  10^-5. We adopt here their approach and derive a relation between gas mass and intensity (4)where I_{CO   3 − 2} is the integrated CO 3 − 2 line intensity in K km s^-1. Since most of the CO mass will reside in the outer disc, the line flux will be dominated by cold gas. Hence, an average gas excitation temperature of T_ex = 50 K is assumed in the following estimates. To apply Eq. (4) to our grid data, we need to convert the modeled integrated line fluxes in W/m² to K km s^-1 (formula adapted from Thi et al. 2004) (5)Here, λ = 8.6696  ×  10^-2 cm and k is the Boltzman constant in erg/K. The solid angle is assumed to be Ω = π⁽HPBW/2⁾², with HPBW = 13.7′′ the half power beamwidth of the JCMT at 345 GHz. The formulas for the CO 2−1 line look similar; however, the line is generally a factor 2.6 weaker than the CO 3−2 line reflecting the difference in the Planck function (ν₃₋₂/ν_2 − 1)³.

We insert our modeled fluxes first into Eq. (5) and then the converted intensities into Eq. (4) to derive a gas mass estimate. These derived values underestimate the actual gas mass, especially at the high mass end, where the CO 3−2 line is optically thick. Figure 23 shows that the behaviour is much more complex and that there is a regime, where the simple formula even overestimates the disc gas mass, most likely due to the simplified assumption of a uniform excitation temperature of 50 K for all models. Most notably, the standard deviation across all disc parameters for this method to measure disc gas mass is large, generally more than two orders of magnitude.

5.3. The [O i] 63/CO 2–1 line ratio

We proceed thus to suggest here an alternative approach using the insight on the diagnostic power of the [O i] 63/CO 2−1 line ratio from the previous sections. Figure 24 shows for each bin of [O i] 63/CO 2−1 line ratios the resulting fit using a 3rd order polynomial. The coefficients for these polynomial fits are presented in Table 2. In the following, we check how accurate the disc masses derived from a combination of these two lines are.

For every model, we calculate the [O i] 63/CO 2−1 line ratio and then select the corresponding curve from Fig. 24. This curve in combination with the calculated [O i] 63 μm flux results in a “measured” disc gas mass. This value is then compared to the actual gas mass in the models (Fig. 25). The standard deviation indicates that this method can provide an order of magnitude estimate of the disc gas mass. The method does frequently overestimate the disc gas mass by a factor two. Figure 25 also shows that the relation levels off at disc masses beyond 10^-3   M_⊙ where the [O i] 63 μm line becomes optically thick.

For a large statistical sample in which individual disc parameters are less well known, this method can indeed help to assess gas evolution in discs and to measure the gas dispersal timescale. The errors of this method are much smaller than the ones from the CO 3−2 line alone. Of course, more observational constraints on any object should naturally lead to more accurate measurements.

6. Conclusions

We have calculated a grid of 300 000 disc models using the dust Monte Carlo radiative transfer code MCFOST and the gas chemistry and heating/cooling balance code ProDiMo. These models span an 11-dimensional parameter space and enable us to study disc properties such as temperatures, species masses, gas line fluxes, and dust continuum fluxes in a statistical manner. The goal is to find suitable gas diagnostics that allow us to characterize the main disc parameters such as radial extent, gas mass, and flaring from observational quantities.

The main results found from an extensive statistical analysis of the entire DENT model grid are:

• The[O i] 63 μm traces disc mass but for massive non-flaring, settled disc models without UV excess, the line is optically thick and often suffers from self-absorption in upper cool disc layers.

• The [C ii] 158 μm fine structure line flux is very sensitive to the stellar UV flux and presence of a UV excess.

• The CO low rotational lines trace the outer disc radius.

• Dust settling enhances the water abundance in the surface layers through an efficient radiative association channel. The initial step of H₂ formation proceeds in the gas phase via H⁻.

• Low [O i] 63/145 line ratios (<a few) can be explained with cool gas in the uppermost surface layers instead of 63   μm absorption by cool envelope remnant gas as suggested by Liseau et al. (2006).

• The [O i] 63/CO 2−1 line ratio correlates with several disc properties such as the average O i gas temperature in discs, the outer disc radius, and the UV excess. As such it is a powerful diagnostic to break disc modelling degeneracies.

• A combination of the [O i] 63 μm flux and the [O i] 63/CO 2−1 line ratio can be used for M_gas ≤ 10^-3   M_⊙ to obtain an order of magnitude estimate for the disc gas mass purely from gas observations. The previously used conversion of a CO submm line flux alone generally leads to larger uncertainties.

---

Acknowledgments

F.M., C.P., W.F.T., and J.C.A. acknowledge PNPS, CNES and ANR (contract ANR-07-BLAN-0221) for financial support. C.P. acknowledges funding from the European Comission seventh Framework Program (contracts PIEF-GA-2008-220891 and PERG06-GA-2009-256513). I.K. acknowledges funding from an NWO MEERVOUD grant.

References

Appendix A: The DENT grid

The general outline of our numerical pipeline is described in Woitke et al. (2010) (see their Fig. 1). In the first step, MCFOST is used to compute the dust temperature T_dust(r,z) and frequency dependent radiation field J_ν(r,z) from radiative transfer. This information is passed on in the second step to ProDiMo to compute the detailed gas chemistry n_mol(r,z) and temperature T_gas(r,z). The resulting gas structure together with the non-LTE level population numbers of the relevant species is inserted back into MCFOST to do the line radiative transfer and calculate the line profiles and integrated fluxes (step 3). In a second radiative transfer step, dust spectral energy distributions (SED) are calculated for all models (step 4).

A.1. Parametrized disc models

For the shape and mass distribtion of the gas in the disc, we use the following parametric description (A.1)between an inner and outer disc radius, R_in and R_out, respectively, with sharp edges. ρ(r,z) is the local gas mass density. The constant ρ₀ is adjusted such that the integrated disc mass 2π ^! ρ(x,z)   dz   r   dr equals M_disc. H(r) is the vertical scale height of the disc, assuming to vary with radius as (A.2)H₀ is the reference scale height at reference radius r₀. ϵ is the column density powerlaw index and β the flaring power.

The dust grains are assumed have a unique powerlaw size distribution (A.3)between minimum grain radius a_min and maximum grain radius a_max. The free constant in Eq. (A.3) is adjusted to result in the specified dust/gas mass ratio ρ_d/ρ (see Sect. 4.6 in Woitke et al. 2009a), which is assumed to be constant throughout the model volume (e.g. no dust settling). The dust absorption and scattering opacities are calculated by applying Mie theory with optical constants for astronomical silicate from Draine & Lee (1984). The DENT grid comprises also models with dust setlling, although that particular parameter is only sampled by the two extremes of non-settled models and models in which grains larger than 0.05   μm have settled to a smaller scale height than the gas, H(r,a) ≈ H(r)(a/0.05   μm)^-0.25 (for s = 0.5).

Table A.1 summarizes the free parameters in the DENT grid. For completeness, it also contains the fixed parameters r₀, H₀, χ_ISM, ζ_CR, f_PAH, α, a_max, ρ_gr, and v_turb.

Effective temperature and luminosity are related parameters and are consistently chosen from the evolutionary tracks of pre-main sequence stars (Siess et al. 2000) for solar metallicities. Table A.2 lists the parameter combinations chosen for the representative five stellar masses 0.5, 1.0, 1.5, 2.0, and 2.5 M_⊙. Disc models are calculated at four stellar evolutionary times, 1, 3, 10, and 100 Myr. The stellar input spectra are interpolated from Phoenix stellar atmosphere models (Brott & Hauschildt 2005) with solar metallicities.

At UV wavelengths, however, where stellar activity and mass accretion create an UV excess with respect to the classical stellar atmosphere models, we switch to a powerlaw UV input spectrum with spectral intensity [erg/cm²/s/Å/sr] (A.4)The flux is scaled to yield the prescribed fractional UV luminosity of f_UV = L_UV/L _⋆  = 0.1 (high UV) or 0.001 (low UV) with the UV luminosity L_UV being integrated from 912 Å to 2500 Å as introduced by Woitke et al. (2010). Very few stars have a well known f_UV (given the wavelength range defined above and the gaps in observational data of e.g. STIS and FUSE). One of the few stars is ET Cha, which has a value of 0.02 (Woitke et al. 2011).

A.2. The observables

The quantities that can be compared to observations such as the gas line emission and dust SED are computed in step 3 and 4 of the grid modelling procedure (see introduction to Sect. 2). All quantities are computed for 5 different inclinations that are uniformly distributed, 0^o, 41.41^o, 60^o, 75.52^o, and 90^o (edge-on).

The dust SED is computed at 57 wavelength points between 0.1 and 3500 μm. The line profiles and integrated fluxes of four species are computed: C ii, O i, CO, and H₂O (ortho and para). Table A.3 provides an overview of the respective atomic/molecular line data.

A.3. Disc stability

We choose here the Toomre criterium to check our disc models against gravitational instabilities. The Toomre Q parameter is defined as (A.5)where c_s is the sound speed, Σ the surface density, G the gravitational constant and the epicyclic frequency κ is defined as (A.6)with Ω the orbital frequency. The disc is unstable at a certain distance r from the star if Q(r) < 1. In the following, we derive an expression that provides a simple analytical way of checking “global” disc instability. For that, we make a few simplifying assumptions

• 1.

  The discs rotate with keplerian frequency, hence κ(r) = Ω_kep(r).

• 2.

  The midplane gas temperature can be approximated by a radial power law T = T₀(r/R_in)^-0.5.

We will discuss the impact of these simplifications further below.

Inserting the power laws for temperature and surface density into Eq. (A.5), we obtain (A.7)Here T₀ and Σ₀ denote the gas temperature and surface density at the inner radius of the disc R_in. Since for our values of ϵ (0.5, 1.0, 1.5) Q(r) decreases with radius, disc models that obey Q > 1 at the outer radius are also stable at smaller radii. Hence, in the following, we derive an expression for Q at r = R_out, which can be used to determine the global stability of a disc model.

The surface density and temperature at the inner radius can be derived from the total disc gas mass M_disc and the luminosity and effective temperature of the central star L _∗ , T_eff (Tuthill et al. 2001; van der Plas et al. 2011) Using these, we can rewrite Eq. (A.7) entirely in terms of the input parameters of the DENT grid (A.10)We apply this criterium to all 322 030 disc models and find that only 16 422 of them are unstable according to the above derived criterium. Hence, 94% of all models are not affected by gravitational instabilities. This percentage changes by less than 1% if we (a) choose a more realistic temperature profile for optically thick discs T ~ (r/R_in)^-0.25 or (b) change T₀ by a factor 2. The fact that Q(r) decreases with r still holds for T ~ (r/R_in)^-0.25. Overall, we find that models with high gas disc masses – as expected – are more subject to instabilities. However, even in our highest mass bin the percentage of unstable models is only  ~ 30%. The plots and conclusions of the analysis of the DENT grid are not affected by the small number of models that could become unstable.

A.4. Methodology

The high number of disc models makes it impossible to carry out an individual study of their chemical and thermal structure. Hence in the following, we study the statistics of species gas temperatures and masses and how they depend on certain grid parameters such as the surface density distribution and settling.

For the statistical analysis, we define a scalar quantity X_i of any model i. It can be a parameter, or can be a measured value like the total mass of atomic oxygen M(O) of model i. Next, we select an interval on which we perform the statistics  [X_k − dx/2,X_k + dx/2], where X_k are sampling points for X. The total number of models (selected i) in that selected range are N_sel. The expected value and standard deviation of the distribution are then defined as

All Tables

All Figures