© ESO, 2016

1. Introduction

Disk models are widely used by the community to analyse and interpret line and continuum observations from protoplanetary disks, such as photometric fluxes, low- and high-resolution spectroscopy, images and visibility data, from X-ray to centimetre wavelengths. Historically, disk models could be divided intocontinuum radiative transfer models, such as HOCHUNK3D (Whitney et al. 2003), MC3D (Wolf 2003), RADMC (Dullemond & Dominik 2004a), TORUS (Harries et al. 2004), MCFOST (Pinte et al. 2006) and MCMax (Min et al. 2009), to explore the disk shape, dust temperature and grain properties, and thermo-chemical models, see e.g. (Henning & Semenov 2013) for a review, and Table 1. The thermo-chemical models usually include chemistry and UV and X-ray physics to explore the temperature and chemical properties of the gas, with particular emphasis on the outer disk as traced by (sub-)mm line observations.

However, this distinction is becoming more and more obsolete, because new disk models try to combine all modelling components and techniques, either by developing single, stand-alone modelling tools likeProDiMo (Woitke et al. 2009) and DALI (Bruderer et al. 2014), or by coupling separate continuum and gas codes to achieve a similar level of consistency.

Observational data from protoplanetary disks obtained with a single observational technique in a limited wavelength interval can only reveal certain information about the physical properties at particular radii and at particular vertical depths in the disk. Therefore, in order to derive an overall picture of protoplanetary disks, it is essential to combine all observational data, and to make consistent predictions for all continuum and line observables in a large range of wavelengths on the basis of a single disk model.

However, this holistic modelling approach does not come without a price. The number of free parameters in such models is large (around 20), and the computational time required to run one model can exceed days, weeks or even months (e.g. Semenov et al. 2006). These limitations have resulted in quite limited parameter space being explored in such models.

Therefore, previous chemical models have not fully explored the role of disk shape and dust opacities. In Table 1, we list assumptions made in differnt thermo-chemical disk models about the shape of the disk, the dust size distribution, opacities, dust settling, and polycyclic aromatic hydrocarbons (PAHs). The selection of models is not exhaustive in Table 1, for a more comprehensive overview of modelling techniques and assumptions see Table 3 in (Henning & Semenov 2013). Table 1 shows the diversity of modelling assumptions currently used by different disk modelling groups. These models often focus on the outer disk, consider small dust particles, and use different approaches for dust settling and PAHs.

All these assumptions have crucial impacts on the modelling results, not only with regard to the predicted continuum observations, as known from SED fitting, but also on chemical composition and line predictions. Therefore, to compare modelling results from a large number of protoplanetary disks, a set of consistent standard modelling assumptions is required.

In this paper, we explore what could be a minimum set of physical assumptions about the star, the disk geometry, the dust and PAH opacities, dust settling, gas and ice chemistry, gas heating and cooling, and line transfer, needed to capture the most commonly observed multi-wavelength properties of Class II and III protoplanetary disks. We aim at a consistent, coupled modelling of dust and gas, and we want to predict the full suite of observations simultaneously and consistently from one model.

In Sect.2 we introduce our aims and basic approach, in Sect. 3 we describe the details of our model, and in Sect. 4 we cross-check and verify the implementation of our assumptions into our three main modelling tools ProDiMo, MCFOST, and MCMax. In Sect. 5, we show first results for a simple model of a Class II T Tauri star, and systematically study the impact of all modelling parameters on the various predicted observables. We summarise these results and conclude in Sect. 6.

In the Appendices, we collate a number of auxiliary information. We explain our fitting routine of stellar parameters including UV and X-ray properties (Appendix A), and describe our assumptions concerning interstellar UV and IR background radiation fields. We compare some results obtained with the simplified treatment of PAHs (see Sect. 3.8) against models using the full stochastic quantum heating method in Appendix B. Appendix C compares the results obtained with time-dependent chemistry against those obtained in kinetic chemical equilibrium. We detail how certain observable key properties are computed from the models in Appendix D. Appendix E discusses the behaviour of optically thick emission lines. In Appendix F, we discuss the convergence of our results as function of the model’s spatial grid resolution.

Two forthcoming papers will continue this paper series, to study the impacts of chemical rate networks (Kamp et al., in prep., Paper II) and element abundances (Rab et al. 2015, Paper III) on the resulting chemical abundances and emission lines. Kamp et al. introduce a simplified, small chemical rate network and a more exhaustive, large chemical network, henceforth called the small DIANA chemical standard and the large DIANA chemical standard, respectively.

2. The DIANA project

The European FP7 project DiscAnalysis (or DIANA) was initiated to bring together different aspects of dust and gas modelling in disks, alongside multiwavelength datasets, in order to arrive at a common set of agreed physical assumptions that can be implemented in all modelling software, which is a precondition to cross-correlate modelling results for different objects.

The DIANA goal is to combine dust continuum and gas emission line diagnostics to infer the physical and chemical structure of Class II and III protoplanetary disks around M-type to A-type stars from observations, including dust, gas and ice properties. The project aims at a uniform modelling of a statistically relevant sample of individual disks with coherent observational data sets, from X-rays to centimetre wavelengths.

In protoplanetary disks, various physical and chemical mechanisms are coupled with each other in complicated, at least two-dimensional ways. Turbulence, disk flaring, dust settling, the shape of the inner rim, UV and X-ray irradiation, etc., lead to an intricate interplay between gas and dust physics. Therefore, consistent gas and dust models are required. Each of the processes listed above can be expected to leave specific fingerprints in form of observable continuum and line emissions that can possibly be used for their identification and diagnostics.

However, to include all relevant physical and chemical effects in a single disk model is a challenging task, and we have decided to include only processes which we think are the most important ones, clearly with some limitations. We also want to avoid approaches that are too complicated and pure theoretical concepts that are not yet verified by observations. The level of complexity in the models should be limited by the amount and quality of observational data we have to check the results. The result of these efforts are our disk modelling standards that we are proposing to the community in this paper series.

New challenges for disk modelling have emerged with the advances in high-resolution imaging (e.g. ALMA, NACO, SPHERE and GPI). These observations show evidence for non-axisymmetric structures such as spiral waves, warps, non-aligned inner and outer disks, and horseshoe-like shapes in the sub-mm. In the future, 3D models are clearly required to model such structures, but these challenges go beyond the scope of this paper. This paper aims at setting new 2D disk modelling standards as foundation for the DIANA project, sufficient to reproduce the majority of the observations, simple to implement, yet physically established, and sufficiently motivated by observations. We will offer our modelling tools and collected data sets to the community¹.

3. Standard disk modelling approach

3.1. Stellar parameters and irradiation

To model Class II protoplanetary disks, we need to specify the stellar and interstellar irradiation at all wavelengths. This requires determining the photospheric parameters of the central star, i.e. the stellar luminosity L_⋆, the effective stellar temperature T_⋆ and the stellar mass M_⋆. From these properties, the stellar radius R_⋆ and the stellar surface gravity log (g) can be derived. Fitting these stellar properties to photometric observations is essential for modelling individual disks, which requires knowing the distance d and determining the interstellar extinction A_V. More details about our procedure for fitting the photospheric stellar properties are explained in Appendix A.

Young stars are known to be strong UV and X-ray emitters. This additional, non-photospheric, high-energy disk irradiation can be neglected for pure dust continuum models, but is essential for the modelling of the chemistry and energy balance of the gas in protoplanetary disks. In Appendices A.2 and A.3, we explain how observed UV spectra and measured X-ray data can be used to prescribe this additional high-energy irradiation in detail. If such detailed data is not available, we propose a 4-parameter prescription, the relative UV luminosity f_UV = L_UV/L_⋆, a UV powerlaw index p_UV, the X-ray luminosity L_X and the X-ray emission temperature T_X, see Appendices A.2 and A.3 for details. All stellar irradiation components (photosphere, UV, X-rays) are treated by one point source in the centre of the disk.

In addition to the stellar irradiation, the disk is also exposed to interstellar irradiation: the interstellar UV field, infrared background radiation, and the 2.7 K cosmic background. All these types of irradiation are treated by an additional isotropic irradiation, approaching the disk from all sides (see Appendix A.4).

3.2. Disk mass and column density structure

The gas column density structure Σ(r) [g/cm²] is assumed to be given by a radial powerlaw with index ϵ, modified by an exponential tapering off factor $\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\propto }{\mathit{r}}^{\hspace{0.17em}\mathrm{-}\mathit{ϵ}}\exp \left(\mathrm{-}{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{tap}}}\right)}^{\mathrm{2}\mathrm{-}\mathit{\gamma }}\right)\mathit{,}$(1)where r is the radius and R_tap the tapering-off radius. This approach can naturally explain the often somewhat larger spectral appearance of protoplanetary disks in (sub-)mm molecular lines as compared to millimetre continuum images. If the disk has a tenous continuation, the lines remain optically thick to quite large radii, where the optically thin continuum signal already vanishes in the background noise. For example, the CO J = 2 → 1 line at 1.3 mm only requires a hydrogen nuclei column density of about N_H = Σ/(1.4 amu) ≈ 10 ²¹ cm^-2 in our models to become optically thick, whereas the 1.3 mm continuum requires N_H ≈ 2 × 10 ²⁴ cm^-2. In our reference model (see Sect. 5.1), these column densities correspond to radii of about 450 AU and 20 AU, respectively.

While a tapering-off column density structure according to Eq. (1), with constant dust/gas ratio, seems sufficient to explain the different apparent sizes of gas and dust in many cases, for example (Isella et al. 2007; Tilling et al. 2012) for HD 163296 and (Panić et al. 2009) for IM Lupi, recent high S/N ALMA observations of TW Hya (Andrews et al. 2012) and HD 163296 (de Gregorio-Monsalvo et al. 2013) suggest that the radial extension of millimeter-sized grains can be significantly smaller, with a sharp outer edge, i.e. a varying dust/gas ratio. However, since the predicted timescales for radial migration are in conflict with current observations (e.g. Birnstiel & Andrews 2014), and not much is known quantitatively about migration yet, we assume a constant dust/gas ratio in this paper, to avoid the introduction of additional free parameters.

The default choice for the tapering-off exponent is γ = ϵ (self-similar solution, Hartmann et al. 1998). However, here the two exponents are kept independent, in order to avoid ϵ to be determined by γ in cases where high-quality sub-mm image data allow for a precise determination of γ. We think that the inner disk structure should rather be constrained by observations originating from the inner regions, for example near-IR excess, IR interferometry, ro-vibrational CO lines, etc.. Radial integration of Eq. (1), from R_in to R_out, results in the total disk mass M_disk, which is used to fix the proportionality constant in Eq. (1). The inner rim is assumed to be sharp and positioned at R_in. The outer radius R_out ≫ R_tap is chosen large enough to ensure that Σ(R_out) is small enough to be neglected, e.g. N_H(R_out) ≈ 10²⁰ cm^-2.

3.3. Vertical gas stratification

We assume a Gaussian vertical gas distribution with parametric gas scale height H_g as function of r as $\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{\propto }\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}_{\mathrm{g}}\mathrm{(}\mathit{r}{\mathrm{)}}^{\hspace{0.17em}\mathrm{2}}}\right) \mathrm{with} {\mathit{H}}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{H}}_{\mathrm{0}}\hspace{0.17em}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\hspace{0.17em}\mathit{\beta }}\mathit{,}$(2)where ρ(r,z) is the gas mass density in cylindrical coordinates, H₀ is the reference gas scale height at radius r₀, and β is the flaring exponent. Vertical integration of the gas density ρ(r,z) results in Σ(r) which is used to fix the proportionality constant in Eq. (2). We have chosen this simple approach to be most flexible with our fits of near-IR excess, far-IR excess, visibilities and gas lines. An alternative approach would be to assume vertical hydrostatic equilibrium, either using the calculated dust or gas temperature as input, but these models have some issues reproducing T Tauri disk observations, see Sect. 5.2.2. Also, such models take about 5× to 100× more computational time to complete, because an iteration between radiative transfer, gas physics, and vertical structure is required. This approach is not appropriate when considering the calculation of a large number of models, as is required when fitting observations.

3.4. Dust size distribution

We assume a powerlaw dust size distribution f₀(a)[cm^-4] as function of particle radius a [cm] as ${\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }{\mathit{a}}^{\mathrm{-}{\mathit{a}}_{\mathrm{pow}}} \mathrm{with} \mathit{a}\mathrm{\in }\mathrm{\left[}{\mathit{a}}_{\min }\mathit{,}{\mathit{a}}_{\max }\mathrm{\right]}\mathit{.}$(3)Equation (3) prescribes the dust size distribution function in the disk “before settling”. Dust settling concentrates the larger grains toward the midplane, and therefore, the local dust size distribution f(a,r,z)[cm^-4] will, in general, deviate from f₀(a). Since dust settling only re-distributes the dust particles vertically in a given column, vertical integration over that column must again result in f₀(a) = ∫f(a,r,z) dz/∫dz. The local dust mass density [g/cm³], before settling, is given by $\mathit{\rho }\mathrm{\times }\mathit{\delta }\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{\rho }}_{\mathrm{d}}{{}^{\mathrm{\int }}\mathit{a}_{\min }}_{{\mathit{a}}_{\max }}^{}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\hspace{0.17em}{\mathit{a}}^{\mathrm{3}}\hspace{0.17em}\mathrm{d}\mathit{a}$, where ρ_d is the dust material density, ρ the gas density and δ the assumed unsettled dust/gas mass ratio. This condition is used to fix the proportionality constant in Eq. (3).

The minimum and maximum dust size, a_min and a_max, are set by the following simple considerations: sub-micron sized particles are directly seen in scattered light images, high above the midplane. They are abundant in the interstellar medium (larger grains up to a few μm seem to already exist in dense cores, see Lefèvre et al. 2014), and disks are primordially made of such dust. Millimetre sized grains do also exist in protoplanetary disks, as indicated by the observed SED slope at mm-wavelengths. Therefore, a powerlaw covering the entire size range seems to be the most simple, straightforward option. We will explore the effects of a_min, a_max and a_pow on the various continuum and line observations in Sect. 5.2.1.

3.5. Dust settling

Dust settling is included according to Dubrulle et al. (1995), assuming an equilibrium between upward turbulent mixing and downward gravitational settling. The result is a size and density-dependent reduction of the dust scale heights H_d(r,a) with respect to the gas scale height H_g(r), $\begin{array}{ccc}& & \\ & & \end{array}$where Ω(r) is the Keplerian orbital frequency, γ₀ ≈ 2 for compressible turbulence, and τ_f(r,a) is the frictional timescale in the Stokes regime. ρ_mid(r) is the midplane gas density, and c_T the midplane sound speed. To avoid iterations involving the midplane temperature as computed by dust radiative transfer, we use c_T(r) = H_g(r) Ω(r) here, where H_g(r) is the gas scale height from Eq. (2). α_settle is the dimensionless viscosity parameter describing the strength of the turbulent mixing. The l.h.s. of Eq. (4) is smoothly limited to a maximum value of one by y² → y²/ (1 + y²) with y = H_d(r,a) /H_g(r). Technically, in every disk column, Eq. (4) is computed for a number of (about 100) dust size bins. Starting from the unsettled dust size distribution, the dust particles in each size bin are re-distributed in z-direction according to f(a,r,z) ∝ exp(−z²/ [2H_d(r,a) ²]), building up a numerical representation of the local dust size distribution function at every point in the disk f(a,r,z).

We consider dust settling as a robust physical effect that should occur rapidly in any disk (Dullemond & Dominik 2004b), with the dust grains relaxing quickly toward a vertical equilibrium distribution as described by Eq. (4). Therefore, we think this important effect should be included in radiative transfer as well as in thermo-chemical disk models. Equations (4) and (5) offer an easy-to-implement, yet physically well-justified method to do so, with just a single parameter α_settle.

3.6. Radial zones, holes, and gaps

There is increasing evidence that protoplanetary disks are frequently sculptured by the planets forming in them, which results in the formation of various shape defects, in particular disk gaps (e.g. Forrest et al. 2004; Andrews et al. 2011; Kraus et al. 2013). The gaps are apparently mostly devoid of dust, but may still contain gas as traced by CO rotational and ro-vibrational emission lines (e.g. Bruderer 2013; Carmona et al. 2014). Such objects are classified as transitional disks, with a strong deficiency of near-IR to mid-IR flux. Understanding the SEDs of transitional disks requires to position the inner wall of the (outer) disk at much larger radii than expected from the dust sublimation temperature or the co-rotation radius (e.g. Espaillat et al. 2014). The physical mechanisms responsible for gap formation and disk truncation are still debated, for example planet formation and migration (Lin & Papaloizou 1986; Trilling et al. 1998; Nelson et al. 2000; Zhang et al. 2014), and/or photo-evaporation winds (Font et al. 2004; Alexander et al. 2006; Gorti & Hollenbach 2009), but one observational fact seems to have emerged: disk shape defects are common. In fact, single-zone, continuous protoplanetary disks (e.g. FT Tau, see Garufi et al. 2014) could be a rather rare class (e.g. Maaskant et al. 2014), and for the modelling of individual protoplanetary disks we need an additional option. For an archetypal shape defect, we consider two distinct radial disk zones, with a gap in-between. In such a case, all disk shape, dust and settling parameters come in two sets, one for the inner zone, and one for the outer zone.

3.7. Standard dust opacities

As we will show in this paper, the assumptions about the dust opacities have a crucial impact not only on the predicted continuum observations, but also for chemistry and emission lines. Therefore, the authors of this paper have agreed on a new common approach, which includes a number of robust facts and requirements that are essential to model disks. We will explain these new standard dust opacities carefully in the following, because we think that the dust in disks is different from the dust in the interstellar medium and standard dust opacities so far only exist for the interstellar medium (e.g. Draine & Lee 1984; Laor & Draine 1993).

Our assumptions for the dust opacity treatment are guided by a study of multi-wavelength optical properties of dust aggregates (Min et al. 2016) where the Discrete Dipole Approximation (DDA) is used to compute the interaction of light with complexly shaped, inhomogeneous aggregate particles. These opacity calculations are computationally too expensive to be applied in complex disk models, but Minet al. have developed a simplified, fast numerical treatment that allows us to reproduce these opacities reasonably well.

We consider a mixture of amorphous laboratory silicates (Dorschner et al. 1995, Mg_0.7Fe_0.3SiO₃) with amorphous carbon (Zubko et al. 1996, BE-sample). Pure laboratory silicates are “glassy” particles with almost negligible absorption cross sections in the near-IR, just where T Tauri stars are most luminous. Such particles would rather scatter the incident light away from the disk, which would lead to substantial problems in explaining the near-IR excess of T Tauri stars. Therefore, the inclusion of a conductive, hence highly opaque, albeit featureless material in the near-IR is necessary. However, it is unclear whether amorphous carbon, metallic iron, or e.g. troilite (FeS) should be used. Our simplified material composition is inspired by the past disk modelling expertise of the team, and by the solar system composition proposed by Min et al. (2011).

The dust grains are assumed to be composed of 60% silicate, 15% amorphous carbon, and 25% porosity, by volume, well-mixed on small scales. The effective refractory index of this porous material is calculated by applying the Bruggeman (1935) mixing rule. In contrast, just adding opacities (assuming separate grains of pure materials with the same size distribution and equal temperatures) does not seem physically justified at all, and does not reproduce the optical properties of aggregate particles equally well (Min et al. 2016).

We use a distribution of hollow spheres (DHS) with a maximum hollow volume ratio ${\mathit{V}}_{\mathrm{hollow}}^{\max }\mathrm{=}\mathrm{0.8}$. This approach avoids several artefacts of Mie theory (spherical resonances) and can account for the most important shape effects, see details in (Min et al. 2016). In combination with our choice of the amorphous carbon optical constants from Zubko et al. (1996, BE-sample), this approach captures the “antenna-effect” observed from the aggregate particles, where irregularly shaped inclusions of conducting materials result in a considerable increase of mm-cm absorption opacities.

Table 3 summarises our standard choices of dust parameter values, and Fig. 2 shows the resulting opacities including scattering and extinction. Figure 3 shows the dependencies of the dust absorption opacities on the remaining free dust size and material parameters. We will continue to discuss these results and their impact on predicted continuum and line observations in Sect. 5.2.1. Our standard dust opacities feature

• a FUV-dust extinction opacity of about 1000 cm²/g(dust), which is about 100 times less than in standard ISM models with MRN size distribution (Mathis et al. 1977);

• a dust albedo of 64% at FUV wavelengths, and 58% at 1 μm;

• a dust absorption cross section of about 5.8 cm²/g(dust) at 850 μm, which is about a factor of 1.6 larger than the value of 3.5 cm²/g(dust) used by Andrews & Williams (2005) to determine disk masses from sub-mm fluxes²;

• a millimetre dust absorption slope of about 1; and

• a centimetre dust absorption slope of about 1.5.

The dust opacities described above have been determined from the well-mixed dust size distribution function f₀(a), see Sect. 3.4. Since we have dust settling in the disk, the dust opacities need to be computed from the local settled dust size distribution function f(a,r,z), and will hence not only depend on λ, but also on the location in the disk (r,z). Simply put, dust settling leads to a strong concentration of the (sub-)mm opacity in the midplane, in particular in the outer disk regions, but only to a mild reduction of the UV opacities in the upper and inner layers. According to Fig. 3, the main effects are:

• 1)

  Large a_min values reduce the optical and UV opacities, and destroy the 10 μm/20 μm silicate emission features.

• 2)

  Large a_max values reduce the UV, optical, near-IR and far-IR opacities considerably, because the available dust mass is spread over a larger size range, and the big particles do not contribute much to the dust opacities at those wavelengths. Here, increasing a_max has similar consequences as lowering the dust/gas mass ratio. At longer wavelengths, a_max determines where the final transition to the Rayleigh-limit takes place (λ_trans ≈ 2πa_max), beyond which (λ>λ_trans) the opacity changes to a steeper slope. For larger a_max, this transition occurs at longer wavelengths. The choice of our reference value of a_max = 3 mm is motivated by cm-observations which usually do not show such breaks. However, Greaveset al. (in prep.) report on first evidence of such a steepening toward cm wavelengths after removal of the free-free emission component, based on new Green Bank Telescope data up to 8.6 mm.

• 3)

  The powerlaw size index a_pow determines the mixing ratio of small and large grains. Since the smaller particles are responsible for the short wavelength opacities, and the large grains for the long wavelength opacities, a_pow determines the general opacity slope, and the mm and cm-slopes in particular.

• 4)

  A large volume fraction of amorphous carbon reduces the 10 μm, 20 μm silicate emission features, fills in the opacity deficits of the major solid-state silicate resonances (up to about 8 μm), and flattens the absorption opacities at millimetre and centimetre wavelengths.

In order to facilitate the adoption of these opacities in other work, a Fortran-90 package to compute the DIANA standard dust opacities³.

3.8. PAHs

Polycyclic aromatic hydrocarbon molecules (PAHs) play an important role in our disk models via (i) continuum radiative transfer effects; (ii) photoelectric heating of the gas; and (iii) chemical effects. The chemical effects include the charging of the PAHs, the release and consumption of free electrons via photo-ionisation and recombination, and further effects due to charge exchange reactions. We study (i) and (ii) for neutral PAHs in this paper, but do not include the chemical effects as we are using here the small DIANA chemical standard. In contrast, the large DIANA chemical standard (Paper II) has the PAHs included in the selection of chemical specimen, and therefore accounts for (iii) as well.

PAHs are observed via their strong mid-IR emission bands in many Herbig Ae/Be stars (e.g. Maaskant et al. 2014), whereas detection rates in T Tauri stars are much lower (Geers et al. 2006), possibly because T Tauri stars generate much less blue and soft UV stellar radiation to heat the PAHs. PAHs in Herbig Ae/Be disks seem to have sizes of at least 100 carbon atoms (Visser et al. 2007). The PAH abundance in the disk is assumed to be given by the standard abundance in the interstellar medium (Tielens 2008), modified by factor f_PAH$\frac{{\mathit{n}}_{\mathrm{PAH}}}{{\mathit{n}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-7}\hspace{0.17em}{\mathit{f}}_{\mathrm{PAH}}\hspace{0.17em}\frac{\mathrm{50}}{{\mathit{N}}_{\mathrm{C}}}\mathrm{·}$(6)f_PAH = 1 corresponds to the interstellar medium (ISM) standard⁴. Here, n_PAH [cm^-3] is the PAH particle density, n_{⟨ H ⟩} [cm^-3] is the hydrogen nuclei density and N_C is the number of carbon atoms in the PAH. The actual PAH abundance in disks is disputed (e.g. Geers et al. 2006; Visser et al. 2007). Values f_PAH ≈ 0.1 or lower seem typical in Herbig Ae disks (Geers et al. 2006).

We assume N_C = 54 carbon atoms and N_H = 18 hydrogen atoms (“circumcoronene”) in the reference model, resulting in a PAH mass of 667 amu and a PAH radius of 4.87 Å (Weingartner & Draine 2001). N_C and f_PAH are free model parameters, as well as a decision whether to select the neutral or charged PAH opacities. Circumcoronene IR and UV spectra have been directly measured by Bauschlicher Jr. & Bakes (2000). However, in the DIANA framework, we use “synthetic” PAH opacities of neutral and charged PAHs are calculated according to Li & Draine (2001) with updates from Draine & Li (2007), including the “graphitic” contribution in the near-IR and the additional “continuum” opacities of charged PAHs.

In comparison to the low UV opacities of evolved dust in disks (Sect. 3.7), PAHs can easily dominate the blue and UV opacities, see Fig. A.2. This happens in the well-mixed case for f_PAH ≳ 0.1 in our disk models. The dominance of the PAH opacities in the UV is even stronger in the upper disk regions because of dust settling (we assume that PAH molecules do not settle). For Herbig Ae disks, where the maximum of the stellar radiation is released around 400 nm, f_PAH ≳ 1 would imply that the stellar photons are predominantly absorbed by the PAHs rather than by the dust. The absorbed energy would then be re-emitted via the strong PAH mid-IR resonances, and it is this mid-IR PAH emission that would predominantly heat the disk. Furthermore, the stellar UV usually reaches the line forming regions in a disk indirectly, via scattering on dust particles from above. For f_PAH ≳ 0.1, the PAHs would effectively shield the disk from UV radiation, because UV scattering by PAHs is extremely inefficient. These two effects have large implications on our models for both the internal dust and gas temperature structure in a disk.

The treatment of PAHs in the Monte Carlo programs MCFOST and MCMax is standard, using a quantum heating formalism with stochastic PAH temperature distribution. This mechanism was first proposed for small grains by Desert et al. (1986) and later applied to PAHs in the interstellar medium by Manske & Henning (1998), Guhathakurta & Draine (1989), Siebenmorgen et al. (1992). The Monte Carlo programs offer additional options to take into account e.g. a PAH size distribution and an internal determination of the PAH charge, by balancing the basic photo-ionisation and recombination rates (Maaskant et al. 2014). However, these options involve some quite specific simplifications, for example no negative PAHs, no charge exchange reactions, and an assumed electron concentration, which we need to avoid for reasons of consistency for the DIANA modelling efforts.

While looking for a fast, simplified, and robust way to treat the most important effects of PAHs equally well in all our disk models, we discovered that a simplified treatment of the PAHs in radiative equilibrium, according to $\begin{array}{ccc}& & \\ & & \\ & & \end{array}$leads to quite accurate results in comparison to the stochastic PAH treatment, see Appendix B and Fig. A.3. Here, S_ν is the source function, ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}$ and ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{sca}}$ are the absorption and scattering opacities (of dust and PAHs as annotated), B_ν(T) is the Planck function and J_ν is the mean intensity. The scattering term ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{sca}}{\mathit{J}}_{\mathit{\nu }}$ is here simplifyingly written for isotropic scattering. Equations (8) and (9) express the conditions of radiative equilibrium with separate dust and PAH temperatures, T_dust and T_PAH, respectively.

The quantum heating formalism is appropriate for ISM conditions, where PAHs are only sometimes heated by rare FUV photons. In contrast, we show in Appendix B that the PAHs in the inner regions of protoplanetary disks, which are responsible for the observable mid-IR PAH emission features, are situated in an intense optical and infrared radiation field created by the star and by the dust and the PAHs in the disk, which keeps the stochastic PAH temperature distribution high and quite narrow, i.e. close to the analytical treatment in radiative equilibrium. Li & Mann (2012) found similar results for nano grains acquiring an equilibrium temperature when exposed to intense starlight.

3.9. Chemistry and heating/cooling balance

Based on the results of the continuum radiative transfer as described in the previous sections, the gas phase and ice chemistry is calculated in kinetic chemical equilibrium, coupled to the gas heating/cooling balance. These parts of the model have been described elsewhere, see (Woitke et al. 2009) for the basic model, (Thi et al. 2011) for continuum radiative transfer, (Woitke et al. 2011) for updates concerning non-LTE treatment and heating and cooling, and (Aresu et al. 2011) for X-ray heating and chemistry, and are not discussed in this paper. In this paper we use the small chemical network, as proposed in Paper II. We carefully select 100 gas phase and ice species (see Paper II), and take into account altogether 1288 reactions. All 1065 gas phase and UV reactions among the selected species are taken into account from the UMIST 2012 database (McElroy et al. 2013), including the old collider reactions. We replace the treatment of photo-reactions by individual photo cross-sections from the Leiden Lamda database (Schöier et al. 2005) where possible, add 145 X-ray reactions (Aresu et al. 2011), 40 ice adsorption and thermal, UV photo and cosmic ray desorption reactions, and 38 auxiliary reactions including those of vibrationally excited molecular hydrogen H$\begin{array}{c}\mathit{\star }\\ \mathrm{2}\end{array}$, see details in Paper II. We take the (gas + ice) element abundances from Table 2 in Paper III. Appendix C discusses the validity of our approach to use the time-independent solution of our chemical rate network to compute the chemical composition of the disk and the gas emission lines.

3.10. Line radiative transfer

After the continuum radiative transfer, gas and dust temperature structure, chemistry and non-LTE level populations have been determined, a formal solution of line and continuum radiative transfer is carried out in 3D, using a bundle of 356 × 144 parallel rays towards the observer at distance d under inclination angle i, see (Woitke et al. 2011, Appendix A.7 therein) for details. These computations result in observable quantities like line fluxes, line velocity-profiles, molecular line maps and channel maps.

4. Model implementation and verification

The DIANA standard modelling assumptions summarised in Sect. 3 (stellar and interstellar irradiation, disk shape, dust opacities, dust settling, treatment of PAHs) have been implemented into MCFOST (Pinte et al. 2006, 2009) MCMax (Min et al. 2009) and ProDiMo (Woitke et al. 2009). The independent implementation of our modelling assumptions has allowed us to perform stringent checks on our computational methods and numerical results. Figure 2 shows a validation of our dust opacity implementation. Figure A.2 compares the assumed gas calculated settled dust densities, the resulting dust and PAH temperatures, as well as the SEDs. Apart from some minor temperature deviations in the optically thick midplane regions, which are irrelevant for the predicted observations, we achieve an excellent agreement concerning the physical state of the disk and all predicted observations. In particular, the upper right part of Fig. A.2 shows that we obtain very similar SED and spectral shape of the PAH features no matter whether we useProDiMo, MCFOST, or MCMax.

Further verification tests (not shown here) have been undertaken for disk models with gaps, where the numerical resolution of the inner wall of the outer disk is particularly important, and for MC → ProDiMo “chain models”. In these chain models, we use the Monte Carlo codes to compute the disk structure, the dust and PAH temperatures, and the internal radiation field J_ν(r,z), and then pass these results on to ProDiMo to compute the gas temperature structure, the chemical composition of ice and gas, and the emission lines.

The advantages of using the MC → ProDiMo chain models are (i) the Monte-Carlo technique is computationally faster; (ii) the temperature iteration scheme is more robust, in particular at high optical depths; and (iii) the Monte-Carlo technique allows for a more detailed implementation of radiation physics, in particular anisotropic scattering, PAHs with stochastic quantum heating, and polarisation. For the effects of anisotropic scattering, see Fig. G.2 in Appendix G. These more sophisticated options are not used in this paper, in order to facilitate comparisons to the results obtained with pureProDiMo. The pre-existing interface between MCFOST andProDiMo (Woitke et al. 2010) has been generalised and implemented in MCMax, such that now all MCFOST and MCMax users are able to use ProDiMo to predict chemical and line results on top of their continuum models.

Appendix F discusses the numerical convergence of our models as function of numerical grid resolution, for both the pure ProDiMo and the MC → ProDiMo chain models. Our conclusion here is that we need about 100 × 100 grid points in both ProDiMo and MC models, to achieve an accuracy <10% for all continuum observables and line flux predictions.

5. Results

The results of our disk models are presented in the following way. We first introduce a simple single-zone reference model in Sect. 5.1 which roughly fits a number of typical continuum and line observations of Class II T Tauri stars. In the following two sub-sections, we then study the impact of our model parameters on the various continuum and line observables by looking at how our model predictions change with respect to the reference model. In Sect. 5.2, we study the impact of selected model parameters on all observables at a time, and in Sect. 5.3, we discuss particular observables separately.

5.1. The reference model

Table 3 summarises our model parameters, and lists the values used for the reference model. The resulting spectral energy distribution (SED), visibilities, line observations, and some integrated properties are shown in Fig. 4. Concerning the integrated properties, we calculate the mean gas temperature in the disk ⟨ T_gas ⟩, the mean dust temperature ⟨ T_dust ⟩, the near-IR excess, the 10 μm SED amplitude, the mm-slope and the cm-slope as explained in Appendix D, see Eqs. (D.1) to (D.8). The reference model is characterised by

• a near-IR excess of about 0.12 L_⊙;

• clearly visible silicate dust emission features around 10 μm and 20 μm;

• a descending SED-slope dlog (νF_ν)/dlog λ< 0 beyond 20 μm, as is typical for continuous (i.e. non-transitional) T Tauri disks;

• a 1.3 mm continuum flux of 60 mJy with an apparent radius (semi-major axis) of about 100 AU (0.75″ at a distance of 140 pc) – typical observed values are about 20−200 mJy, and 0.25′′−1.4′′, see (Guilloteau et al. 2011);

• a mm-slope of about 2.4 – typical observed values are about 1.9−2.7, see (Ricci et al. 2010, 2012);

• a [OI] 63 μm line flux of 25 × 10^-18 W/m² – typical observed values for non-outflow sources are about (3−50) × 10^-18 W/m², see (Howard et al. 2013);

• a ¹²CO J = 2 → 1 line flux of about 15 Jy km s^-1, and a ¹²CO/¹³CO line ratio of about 5 – for typical values, see (Williams & Best 2014);

• an apparent radius (semi-major axis) in the ¹²CO J = 2 → 1 line of about 450 AU (3.5″ at 140 pc), typical observed values are about 1″−5″, see (Williams & Best 2014);

• weak CO ro-vibrational lines with a broad, box-shaped emission profile mostly emitted from the far side of the inner rim, which are not very typical with respect to observations. The “central nose” on top of the averaged line profile is a contribution from low-J lines which are also emitted from more distant disk regions, see Fig. 17. For typical CO ro-vibrational observations, see Sect. 5.3.7.

All other emission lines in the IR to far-IR spectral region are rather weak, which would likely result in non-detections with current instruments (for example o-H₂O 63.3 μm, CO J = 18 → 17, o-H₂ 17.03 μm), maybe except for the optical [OI] 6300 Å line (model flux 7 × 10^-18 W/m²). The CO fundamental ro-vibrational lines are also rather weak in the reference model (of order 2 × 10^-18 W/m² at FWHM = 130 km s^-1), which is below the detection limit of e.g. the CRIRES spectrograph (about 10^-18 W/m² at FWHM = 20 km s^-1).

5.2. Impact of selected model parameters

5.2.1. Dust size and opacity parameters

It is important to realise that dust grains in protoplanetary disks are likely to be very different from the tiny dust particles in the diffuse interstellar medium (ISM) for which the astronomical silicate opacities have been constructed by Draine & Lee (1984), only considering λ< 1 μm, and using MRN (Mathis et al. 1977) size parameters. In contrast, we expect the dust grains in disks to be much larger, up to mm-sizes, which reduces the UV dust opacities by a large factor (about 100) depending on parameters a_min, a_max and a_pow, see Table 2 and Fig. 3.

This simple and straightforward fact distinguishes our disk models from other chemical models (Table 1). In our models, the dust is much more transparent in the UV, allowing the UV to penetrate deeper into the disk, which increases the importance of molecular self-shielding, and reduces the importance of X-rays relative to the UV.

Table 2 shows that a disk-typical dust size distribution can be expected to have additional substantial impacts on disk chemistry, see also Vasyunin et al. (2011). The total grain surface area per H nucleus, important for surface chemistry, H₂ formation and photoelectric heating, is reduced by a factor of about 250 with respect to the MNR dust model, and the dust particle concentration n_d/n_{⟨ H ⟩} is only of order 10^-14, where ${\mathit{n}}_{\mathrm{d}}\mathrm{=}{}^{\mathrm{\int }}\mathit{f}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{a}$ is the (unsettled) total dust particle density. This implies, for example, that even if every dust grain was negatively charged once, there would be hardly any effect on the midplane electron concentration.

Figure 5 shows the results of the model when switching to uniform 0.1 μm sized dust particles and astronomical silicate opacities. The SED is now featured by stronger 10 μm/20 μm silicate emission features, higher far-IR continuum fluxes, and a quite sudden kink around 200 μm, followed by a steeper decline toward millimetre wavelengths. The apparent size of the disk is smaller at 1.6 μm, but larger at 10 μm and 1.3 mm. The disk is now warmer in dust, but cooler in gas, in fact mean dust and gas temperatures are more equal. Most emission lines in the model with uniform 0.1 μm dust particles show weaker fluxes, by up to a factor of ten. However, the CO J = 10 → 9 line doesn’t follow this general trend.

Understanding the impact of the dust size and material parameters on gas temperature and emission lines can be isolated to the effect of a single quantity, namely the dust optical depths at UV wavelengths τ_UV, see Sects. 5.3.4, 5.3.5 and 5.3.7. All dust parameter alterations that result in lower τ_UV will generally lead to a deeper penetration of UV into the disk, causing an increase of the thickness of the warm molecular gas layer, and this leads to stronger emission lines at optical to far-IR wavelengths. The impact is less pronounced on (sub-)mm lines, see Sect. 5.3.6, although secondary temperature and chemical effects are important to understand the (sub-)mm line ratios.

5.2.2. Hydrostatic disk models

Figure 6 summarises the results of two variants of (1+1)D hydrostatic models, where the vertical disk extension at any radius is computed from the condition of hydrostatic equilibrium, which requires an iterative approach of radiative transfer, chemistry and gas heating/cooling balance, see (Woitke et al. 2009) for details. The lower model is the proper hydrostatic solution, where the pressure p = ρkT_gas/μ is calculated according to the local gas temperatures T_gas and mean molecular weight μ resulting from chemistry. The upper model is a simplified version thereof, where the dust temperature is used instead (T_gas ≈ T_dust), and the mean molecular weight is assumed to be constant (μ ≈ 2.3 amu). The observable properties of these hydrostatic models are as follows.

• The SEDs of both types of hydrostatic models cannot explain the observed levels of near-IR excess for T Tauri stars, because the inner rim is quite low.

• Between 20μm and 50μm, the SED displays an  increasing  slope, which is caused by the strong flaring of the outer disk (see also Fig. 2 in Meijerink et al. 2012).

• The models are substantially warmer, both in gas and dust, as compared to the reference model, again because of the flaring of the outer disk.

• The far-IR to mm emission lines are all stronger as compared to the reference model. The [OI] 63.2 μm flux is (7−10) × 10^-17 W/m², which is out 4 times stronger as in the reference model, and quite high with respect to observations. The ¹²CO/¹³CO line ratio is as large as 8.

• The ro-vibrational CO ro-vibrational lines are weaker, because of the low inner rim, and have complicated multi-component line profiles.

The two variants of hydrostatic models have quite similar observable properties, including the spectral lines, although the disk gas density structure is remarkably different. The proper hydrostatic model displays larger amounts of extended hot gas high above the disk in the inner regions. However, this highly extended hot gas is purely atomic, hence it does not emit in molecular lines, one has to use atomic tracers to detect it, such as [OI] 6300 Å or maybe [NeII] 12.82 μm. In summary, hydrostatic passive disk models have some issues explaining the observed SED and line properties of T Tauri stars.

5.2.3. Disk flaring and/or dust settling?

Little disk flaring and strong dust settling have similar effects on the SED, see Fig. 7. The β = 1.05 and α_settle = 10^-4 models have practically indistinguishable SEDs beyond 20 μm, with a more steeply decreasing slope as compared to the reference model around ~ 50 μm. The reason for that steeper slope is that a more self-shadowed dust configuration leads to less interception of star light per disk radius interval, hence to a steeper decline of the dust temperatures as function of radius (e.g. Beckwith et al. 1990; Chiang & Goldreich 1997). We can see from Fig. 7 that indeed the mass averaged dust temperature ⟨ T_dust ⟩, which is dominated by the outer disk regions, has fallen from about 19 K in the reference model to about 13 K in both cases. Lacking disk flaring and strong dust settling both cause very low dust temperatures in the midplane, of order 4 K already at r = 150 AU in the α_settle = 10^-4 model, only limited by CMB and other background radiation. These very cold disk models have particular mm and cm-properties, because the cold dust is not entirely emitting in the Rayleigh-Jeans limit even at millimetre wavelengths, see Sect. 5.3.3.

Since the scale height is anchored at r = 100 AU in the model, little flaring (the β = 1.05 model) implies a tall inner disk which causes a strong near-IR excess. Assuming strong dust settling instead avoids these artefacts, because the impact of dust settling is strongest where the gas densities are lowest, i.e. in the outer regions, whereas the inner regions are only affected a little.

When looking at the impact on gas temperature and emission lines, however, the two models show just opposite effects. Lacking disk flaring moves the gas into the disk shadow, causing lower gas temperatures and weaker emission lines. Dust settling, in contrast, leaves the gas bare and exposed to the stellar UV radiation, leading to higher gas temperatures and stronger gas emission lines in general.

Therefore, in order to diagnose lacking disk flaring and/or strong dust settling, the far-IR SED slope around 50 μm is crucial, but to distinguish between disk flaring and dust settling, the simultaneous observation of far-IR gas lines is the key.

5.3. Parameter impact on selected observables

5.3.1. SED

Figure 8 shows the impact of our model parameters on the calculated spectral energy distribution (SED). Some parameter dependencies have already been discussed and explained in Sects. 5.2.1 and 5.2.3, but we repeat the essence here to give a comprehensive overview of all important effects.

The dust mass M_dust shifts the SED up and down at long wavelength, where the disk is predominantly optically thin. Its influence diminishes at λ ≲ 100 μm, but even at 20 μm, where the disk is massively optically thick, a change of M_dust still produces noticeable changes. This is because more mass increases the height at which the disk becomes radially optically thick, which has similar consequences as increasing the scale height.

The reference scale height H₀ affects the SED at all shorter wavelength λ ≲ 200 μm, but not the Rayleigh-Jeans tail of the SED. Larger scale heights mean to intersect more star light, and to produce a warmer disk interior which re-emits more thermal radiation from all optically thick disk regions.

The flaring index β rotates the SED around a point at λ ≈ 20 μm here, depending on the model and on the choice of the reference radius H₀ in Eq. (2). Large β values mean that we have a flared disk with a low inner rim but with tall outer regions, which produce less near-IR but more far-IR excess. Small β lead to a “self-shadowed” disk structure with very cold dust in the outer parts.

The inner radius R_in regulates the maximum temperature of the dust grains at the inner rim. Larger R_in therefore result in less near-IR emission (“transitional disks”). However, the total amount of excess luminosity is not changing. For large R_in, the luminosity excess merely shifts from the near-IR to the mid-IR region, and beyond.

Dust settling, with parameter α_settle describing the strength of the turbulent mixing, almost exclusively affects the long wavelength parts of the SED (λ ≳ 20 μm). According to the Dubrulle prescription (see Eq. (5)), dust settling is much more effective at large radii where the densities are low, in which case the dust grains cannot be easily dragged along turbulent gas motions. Consequently, the outer disk parts become flat as seen in dust, although the gas still extends high up. Therefore, strong settling has similar consequences as lacking disk flaring at long wavelengths. For a well-mixed dust/gas mixture, the mm-grains tend to cover all spectral features produced by the small grains with their flat, greyish opacity, washing out the 10 μm and 20 μm silicate emission features. Dust settling removes the large grains from the disk surface, and therefore amplifies the silicate emission features, which seems necessary in many cases to re-produce the observed shape of the silicate emission features.

The maximum dust particle size a_max has a similar influence as M_dust at long wavelengths. Increasing a_max effectively means to put more dust mass into very large particles which have almost no opacity at shorter wavelengths. However, beyond about 1 mm in this model, where the largest particles do provide the dominating opacities, the SED starts to change slope depending on the value of a_max.

The dust size distribution powerlaw index a_pow regulates the mixture of small and large dust particles in the disk. It thereby changes, in particular, the mm and cm-slopes. Larger a_pow values also amplify the 10 μm and 20 μm silicate emission features, because the grayish opacity of the large grains is mostly removed from the model.

The volume fraction of amorphous carbon has a surprisingly large impact on the SED at all wavelengths. As discussed in Sect. 3.7, pure laboratory silicates are very effective scatterers, keeping the stellar radiation out of the disk, but they will hardly absorb it. Therefore, disks made of pure silicates are much cooler and emit less near-IR and far-IR excess. The fraction of amorphous carbon also changes substantially the mm and cm slopes through opacity effects.

The surface density powerlaw index ϵ has practically no influence on the SED, same with the tapering-off radius R_tap (not depicted), the outer radius R_out (not depicted, see also Bouy et al. 2008) and the minimum dust size a_min (as long as a_min ≲ 0.5 μm, not depicted). The dependencies of the SED on those non-depicted parameters are less than those shown for ϵ.

From the shape of the SED changes caused by the nine parameters depicted in Fig. 8, one can easily imagine how degenerate pure SED fitting can be (e.g. Robitaille et al. 2007), just consider, for example, a combination of lower dust mass with more amorphous carbon.

5.3.2. Visibilities

Figure 9 shows the impact of model parameters on the calculated visibilities at 1.6 μm (e.g. PIONIER), at 10 μm (e.g. MIDI), and at 1.3 mm (e.g. CARMA, ALMA). The mm-visibility (see e.g. Guilloteau et al. 2011) probes the apparent spatial extension of the disk, limited by minimum optical depth requirements to produce a detectable signal at those wavelengths. This apparent size is most directly influenced by the tapering-off radius R_tap. The “sharpness” of the outer edge γ is reflected by the steepness of the V²-decline. However, more dust in the outer regions (larger M_dust, smaller ϵ) also increases the optical depths in these regions, which leads to larger apparent sizes as well. Strong dust settling α_settle → 10^-4 moves the grains toward the midplane into the disk shadow where they are substantially cooler (see Sect. 5.3.3 and Fig. 7), so cold that some fraction of the dust grains does not emit in the Rayleigh limit at 1.3mm, thus producing less extended flux, which leads to a smaller apparent size⁵.

The 10 μm visibilities reflect the radial extension of warm dust in the disk surface layer producing the 10 μm silicate emission feature (about 1 AU in the reference model). This extension is larger for warm, e.g. flared disks. In contrast, parameter choices which lead to cooler conditions at 1 AU cause a smaller appearance of the disk at 10 μm. Such parameter choices include smaller scale heights H₀, dust size parameter variations that increase the mean dust size (larger a_min, smaller a_pow), and lacking amorphous carbon. Dust settling plays no significant role here.

The 1.6 μm visibilities are more difficult to understand, see (Anthonioz et al. 2015). They have three components: the star, scattering and emission from the inner rim, and extended scattering. The extended scattering leads to a slight tilt of the V²-curves beyond baselines [kλ] ≳ 100, before V² drops to much lower values at baselines which corresponds to the inner rim of the disk (0.07 AU in the model). At even longer baselines, the interferometer would start to resolve the star (R_⋆ = 0.0097 AU), and the visibility would drop sharply, but such long baselines are currently not accessible, and not included in Fig. 9). The 1.6 μm visibilities hence probe the relative contributions of these three components and their spatial extensions. Most remarkably is the influence of flaring and settling, which powers/suppresses the extended scattering component, and the fraction of amorphous carbon which changes the albedo of the dust particles. Pure silicate dust particles (amC = 0), for example, are almost perfect scatterers at 1.6 μm, leading to a much more pronounced extended scattering component.

Noteworthy, the hydrostatic disk models show a stronger extended scattering component as well (because of the strong flaring of the outer disk), and less contributions from the inner rim, which has a lower wall height.

The inner rim radius R_in directly affects the second component of the 1.6 μm visibilities directly, namely the emission and scattering from the inner rim. Large R_in can also limit the radial extension of the 10 μm emission region from the inside, introducing new small scales in form of the ring thickness and the apparent height of the inner rim wall, with sharp edges, which leads to more complex visibility shapes.

5.3.3. The mm-slope and cm-slope

Figure 10 shows the impact of the model parameters on the observable SED slopes at millimetre and centimetre wavelengths, as defined by Eqs. (D.7) and (D.8), which are important diagnostics of grain growth in protoplanetary disks, see e.g. (Natta et al. 2007) and (Testi et al. 2014). The SED slopes are expected to reflect the dust absorption opacity slopes, with some flattening due to optical depth effects ${\mathit{\alpha }}_{\mathrm{SED}}\mathrm{\approx }\mathrm{2}\mathrm{+}\frac{{\mathit{\beta }}_{\mathrm{abs}}}{\mathrm{1}\mathrm{+}\mathrm{\Delta }}\mathit{,}$(10)Eq. (10) was derived by Beckwith et al. (1990)⁶ for a powerlaw surface density structure Σ ∝ r^{− p} and a vertically isothermal powerlaw temperature distribution T ∝ r^{− q}. This formalism was later relaxed by Ricci et al. (2010, 2012), who determined T(r) according to Chiang & Goldreich (1997) and considered a self-similar disk with tapered outer edge. Ricciet al. found values α_SED ≈ 1.9−2.7 for the Taurus-Auriga region, which are significantly lower than what is expected from small interstellar grains β_abs ≈ 1.7 (Draine 2006), suggesting that the dust in protoplanetary disks must have much smaller β_abs ≈ 0.3−1.0 at (1−3) mm, indicating dust growth.

The standard DIANA opacities have ${\mathit{\beta }}_{\mathrm{abs}}^{\hspace{0.17em}\mathrm{mm}}\mathrm{\approx }\mathrm{1.0}$ and ${\mathit{\beta }}_{\mathrm{abs}}^{\hspace{0.17em}\mathrm{cm}}\mathrm{\approx }\mathrm{1.5}$, thus the SED slopes of the reference model are expected to be ${\mathit{\alpha }}_{\mathrm{SED}}^{\mathrm{mm}}\mathrm{=}\mathrm{3.0}$ and ${\mathit{\alpha }}_{\mathrm{SED}}^{\mathrm{cm}}\mathrm{=}\mathrm{3.5}$ in the optically thin limit Δ = 0. However, the reference model exhibits ${\mathit{\alpha }}_{\mathrm{SED}}^{\mathrm{mm}}\mathrm{\approx }\mathrm{2.4}$ and ${\mathit{\alpha }}_{\mathrm{SED}}^{\mathrm{cm}}\mathrm{=}\mathrm{3.35}$, in agreement with observations, suggesting optical depths corrections of Δ_mm ≈ 1.5 and Δ_cm ≈ 0.1, respectively.

Closer inspection shows, however, that these derived Δ do not agree at all with the expected flattening due to optical depth effects. The radius r₁ where the vertical dust optical depth equals unity ${\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{\left(}{\mathit{r}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{dust}\mathit{,}\mathrm{abs}}\mathrm{\left(}{\mathit{r}}_{\mathrm{1}}\mathit{,z}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{=}\mathrm{1}\mathit{,}$(11)is only r₁ ≈ 6.9 AU at λ = 1.3 mm and r₁ = 0.9 AU at λ = 7 mm in the reference model, which results in tiny corrections, Δ_mm = 0.04−0.12 and Δ_cm = 0.01−0.03, as derived from the equations in Beckwith et al. (1990), depending on what is assumed for the outer radius in our tapered-edged models. At both wavelengths, the expected Δ-corrections are too small, inconsistent with the results obtained from our radiative transfer models.

This conclusion holds for all models computed in this paper. The strongest optical depth effects occur in massive disks (M_disk = 0.1 M_⊙), if the mass is more concentrated toward the centre (ϵ = 1.5), if the disk is small (R_tap = 50 AU), and/or if the dust size and opacity parameters lead to larger mm-opacities. But even in all these cases, the expected Beckwith et al. Δ-corrections for optical depths effects at 1.3 mm stay well below unity, which is insufficient to explain the gentle mm-slopes obtained from our radiative transfer models. The mm-slopes from the computed SEDs are more gentle than expected, and optical depth effects are not the key to explain these discrepancies⁷.

Beckwith et al. (1990) derived Eq. (10) by assuming that all dust grains emit in the Rayleigh limit. If we ignore optical depth effects for a moment, the observable flux F_ν is exactly given by ${\mathit{F}}_{\mathit{\nu }}\mathrm{=}\frac{\widehat{{\mathit{\kappa }}_{\mathit{\nu }}}\mathrm{abs}}{{\mathit{d}}^{\mathrm{2}}}\mathrm{\int }{\mathit{\rho }}_{\mathrm{dust}}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{V}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{M}}_{\mathrm{dust}}}{{\mathit{d}}^{\mathrm{2}}}\hspace{0.17em}⟨{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}⟩\hspace{0.17em}\widehat{{\mathit{\kappa }}_{\mathit{\nu }}}\mathrm{abs}\mathit{,}$(12)where $\widehat{{\mathit{\kappa }}_{\mathit{\nu }}}{}^{\mathrm{abs}}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathit{/}{\mathit{\rho }}_{\mathrm{dust}}$ is the dust absorption coefficient per dust mass (assumed to be constant throughout the disk), $⟨{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}⟩\mathrm{=}{}^{\mathrm{\int }}\mathit{\rho }_{\mathrm{dust}}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{V}􏼮{}^{\mathrm{\int }}\mathit{\rho }_{\mathrm{dust}}\hspace{0.17em}\mathrm{d}\mathit{V}$ is the dust mass averaged Planck function, and M_dust is the total dust mass. The log-log slope α_SED = −∂log F_ν/∂log λ is then given by ${\mathit{\alpha }}_{\mathrm{SED}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }\log \hspace{0.17em}\mathrm{⟨}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\mathrm{⟩}}{\mathit{\partial }\log \mathit{\lambda }}\mathrm{+}{\mathit{\beta }}_{\mathrm{abs}}\mathit{,}$(13)Table 4 shows that the deviations of the Planck derivative from its limiting value of 2 can be substantial. At 1.3 mm, for example, the Planck slope is about 1.35 and not 2, if the grains emit at 10 K. At 7 mm, deviations ≳0.2 dex are still conceivable if the majority of grains would emit at 5 K. Indeed, using these deviations from Rayleigh-Jeans regime, Dutrey et al. (2014) report on vertical mean dust temperatures of 8.5 K at 300 AU in the disk of GGTau.

The dependencies of α_SED on the disk temperature structure, according to Eq. (13), can explain the results obtained from our radiative transfer models. The mean dust temperature according to Eq. (D.2) is ⟨ T_dust ⟩ ≈ 19 K in the reference model, but this is a linear mean, and the Planck function is highly non-linear at low temperatures ⟨ B_ν(T) ⟩ ≪ B_ν( ⟨ T ⟩). Simply put, a considerable part of the dust in the disk is so cold that it does not contribute significantly to the 1.3 mm flux. The minimum dust temperature in the reference model is about 4.5 K. The cold dust over total dust mass fraction is 0.12 (for T_dust< 7 K), 0.31 (for T_dust< 10 K), 0.58 (for T_dust< 15 K), and 0.75 (for T_dust< 20 K). It is about this fraction, with efficiencies according to Table 4, that is missing in the observable flux, causing the deviations from α_SED = 2 + β_abs.

These temperature effects explain the qualitative behaviour of the SED millimetre and centimetre slopes in the models as shown in Fig. 10. The dust size and material parameters impact the SED slope directly via changing the dust absorption slope β_abs, compare Fig. 3. These parameters have the strongest impact on the SED slope. The green disk shape parameters have an influence on the dust temperature structure in the disk T_dust(r,z) and, therefore, have an indirect influence on the SED slope. Figure 11 shows that all disk shape parameters that lead to very cold midplane conditions in the disk around r ≈ 50 AU are well correlated with the more gentle mm-slopes in the models, in particular little flaring (β = 1.05) and a steep column-density structure (ϵ = 1.5). However, dust settling has in fact an even stronger impact (α_settle = 10^-4, leftmost point in Fig. 10). By settling, the bigger grains are moved to the cooler midplane, and concentrating the grains toward the midplane also increases the shadow formation there, making the midplane even cooler.

Our conclusion is that, indeed, the dust absorption opacity slope β_abs mostly determines the SED slope. However, α_SED flattens significantly for cold disks, an effect that has not been reported so far and that seems more important that the classical Δ-correction for optical depth effects (Beckwith et al. 1990).

5.3.4. The [OI] 63.2 μm emission line

We now turn our attention from continuum observations to gas emission lines. We have selected four representative emission lines for the discussion in this paper, which are frequently observed, and which emerge from different disk regions, see Fig. 12. All selected lines have chemically robust carriers, namely the O-atom or the CO-molecule, which are not critically dependent on chemical details. These lines are rather influenced by the shape of the disk and gas temperature distribution in the disk surface layer.

The [OI] 63.2 μm line is usually the strongest disk emission line throughout the electromagnetic spectrum, compare Fig. 4, and has been detected in 84% of T Tauri stars with disk dust masses >10^-5 M_⊙ by the Herschel open time key program GASPS (Dent et al. 2013), including outflow sources (see Howard et al. 2013). With an excitation energy of 227 K and a critical density of ~6 × 10⁵ cm^-3, this line originates in disk layers even above the CO containing molecular layers, see Fig. 12. In order to excite this line, the gas needs to be as warm as ≳50 K (Kamp et al. 2010). These conditions are only present in outflows and in the disk up to radial distances ≲100 AU, above the molecular layers which are too cold because of molecular line cooling.

The impact of the model parameters on the predicted [OI] 63.2 μm line flux is shown in Fig. 13. According to the physical excitation mechanism explained in the previous paragraph, the line flux increases with all parameters that directly trigger the heating in the uppermost disk layers, marked in blue in Fig. 13, namely the stellar UV excess f_UV, the X-ray luminosity L_X, the PAH concentration f_PAH and the efficiency of exothermic reactions γ_chem. The disk shape parameters (marked in green in Fig. 13) also play an important role. For a self-shadowed disk (e.g. for flaring index β = 1.05), the distant oxygen gas is situated in the disk shadow, and the [OI] 63.2 μm line is substantially suppressed. The line is massively optically thick (τ_line ≈ 100 in the reference model at 100 AU), so what counts is the size of the disk surface area with gas temperatures T_gas ≳ 50 K, which depends on the disk mass and shape parameters H_g(100 AU), β and ϵ. In comparison, dust size parameters and inclination play no significant role for this line.

Figure 13 shows that 5 different model parameters are able to change the [OI] 63.2 μm line flux by at least a factor of 2 within their reasonable ranges of values, therefore, finding clear correlations of the [OI] 63.2 μm line flux with one of these parameters seems quite unlikely, which could explain why Meeus et al. (2012) and Aresu et al. (2014) have reported on negative results concerning such correlations.

Since the [OI] 63.2 μm line is optically thick, one would not expect the gas mass nor the gas/dust ratio to be important. However, for this particular line, there is an interesting energy conservation mechanism at work. The emission of [OI] 63.2 μm line photons is the dominant cooling process in the line emitting regions, and therefore, the line luminosity must roughly equal the integrated heating rate in the [OI] 63.2 μm line emitting volume. Consequently, the gas temperature in this volume can be expected to relax towards an equilibrium value where the [OI] 63.2 μm line cooling balances the total heating, henceforth denoted as “self-regulation mechanism”. The heating in the line emitting region is provided by a number of physical processes that absorb and thermalise fractions of the incoming UV and X-ray photon energies, namely X-ray Coulomb heating, heating by neutral carbon photo-ionisation, PAH heating via photo-effect, and heating by exothermic chemical reactions driven by UV and X-ray reactions. For larger disk masses, the solid angle (as seen from the star) of the distant disk regions that absorb the UV and X-ray photons increases, therefore, the [OI] 63.2 μm line flux increases. The impact of the gas/dust ratio is actually more significant, as there is a competition between dust versus gas absorption of UV photons. For larger gas/dust ratios, fewer UV photons are absorbed by the dust (and re-emitted as continuous far-IR radiation), hence more of the available UV flux is converted into gas heating.

As a consequence of these processes, more gas (less dust) generally leads to an increase of the size of the disk surface area where T_gas ≳ 50 K. More specifically, an increase of the gas mass causes the [OI] 63.2 μm line emission region to shift upwards (Kamp et al. 2010), which captures more of the impinging UV and X-ray photons. However, because of the self-regulating energy conservation mechanisms explained above, effects are rather modest. By varying the gas mass by a factor of 10, the [OI] 63.2 μm line only changes by factors of a few.

5.3.5. CO high-J emission lines

As an example for high-J CO emission lines, we have selected the CO J = 18 → 17 line at 144.8 μm with an excitation energy of 945 K and a critical density of ~2 × 10⁶ cm^-3. Due to its high excitation energy, gas temperatures T_gas ≳ 200 K are required to excite this line, and these conditions are only present in the CO surface layer inside of r ≲ 5 AU in the reference model. The line is optically thick, but approaches τ_line = 1 at the outer boundary of the line emitting region. The line is rather weak in the reference model, a factor of about 30 lower than the Herschel/PACS detection limit. In fact, detection rates of this line for T Tauri stars are rather low, about 41% of the 34 objects selected by GASPS (Dent et al. 2013)⁸, most of them being identified as outflow sources (see Howard et al. 2013).

The “impactogram” of the CO J = 18 → 17 line (Fig. 14) shows many similarities to the behaviour of the [OI] 63.2 μm line (Fig. 13), for example the direction of effects, but the scaling of the abscissa is different. We are now reporting on effects that can potentially change a line flux by one order of magnitude.

Since both line and continuum are optically thick, and level populations are close to LTE, we can use Eq. (E.3), see Appendix E, to estimate the line flux. This demonstrates that for these optically thick far-IR lines, it is the difference between gas and dust temperatures (in the disk surface layers between about 1 AU and 5 AU) that determines the line flux.

Besides the parameters directly involved in the gas heating (blue), we can conclude from Fig. 14 that also the dust size parameters play an essential role. All parameter changes that imply an increase of the mean dust particle size (orange) lead to a reduction of the dust UV opacity, see Sect. 5.2.1, and hence to an increase of the thickness of the layer where the gas temperature is substantially larger than the dust temperature. Precisely speaking, as will be explained in Sect. 5.3.7 and Fig. 20, it is the ratio τ_line/τ_UV which is important. Since the gas temperature structure is more or less fixed to τ_UV, more gas, or lower UV continuum optical depths, both result in larger CO J = 18 → 17 line fluxes. This also explains the dependence on dust settling.

Interestingly, when the inner disk radius is increased in the model to 10 AU, i.e. well beyond the outer radius of the line emitting region in the reference model, the line flux increases. In this case, instead of radially continuous line emission, the line is preferentially emitted from the distant inner wall, and the larger emitting area of that wall seems more relevant than the radial dilution of the impinging UV and X-ray photons, as long as this high-energy irradiation is sufficient to cause gas temperatures ≳200 K at the inner wall.

5.3.6. (sub-)mm CO isotopologue lines

The CO J = 2 → 1 isotopologue lines have excitation energies of about 17 K and critical densities ~7 × 10⁴ cm^-3. The impact of the model parameters on our computed CO isotopologue line fluxes is shown in Fig. 15. In the following, we discuss these results obtained by our 3D diagnostic line transfer computations (Sect. 3.10), by means of a few simplified equations, to demonstrate and understand the main dependencies we find. The ¹²CO line is optically thick in all models, the ¹³CO is optically thick in most models, and the C¹⁸O line is borderline optically thin, meaning that the radial distances R(τ_line = 1), up to which the lines are optically thick, depends on isotopologue. In the reference model, R(τ_line = 1) ≈ 450 AU for the ¹²CO line, ≈210 AU for the ¹³CO line, and ≈110 AU for the C¹⁸O line. The level populations connected to the CO isotopologue lines in the (sub-)mm regime result to be close to LTE, and the continuum is optically thin. Under these circumstances, the line fluxes are approximately given by ${\mathit{F}}_{\mathrm{line}}\mathrm{\approx }\mathrm{2}\mathrm{\Delta \Omega }\hspace{0.17em}\mathrm{\Delta }\mathit{\nu }\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left[}{\mathit{T}}_{\mathrm{gas}}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{line}}\mathrm{=}\mathrm{1}\mathrm{\left)}\mathrm{\right]}\mathit{,}$(14)see Eq. (E.2) and explanations in Appendix E. The mean gas temperature in the CO (sub-)mm line emission regions (see Fig. 12) results to be ⟨ T_gas ⟩ ≈ T_gas(τ_line = 1) ≈ 20−35 K, and varies only little throughout the presented models (all isotopologues). This temperature only affects the line fluxes in a linear way in the Rayleigh-Jeans limit. Therefore, the discussion of the CO isotopologue line fluxes ${\mathit{F}}_{\mathrm{line}}\mathrm{\propto }\mathrm{\Delta \Omega }\mathrm{\approx }\frac{\mathit{\pi }{\mathit{R}}^{\mathrm{2}}\mathrm{(}{\mathit{\tau }}_{\mathrm{line}}\mathrm{=}\mathrm{1}\mathrm{)}\cos \mathrm{\left(}\mathit{i}\mathrm{\right)}}{{\mathit{d}}^{\mathrm{2}}}\mathit{,}$(15)simplifies to a discussion of the radius up to which the CO lines are optically thick R(τ_line = 1) as follows.

• The ¹²CO line fluxes depend only little on all parameters which do not change R(τ_line = 1) significantly, in particular f_UV, L_X, f_PAH, dust size and settling parameters, scale height H₀ and flaring parameter β, as well as parameters only relevant for the inner disk, such as R_in.

• Most important are those model parameters which directly determine R(τ_line = 1), these are the tapering-off radius R_tap and the tapering-off exponent γ.

• The dependency on γ is remarkable. The reference model has a tapering-off radius of R_tap = 100 AU, a mean CO emission temperature of T_gas(τ_line = 1) ≈ 25 K, and a ¹²CO J = 2 → 1 line flux of 15.7 Jy km s^-1 (1.2 × 10^-19 W/m²). According to Eq. (14), this line flux corresponds to an emitting radius of R(τ_line = 1) = 450 AU, which is possible only because we assume a smoothly decreasing surface density structure beyond R_tap with γ = 1 in the reference model. If the disk has a much sharper outer edge (γ = −0.5) the mean CO emission temperature increases by 10 K, but the flux is down to 4.2 Jy km s^-1 (3.3 × 10^-20 W/m²), which then corresponds to a radius of only 200 AU. These numbers are in very good agreement with the radii R(τ_line = 1) that we can directly measure in the models, see Fig. 12.

• There are a number of model parameters which indirectly change R(τ_line = 1), among them the disk mass M_disk, the gas/dust ratio, and the column density powerlaw exponent ϵ. These parameters change the amount of CO gas in the outer regions, which implies that R(τ_line = 1) changes, too. These indirect influences on the CO line fluxes are stronger in models with exponential tapering-off. Models with a sharp outer edge do not show much of those effects.

To summarise, the ¹²CO (sub-)mm lines probe the conditions in the tapering-off distant gas above and around the disk, well beyond the radial zones emitting the continuum. These conclusions are robust, almost purely geometrical, because CO is such an abundant and robust chemical constituent of the disks. There is no need for any sophisticated chemical effects to understand these lines.

In contrast, the CO isotopologue line ratios (see r.h.s. of Fig. 15) are more difficult to understand. It has been suggested, e.g. by Williams & Best (2014) and Miotello et al. (2014), that the CO isotopologue line ratios are an excellent probe of the disk mass, in particular if rare isotopes like C¹⁸O and C¹⁷O can be observed with high sensitivity⁹. We see the dependence of isotopologue line ratios on disk mass clearly in our models, too, but we also see other parameter dependencies that can be equally important. When looking at line ratios, the major dependencies (like F_line ∝ R²) cancel, though not completely, and other secondary temperature and chemical effects come into play. The (sub-)mm CO isotopologue line ratios show the following effects:

• The dependence of the ¹²CO /¹³CO line ratio on disk mass is almost negligible in our models, consistent with Williams & Best (2014), because both lines are optically thick. The dependence of the ¹³CO/C¹⁸O line ratio on M_disk is more pronounced, because the C¹⁸O line is borderline optically thin.

• The asymmetric dependence of the ¹³CO/C¹⁸O line ratio on M_disk indicates that we need disk masses as low as 10^-3 M_⊙ to fully be in the optically thin limiting case for the C¹⁸O line (the ¹³CO line is optically thick in most models anyway), also evident from the missing dependence on inclination.

• The impact of the sharpness of the outer disk edge γ is remarkable. For a very sharp outer edge (γ = −0.5) we come close to a configuration where all three CO isotopologue lines remain optically thick until the disk ends abruptly, see Fig. 16, in which case the line ratios would be unity. However, there are secondary temperature-effects (the ¹³CO and C¹⁸O lines are emitted from deeper layers which are cooler) which always keep the isotopologue line ratios >1.

• The dependencies on the gas/dust ratio are much more pronounced than on M_disk, which shows that it is not simply the gas mass that counts. In fact, the two models denoted by M_disk = 0.001 M_⊙ and gas/dust = 10 have the same gas mass, they only differ in terms of their total dust mass (see footnote below Table 3). The more dusty model (gas/dust = 10) has a taller midplane shadow, hence cooler conditions even at higher disk layers, and more CO ice. Therefore, the ¹³CO line becomes borderline optically thin and its flux drops quickly, leading to a larger ¹²CO /¹³CO line ratio.

• Figure 15 shows a strong impact of disk flaring (parameter β) on the CO isotopologue line ratios, which again can be attributed to CO ice formation. The β = 1.05 model produces a self-shadowed disk configuration with very cold midplane conditions, favouring CO ice formation. The resulting lower gaseous CO column densities weaken in particular the ¹³CO line, hence the ¹²CO /¹³CO line ratio becomes larger.

• There are also noticeable dependencies on dust size parameters and dust settling, which are again related to CO ice formation. All dust size and opacity parameters which lead to an increase of the dust opacities around 1 μm lead to a more pronounced midplane shadow, cooler conditions around (100−300) AU, more CO ice, hence larger ¹²CO /¹³CO line ratios.

To summarise, the analysis of the (sub-)mm CO isotopologue line ratios requires sophisticated modelling where the sharpness of the outer disk edge, the strength and thickness of the disk midplane shadow impacting the CO ice formation, and vertical temperature effects all play a significant role. Spatially resolved line data with high S/N (e.g. ALMA) are required to disentangle these effects.

5.3.7. CO υ = 1 → 0 emission lines

The fundamental ro-vibrational CO υ = 1 → 0 emission lines are regularly detected in T Tauri stars (e.g. Najita et al. 2003; Salyk et al. 2011; Brown et al. 2013), in Herbig Ae/Be stars (Brittain et al. 2003; Blake & Boogert 2004; van der Plas et al. 2015), as well as in transition disks (Goto et al. 2006; Salyk et al. 2009; Pontoppidan et al. 2008). Surveys of CO ro-vibrational emission lines in young stars (e.g. Brown et al. 2013) detect ¹²CO emission in about ~80% of the objects. The line profiles are generally double-peaked, but, in many objects, in particular in T Tauri stars, line profiles can also be singly peaked which is usually interpreted in terms of a slow (few km/s) molecular disk wind (Pontoppidan et al. 2011; Bast et al. 2011; Brown et al. 2013). ¹²CO integrated line fluxes observed are of order 10^-17 W/m², and line widths are 12−200 km s^-1, mostly reflecting the expected distribution of inclination. Disks around early-type Herbig AeBe stars have narrower line profiles, on average, than T Tauri stars Brown et al. (2013).

The lines have excitation energies (3000−6000) K and critical densities¹⁰ of order (10¹²−10¹⁴ ) cm^-3, depending on temperature and depending on whether the gas is H₂-rich or atomic, see (Thi et al. 2013) for details. The low-J fundamental CO lines are always massively optically thick in our models, with vertical line centre optical depths of order 10⁴−10⁷ for r< 10 AU in the reference model. Figure 12 shows that the CO υ = 1 → 0 lines are vertically emitted by a region that extends radially to about 1 AU in the reference model, from a thin horizontal layer at the top of the warm molecular layer. In addition, the line can also be emitted from the far side of the directly illuminated inner rim, if the disk is seen under an inclination angle i> 0°. This contribution is not accounted for in Fig. 12, and may in fact domiante under certain circumstances.

In LTE, the CO gas will emit substantially in the fundamental lines if (i) the gas temperatures exceed about 500−1000 K; and (ii) the gas is significantly warmer than the local dust temperature, according to Eq. (E.3) in Appendix E${\mathit{F}}_{\mathrm{line}}\mathrm{\approx }\mathrm{2}\mathrm{\Delta \Omega }\hspace{0.17em}\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}{\mathit{B}}_{\mathit{\nu }}\mathrm{\right[}{\mathit{T}}_{\mathrm{gas}}\mathrm{(}{\mathit{\tau }}_{\mathrm{line}}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{]}\mathrm{-}{\mathit{B}}_{\mathit{\nu }}\mathrm{[}{\mathit{T}}_{\mathrm{dust}}\mathrm{(}{\mathit{\tau }}_{\mathrm{cont}}\mathrm{=}\mathrm{1}\mathrm{)}{\mathrm{\left]}}^{\mathrm{\right)}}\mathit{.}$(16)These two conditions are always fulfilled at the inner rim, or, more precisely speaking, in the thin hot surface layer that covers the inner rim facing the star, which forms a very thin radial photodissociation region (PDR). These emissions from the inner rim creates a broad, box-like “minimum CO emission profile” for all R-branch and P-branch lines (Fig. 17) with a hot characteristic emission temperature. However, the CO emission lines created this way would be quite faint and could not be detected with current instruments, because of the tiny solid angle ΔΩ occupied by the inner rims of T Tauri stars.

In order to create observable line fluxes, the model must fulfil the above stated two conditions for CO υ = 1 → 0 emission also at larger radii, which only occurs in some models, depending on the parameters important for the gas heating, flaring, and dust shielding (Figs. 18 and G.1). If the extended gas is sufficiently warm, we have a mixture of narrow (extended) cool CO emission with some broad, hot emission from the inner rim, (Fig. 17). The emissions from the extended regions are lacking strong high-J lines not only because the gas is cooler there, but also because the optical depths are smaller for the high-J lines, i.e. the high-J lines probe deeper layers which are cooler, and where the differences between gas and dust temperatures start to vanish (Fig. 20). The cool, extended contributions add a central component to the line profile which is often double-peaked, but not always (Fig. G.1).

The following results have been obtained with the default CO model molecule in ProDiMo with vibrational quantum numbers v ≤ 2 and rotational quantum numbers J ≤ 50, limited to 110 levels, see (Thi et al. 2013) for extended options. Figure 18 shows the impact of the model parameters on the R(10) line flux and characteristic CO emission temperature as measured by R(35)/R(10). Figure G.1 shows how the mean line profile changes when selected model parameters are varied. The main effects are as follows.

• Stronger CO line fluxes are connected with more centrally peaked line emissions coming from larger emitting areas.

• All parameters that increase gas heating will generally increase line flux, in particular the UV excess f_UV (Garufi et al. 2014) and the efficiency of chemical heating γ_chem.

• However, the X-ray luminosity has no significant influence on the fundamental CO emission lines.

• Disk flaring makes the CO lines stronger and cooler.

• A larger inner disk radius R_in leads to a larger solid angle $\mathrm{\Delta \Omega }\mathrm{\propto }{\mathit{R}}_{\mathrm{in}}^{\mathrm{2}}$, hence substantially stronger CO lines. If R_in is increased, the heating UV flux gets radially diluted ∝$\mathrm{1}\mathit{/}{\mathit{R}}_{\mathrm{in}}^{\mathrm{2}}$, but can still provide sufficient gas heating and CO excitation up to about R_in ~ 3 AU, depending on stellar UV luminosity. For even larger R_in, however, the excitation conditions for fundamental CO line emission break down, see Fig. 19.

  Fig. 19 Impact of inner disk radius Rin on CO υ = 1 → 0 line emission fluxes and FWHM. The black line shows the computed R(10) line fluxes for the reference model, when varying Rin between 0.07 AU and 100 AU. The blue line shows the R(10) /R(35) line ratio, and the red line shows the computed R(10) line widths.

• There is a strong impact of the dust size parameters on the CO fundamental line emission. This is due to changes in the dust UV optical depths τ_UV, see Fig. 20. Dust parameter choices which favour larger particles (larger a_min, larger a_max, smaller a_pow) and less amorphous carbon reduce the UV dust opacities, i.e. the UV light can penetrate deeper into the CO line emitting regions, which increases the gas temperatures and CO line fluxes.

• After the inner disk radius, the gas/dust ratio has the largest impact on the fundamental CO emission. The reason for this effect is the same as for the dust size parameters. Larger gas/dust ratios lead to an increase of the τ_line/τ_UV ratio.

To summarise, the impinging UV flux and the local τ_line/τ_UV ratio regulate the strength of the CO fundamental line emissions, but the disk shape determines how much of the stellar UV flux reaches a disk region under consideration. The reference model shows very weak CO fundamental line fluxes with broad, box-like profiles, which would in fact lead to non-detections, but the model parameters can be changed to result in observable line flux levels and more typical CO line profiles.

The strong dependence of the CO υ = 1 → 0 line fluxes on R_in suggests that any disk shape irregularities in the inner regions ≲30 AU, which create highly extended vertical gas columns that are exposed to fresh UV light from the central star, could easily dominate the CO fundamental emission lines. Such irregularities could be, for example, (i) large inner holes; (ii) disk gaps with highly puffed-up secondary inner walls; (iii) spiral waves which are warmer than the surrounding gas, causing those regions to stick out vertically; or (iv) the launching regions of disk winds.

6. Summary and conclusions

The analysis and interpretation of observational data from Class II and III protoplanetary disks is a challenging task. Various hydrodynamical, chemical, dust, and radiative processes are coupled to each other in complicated ways, and current disk modelling groups are using quite different assumptions to setup their 2D radiative transfer and thermo-chemical models concerning disk shape, dust size and opacity parameters, treatment of PAHs, and dust settling.

We have systematically investigated the effects of those assumptions in this paper by have studying the impact of the associated model parameters on the various continuum and line predictions, using a holistic disk modelling approach which allows us to calculate all continuum and line observations on the basis of a single model. The most important effects are:

• Disk shape matters. In particular, scale heights and disk flaringhave large impacts on the gas and dust temperature distribution inthe disk, and hence on SED, visibilities and gas emission lines. Forvery cold, self-shadowed disk configurations, ice formation canstrongly reduce the (sub-mm) emission lines.

• Hydrostatic T Tauri disk models produce only little near-IR excess, but a strong flaring of the outer disk, which leads to a  re-increase  of the spectral flux νF_ν between λ = (20−50) μm, whereas the opposite is typically observed for non-transitional disk.

• Dust size is important, not only for modelling the continuum, but also for modelling chemistry and emission lines, which is often not discussed in thermo-chemical disk models. By extending the dust size distribution to mm sizes, as required to fit the SED, the UV dust opacity is reduced by a factor of about 100, which means that UV photons can penetrate much deeper into the disk with ample effects on chemistry, temperature structure and gas emission lines.

• New dust standard opacities have been developed. Guided by a study of the multi-wavelength optical properties of dust aggregates particles (Min et al. 2016), we have developed a simplified and fast numerical treatment for dust opacities. We propose to use an effective porous mixture of amorphous laboratory silicates with amorphous carbon, a powerlaw size distribution with a distribution of hollow spheres to capture the most important size, material, and shape effects. A Fortran-90 package to compute the DIANA standard dust opacities¹¹.

• Disk flaring and/or dust settling? Dust settling affects primarily the outer disk regions, leading to cooler disks, lower continuum fluxes at mid-IR to cm wavelengths, and smaller apparent disk sizes in the millimetre continuum. These effects on the continuum observables are very similar to those in models with little or lacking disk flaring. However, concerning the gas emission lines, dust settling has just the opposite effect. Dust settling leaves the vertically extended gas bare and exposed to the stellar UV radiation, leading to higher gas temperatures and stronger gas emission lines in general. Thus, we can expect to distinguish between disk flaring and dust settling by observing in particular far-IR emission lines.

• PAHs can have important effects on the disk radiative transfer, even if the mid-IR PAH features are not visible in the SED. If the PAH abundance reaches about 10% of the interstellar standard, the PAH opacities start to catch up with the dust opacities in the UV and blue parts of the spectrum. Since the PAHs are not settled, and have negligible scattering opacities, the PAHs change the ways in which UV photons reach the disk. In fact, the PAHs can effectively “shield” the disk from UV photons. We are proposing a simplified method by treating the PAHs consistently in the continuum radiative transfer assuming radiative equilibrium, with a PAH temperature independent from the dust temperature.

With regard to particular observations, we find the following robust dependencies that can be used for diagnostic purposes.

• The SED mm-slope is more gentle than expected from the dustopacity-slope if the disk is cold in the midplane, which happensconsistently in all of our T Tauri models.

• The [OI] 63.2 μm line is optically thick under all explored circumstances, and probes the tenuous layers above the molecular disk at radii ~(10−100) AU. Since this line provides the most important cooling process in these layers, it is subject to a self-regulation mechanism where the line luminosity must equal the spatially integrated heating rate in that region, in form of various UV and X-ray processes.

• The high-J CO lines are optically thick and triggered by an excess of the gas temperature over the local dust temperature at the upper edge of the molecular layer around 1−10 AU. The lines are stronger for larger τ_line/τ_UV ratios, and hence an excellent tracer of the gas/dust ratio.

• The (sub-)mm ¹²CO lines probe the radial extension and conditions in the tapering-off distant gas around the disk, well beyond the radial zones responsible for the (sub-)mm continuum. These conclusions are robust, almost purely geometrical, because these lines are always optically thick in our models, and because CO is such an abundant and chemically robust constituent of the disks.

• The (sub-)mm CO isotopologue lines of ¹³CO and C¹⁸O require more sophisticated modelling, where the sharpness of the outer disk edge, the strength and vertical extension of the disk midplane shadow, and vertical temperature gradients all play a significant role. The ¹²CO/¹³CO isotopologue line ratio is quite independent from disk mass, because both lines are usually optically thick in our models, but increases significantly in cold disks where the increased efficiency of CO ice formation affetcs the ¹³CO lines more than the ¹²CO lines.

• The CO υ = 1−0 fundamental lines are massively optically thick in all our models, and hence provide a clear diagnostic for the existence of warm (≳500 K) gas in the inner disk regions, with a gas temperature in clear excess to the underlying dust temperature. In our standard T Tauri disk setup, these lines are quite weak, because these conditions are mostly met just at the inner rim. Stronger and more extended CO emission lines are obtained, however, if the impinging UV flux is larger, and/or if the local optical depth ratio τ_line/τ_UV is larger. In practise, such extended emission occurs in models where sub-micron grains are missing, and/or where the gas/dust ratio is larger, leading to more realistic line profiles. Disk shape irregularities, like inner holes with diameters of order several AU, can lead to much stronger CO fundamental line emissions as well.

Our results demonstrate that the various continuum and line observables probe the physical conditions at very different radii and different heights in the disk. Thus, only a combination of suitable multi-wavelength dust and gas observations can break the various degeneracies, for example those in SED modelling, and can lead to more reliable disk diagnosis.

This paper series aims at setting new disk modelling standards for the analysis of multi-wavelength continuum and line observations for protoplanetary disks, with easy to implement, yet physically grounded, and practical assumptions, which are sufficiently motivated by observations. We will continue this series by exploring the effects of chemical networks and rates in Paper II (Kamp et al., in prep.) and element abundances in Paper III (Rab et al. 2015).

We intend to offer our modelling tools and collected data sets to the community at the end of the FP7 DIANA project¹².

---

Acknowledgments

We would like to thank Dr. Pieter Degroote for compiling a database of updated UV to far-IR photometric measurements for all our 80 targets, including accurate filter functions and zero-points, and Dr. Hein Bertelsen for improving the discussion about fundamental CO emission. The research leading to these results has received funding from the European Union Seventh Framework Programme FP7-2011 under grant agreement no 284405. C. Rab and C. Baldovin-Saavedra acknowledge funding by the Austrian Science Fund (FWF), project number P24790, and the Austrian Research Promotion Agency (FFG) under grant agreement FA 538022, respectively. A. Carmona acknowledge support by the Momentum grant of the MTA CSFK Lendület Disk Research Group. The computer simulations were carried out on the UK MHD Consortium parallel computer at the University of St Andrews, funded jointly by STFC and SRIF.

References

Appendix A: Stellar parameters and irradiation

Appendix A.1: Stellar parameters

The photospheric component of the stellar emission is characterised by the effective temperature T_⋆, the surface gravity g, and the stellar (photospheric) luminosity L_⋆. Since these properties are fairly well-known for most of our target objects in the literature, we have decided not to put too much efforts into this issue. Assuming solar abundances for the star, we use standardPHOENIX stellar atmosphere models (Brott & Hauschildt 2005) to fit T_eff, L_⋆ and the interstellar extinction A_V, to our photometric and spectroscopic data (see Fig. A.1). A thorough determination of T_⋆ requires to fit high-resolution optical spectra. For most target objects, this has been done already in the literature, so we can pick T_⋆ from the literature and only fit A_V and L_⋆ for the assumed distance d, using at first an estimate of the surface gravity g.

Once T_⋆ and L_⋆ are determined, we involve pre-main sequence stellar evolutionary models (Siess et al. 2000) to find the stellar mass M_⋆ and the age. The stellar radius and surface gravity are then given by ${\mathit{L}}_{\mathit{\star }}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{R}}_{\mathit{\star }}^{\mathrm{2}}\mathit{\sigma }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}$ and $\mathit{g}\mathrm{=}\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathit{/}{\mathit{R}}_{\mathit{\star }}^{\mathrm{2}}$. We use the resulting value for g to redo the fitting above. This procedure is found to converge very quickly, and – thanks to having fixed T_⋆ from the literature – gives quite unambiguous results, see examples in Table A.1.

Appendix A.2: Stellar UV irradiation

The UV irradiation by the central star is much more difficult to determine, and hampered by the lack of high-quality UV data, especially in the hard FUV region (<130 nm), and for Lyα. We have systematically scanned and collected UV data from IUE, FUSE, HST/STIS, HST/COS, and HST/ACS. These data have been re-binned and collated, using the inverse square of the flux uncertainties as weighting factors, following the idea used by Valenti et al. (2000, 2003). The results of this data collection will be described elsewhere (Dionatos et al. 2015, in prep.).

For many objects, the UV data is poor and incomplete, and we have to use template stars or other tools to complete it. One option is to use a simple powerlaw I_λ ∝ λ^p_UV, or ${\mathit{I}}_{\mathit{\nu }}\mathrm{\propto }{\mathit{\lambda }}^{{\mathit{p}}_{\mathrm{UV}}\mathrm{+}\mathrm{2}}\mathit{,}$(A.1)where the proportionality constant and the UV powerlaw index p_UV can be roughly fitted to the existing data, or in order to fill in gaps in the data. For TW Hya and AB Aur, however, the UV data quality is excellent, and we can directly use the data to determine the UV irradiation, see Fig. A.1 and Table A.1.

In case of the Herbig Ae/Be stars, the UV data seems mostly photospheric in character, still having absorption lines down to wavelengths of about 150 nm and below, but at even shorter wavelengths, the spectrum changes character, is dominated by emission lines and is in excess to photospheric models. The soft part of the UV data can then be used to improve the fits of the stellar parameters and A_V. For AB Aur, we do not need to pick T_⋆ from the literature. We find good agreement between the UV data and our photospheric model with T_⋆ = 9550 K down to about 130 nm, in excellent agreement with (DeWarf et al. 2003).

Appendix A.3: Stellar X-ray irradiation

X-ray data were obtained by using archival and new data from the two X-ray observatories XMM-Newton (Jansen et al. 2001) and Chandra (Weisskopf et al. 2000). The spectra were extracted using the software of either science centres (SAS and CIAO). To get the source spectra, we select a circular extraction region around the center of the emission, while the background area contained a large source-free area on the same CCD. The extraction tools (EVSELECT for XMM and SPECEXTRACT for Chandra) delivered the source and background spectra as well as the redistribution matrix and the ancilliary response files.

The extracted spectra were then fitted with the X-ray emission model XSPEC (Arnaud 1996) assuming a collisionally ionised plasma (VAPEC) and a model for an absorption column (WABS). The element abundance values in the VAPEC models were set to typical values for pre-main sequence stars, as chosen by the XEST project (Güdel et al. 2007). Either a one component (1T), a two component (2T) or a three component (3T) emission model is fitted to the data. Highly absorbed sources or scarce data allow only for 1T fits. The fit delivers an absorption column density towards the source N_H, and a plasma emission temperature T_X and an emission measure EM for each component. Finally, the unabsorbed spectrum is calculated after setting the absorption column density parameter to zero, and the flux is derived by integrating over the energy range 0.3−10 keV.

In cases where no detailed X-ray data are available, a more simple two-parameter approach is used considering a free-free (bremsstrahlung) continuum with total luminosity L_X and a fitted X-ray emission temperature T_X,fit as ${\mathit{I}}_{\mathit{\nu }}\mathrm{\propto }\frac{\mathrm{1}}{\mathit{\nu }}\exp \left(\mathrm{-}\frac{\mathit{h\nu }}{\mathit{k}{\mathit{T}}_{\mathit{X,}\mathrm{fit}}}\right)\mathrm{·}$(A.2)Table A.1 shows that T_X,fit results to be close to the highest X-ray component temperatures found, whereas the mean temperature ⟨ T_X ⟩ (linear mean of component temperatures, weighted by component emission masses) would result in a very bad fit and should  not  be used in Eq. (A.2). The unabsorbed X-ray emission spectrum is finally converted to units of surface intensities I_ν(r = R_⋆) [erg/cm²/ s/Hz/sr], using the previously determined values for d and R_⋆, and merged with the UV data and photospheric model spectrum.

Appendix A.4: Background radiation

Protoplanetary disks are irradiated not only by the central star, but also from the environment. We assume an isotropic interstellar (IS) background radiation field ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{ISM}}$ with 3 components, the IS UV-field (created by distant O-stars), the cosmic background (CMB, approximated by a 2.7 K Planckian), and an infrared background radiation field ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{IR}}$ (created by distant stars and molecular clouds).

The interstellar UV field is approximated by a 20000 K-black-body with a tiny dilution factor W_dil = 9.85357 × 10^-17 such that χ^ISM = 1 corresponds to the integrated standard IS UV field of Draine & Bertoldi (1996)${\mathit{I}}_{\mathit{\nu }}^{\mathrm{ISM}}\mathrm{=}{\mathit{\chi }}^{\mathrm{ISM}}\mathrm{1.71}{\mathit{W}}_{\mathrm{dil}}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathrm{20000}\hspace{0.17em}\mathrm{K}\mathrm{\right)}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{\chi }}_{\mathrm{IR}}^{\mathrm{ISM}}{\mathit{I}}_{\mathit{\nu }}^{\mathrm{IR}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathrm{2.7}\hspace{0.17em}\mathrm{K}\mathrm{\right)}\mathit{.}$(A.3)For the infrared background field ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{IR}}$, we take the spectral shape from (Mathis et al. 1983, their Tables A.3 and B.1), and adjust the parameter ${\mathit{\chi }}_{\mathrm{IR}}^{\mathrm{ISM}}$ such that the integrated background intensity equals $\mathrm{\int }{\mathit{I}}_{\mathit{\nu }}^{\mathrm{ISM}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{=}\frac{\mathit{\sigma }}{\mathit{\pi }}{\mathit{T}}_{\mathrm{back}}^{\mathrm{4}}\mathit{.}$(A.4)A blackbody would accomodate its temperature to T_back is this background radiation field. T_back limits the dust temperatures in the disk, because when including the star as additional light source, the temperatures can only increase. Without the infrared contribution $\mathrm{(}{\mathit{\chi }}_{\mathrm{IR}}^{\mathrm{ISM}}\mathrm{=}\mathrm{0}\mathrm{)}$ the background temperature results to be 2.97 K, but we can increase ${\mathit{\chi }}_{\mathrm{IR}}^{\mathrm{ISM}}$ to simulate a disk in close proximity to, e.g., a star formation region which provides additional IR background radiation, with expected impact on the mm-slope (see Sect. 5.3.3). For example, to achieve T_back = 5 K and 10 K, values for ${\mathit{\chi }}_{\mathrm{IR}}^{\mathrm{ISM}}\mathrm{=}\mathrm{5.1}$ and 92 are required, respectively. The original work of Mathis et al. (1983, i.e. ${\mathit{\chi }}_{\mathrm{IR}}^{\mathrm{ISM}}\mathrm{=}\mathrm{1}$ is derived from stars and molecular clouds at a distance at 10 kpc from the Galactic centre.

When calculating spectral fluxes from a model with background radiation, we need to subtract the background as ${\mathit{F}}_{\mathit{\nu }}\mathrm{=}\mathrm{\int }\mathrm{(}{\mathit{I}}_{\mathit{\nu }}\mathrm{-}{{\mathit{I}}_{\mathit{\nu }}^{\mathrm{ISM}}}^{\mathrm{)}}\hspace{0.17em}\mathrm{d\Omega }\mathit{.}$(A.5)The background subtraction (Eq. (A.5)) is important in particular at long wavelengths where the CMB is bright. For example, at λ = 1 cm, the CMB multiplied by solid angle $\mathrm{\Omega }\mathrm{=}\mathit{\pi }{\mathit{R}}_{\mathrm{out}}^{\mathrm{2}}\cos \mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{/}{\mathit{d}}^{\mathrm{2}}$ is 15 × stronger than the disk signal from the reference model. Equation (A.5) is also necessary to make the results independent of considered image size (as long as the disk is well contained in the image). We note, however, that Eq. (A.5) can lead to negative fluxes, for example an edge-on disk in the UV or whenever the disk appears darker than the background at the considered wavelength, for example the “silhouette disks” in the Orion nebula.

Appendix B: The PAH temperature distribution

Figure A.2 shows some results for the T Tauri reference model introduced in Sect. 5, but, in order to demonstrate the effects, with an artificially large PAH abundance of f_PAH = 10. The dust in the outer disk parts is strongly concentrated towards the midplane, whereas the PAH molecules stay co-spatial with the gas by assumption. In all optically thin regions, and in the surface of the inner rim, the PAH temperatures are higher by a large factor ~ 1.5−2 as compared to the dust temperatures, simply due to the very blue absorption characteristic of the PAHs which facilitates radiative heating. In contrast, in the optically thick midplane of the disk, we find T_PAH ≈ T_dust as expected. The PAH emission features result from the large temperature contrast between PAHs and dust in the optically thin upper and inner disk regions <(1−10) AU, depending on wavelength, see Fig. A.3. In the outer optically thin disk regions, the PAH temperature re-increases due to the interstellar UV irradiation.

The detailed physical modelling of the PAH molecules with the Monte Carlo codes allows for a treatment of the quantum heating by single photons absorption events by the PAHs, followed by radiative cooling, which leads to a stochastic PAH temperature distribution p(T) at every point in the model. In the following, we can thereby check the validity of our simplified treatment of the PAHs in radiative equilibrium, as outlined in Sect. 3.8. The PAH emissivity [erg/s/Hz/sr /PAH-molecule] is given by ${\mathit{ϵ}}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{\int }{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}\mathit{p}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{T,}$(B.1)where ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{PAH}}\hspace{0.17em}\mathrm{[}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\mathit{/}$PAH-molecule] is the PAH absorption cross section and B_ν(T) is the Planck function. The temperature distribution function p(T) is normalised to ∫p(T) dT = 1.

In order to quantify the mean PAH temperature and the width of the PAH temperature distribution, we consider $\begin{array}{ccc}& & \int {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{PAH}}{\mathit{B}}_{\mathit{\nu }}\hspace{0.17em}(\mathrm{⟨}{\mathit{T}}_{\mathrm{PAH}}\mathrm{⟩})\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{=}\mathrm{\int }{\mathit{ϵ}}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{d}\mathit{\nu }\\ & & {\mathit{\kappa }}_{\mathrm{3.3}\hspace{0.17em}\mathit{\mu }\mathrm{m}}^{\mathrm{PAH}}{\mathit{B}}_{\mathrm{3.3}\hspace{0.17em}\mathit{\mu }\mathrm{m}}\mathrm{\left(}{\mathit{T}}_{\mathrm{high}}\mathrm{\right)}\mathrm{=}{\mathit{ϵ}}_{\mathrm{3.3}\hspace{0.17em}\mathit{\mu }\mathrm{m}}^{\mathrm{PAH}}\\ & & {\mathit{\kappa }}_{\mathrm{30}\hspace{0.17em}\mathit{\mu }\mathrm{m}}^{\mathrm{PAH}}{\mathit{B}}_{\mathrm{30}\hspace{0.17em}\mathit{\mu }\mathrm{m}}\mathrm{\left(}{\mathit{T}}_{\mathrm{low}}\mathrm{\right)}\mathrm{=}{\mathit{ϵ}}_{\mathrm{30}\hspace{0.17em}\mathit{\mu }\mathrm{m}}^{\mathrm{PAH}}\\ & & {\mathit{\sigma }}_{\mathit{T}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{T}}_{\mathrm{high}}\mathrm{-}{\mathit{T}}_{\mathrm{low}}\mathrm{)}\mathit{,}\end{array}$i.e. we define a high temperature to match the 3.3 μm PAH emission and a low temperature to match the continuous 30 μ PAH emission. The mean temperature ⟨ T_PAH ⟩ is related to the total frequency integrated PAH emission. Usually, we find T_low> ⟨ T_PAH ⟩ >T_high in the models, see an example in Fig. A.3, lower left part. In the example shown, the relative width of the PAH temperature distribution function is σ_T/ ⟨ T_PAH ⟩ ≈ 15%.

Figure A.3 shows that σ_T/ ⟨ T_PAH ⟩ anti-correlates with ⟨ T_PAH ⟩. Large average PAH temperatures imply a sharply peaked temperature distribution function. The upper right plot shows which spatial disk regions are responsible for the various PAH emission features. The PAH 3.3 μm feature originates in the innermost 1 AU disk regions, whereas the 11.3 μm emission region stretches out to about 10 AU. In all cases, the PAH emissions come from borderline optically thin altitudes where the dust is not yet optically thick, and where ⟨ T_PAH ⟩ ≫ T_dust. All these warm regions are characterised by quite a sharply peaked temperature distribution function σ_T/ ⟨ T_PAH ⟩ ≈ 0.2%−30%, which explains why the fast approximate method, assuming the PAHs to be in radiative equilibrium, works so well for the SED modelling.

Figure B.1 shows, however, that there are substantial differences at larger radii, which are not important for the integral emission of the PAH features. The blue model in Fig. B.1 (assuming equal dust and PAH temperatures) fails completely to predict the PAH emission features in the SED. The black model (PAHs in radiative equilibrium) is good for the SED and sufficient to predict the intensity profiles up to ~20 AU, but only the full stochastic treatment of the PAHs (red model) can predict the intensity profile beyond ~20 AU.

Appendix C: Time-dependent chemistry

To investigate the effects of chemical evolution in the disk on the observable gas emission lines, we have computed additional models where the disk structure, the dust size distribution, dust settling and opacities, the dust and gas temperatures, and the internal radiation field are taken over from the reference model, but at each grid point, the chemical rate network is advanced in time from zero to 5 Myr, starting from initial concentrations typical for the dark cores of molecular clouds, see (Helling et al. 2014) for details.

Figure C.1 shows the resulting line fluxes as function of disk age. After an initial relaxation phase which lasts a few 10⁵ yrs, the gas emission line fluxes become constant in the model, and do not change significantly afterwards. This behaviour is a consequence of the relatively short chemical relaxation timescales in most line forming regions situated well above the icy midplane, compare Fig. 12. The chemical relaxation timescale τ_chem(r,z) is introduced and discussed in (Woitke et al. 2009, see Sect. 8.3 and Fig. 13 therein). The relaxation of the optical and IR lines is indeed very quick; we measure mean values of the chemical relaxation timescale as 1.2 yrs, 90 yrs, 120 yrs, 4000 yrs and 8000 yrs in the line forming regions of [OI] 6300 Å, o-H₂ 17 μm, CO J = 18−17, CO J = 10−9 and [OI] 63 μm, respectively.

However, the CO mm-line fluxes do not fully converge to the results obtained by the time-independent reference model, see arrows on the r.h.s. of Fig. C.1, even for integration times as long as 10 Myr. This is consistent with the very long chemical relaxation timescales ⟨ τ_chem ⟩ we find for those lines, namely 20 Myr, 50 Myr and 100 Myr in the line forming regions of CO J = 2−1, ¹³CO J = 2−1, and C¹⁸O J = 2−1, respectively.

The difference of the CO concentration between the t = 1 Myr disk model and the time-independent reference model is shown in Fig. C.2. All upper disk regions have practically indistinguishable CO concentrations, however, in the dark icy midplane (vertical visual extinction A_V> 10), CO cannot freeze out directly where the dust temperature exceeds about 20 K (at radii ≲30 AU in this model). Instead, the CO is very slowly dissociated in these regions by reactions with He⁺ created by cosmic rays. Most of the liberated oxygen atoms re-form CO, but a small part can form other molecules with higher adsorptions energies like water, and these will freeze out immediately. This mechanism provides a slowly ticking chemical clock in the dark midplane regions of protoplanetary disks, leading to gaseous carbon-to-oxygen ratios C/O ≈ 1 on Myr timescales between the water and the CO ice-lines, and to C/O ≫ 1 on even longer timescales, see (Öberg et al. 2011; Helling et al. 2014).

The (sub-)mm CO isotopologue line fluxes are somewhat affected by these differences, for example remainders of cold gaseous CO in the dark icy midplane can partly re-absorb the CO line photons emitted by the warm surface of the disk on the opposite side which cross the midplane. In this way, the CO isotopologue line fluxes of the time-independent model can be up to 30% larger than those of the time-dependent models. However, these effects are small compared to the impact of disk mass, size and shape, dust size distribution and settling, and the dust/gas ratio, and negligible for most optical and near-IR to far-IR lines.

We conclude that the emission lines of chemically robust gas tracers, such as CO, H₂ and the oxygen atom, are little affected by time-dependent chemical effects in our models. Disk evolution will rather have an impact on those lines via the changing stellar parameters, the changing shape of the disk, and physical gas and dust evolution. See e.g. Akimkin et al. (2013) for time-dependent disk models which include those effects.

Appendix D: Properties derived from the model

The mean gas temperature in the disk ⟨ T_gas ⟩, the mean dust temperature ⟨ T_dust ⟩, the near-IR excess, the 10 μm SED amplitude, the dust absorption mm and cm-slopes β_abs, and the millimetre and centimetre flux-slopes α_SED, are calculated as $\begin{array}{ccc}& & \\ & & \\ & & \mathrm{near}\mathrm{-}\mathrm{IRexcess}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathrm{d}}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{2}\hspace{0.17em}\mathit{\mu }\mathrm{m}}^{\mathrm{7}\hspace{0.17em}\mathit{\mu }\mathrm{m}}\mathrm{(}{\mathit{F}}_{\mathit{\lambda }}\mathrm{-}{\mathit{F}}_{\mathit{\lambda }}^{\mathit{\star }}\mathrm{)}\mathrm{d}\mathit{\lambda }\\ & & \mathrm{10}\hspace{0.17em}\mathit{\mu }\mathrm{mampl}\mathrm{.}\mathrm{=}\frac{{\mathit{F}}_{\mathit{\nu }}\mathrm{\left(}\mathrm{9.6}\hspace{0.17em}\mathit{\mu }\mathrm{m}\mathrm{\right)}}{\sqrt{{\mathit{F}}_{\mathit{\nu }}\mathrm{\left(}\mathrm{6.8}\hspace{0.17em}\mathit{\mu }\mathrm{m}\mathrm{\right)}\mathrm{·}{\mathit{F}}_{\mathit{\nu }}\mathrm{\left(}\mathrm{13.1}\hspace{0.17em}\mathit{\mu }\mathrm{m}\mathrm{\right)}}}\\ & & {\mathit{\beta }}_{\mathrm{abs}}^{\mathrm{mm}}\mathrm{=}\mathrm{-}\frac{\log {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{\left(}\mathrm{0.85}\hspace{0.17em}\mathrm{mm}\mathrm{\right)}\mathrm{-}\log {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{\left(}\mathrm{1.3}\hspace{0.17em}\mathrm{mm}\mathrm{\right)}}{\log \mathrm{0.85}\hspace{0.17em}\mathrm{mm}\mathrm{-}\log \mathrm{1.3}\hspace{0.17em}\mathrm{mm}}\\ & & {\mathit{\beta }}_{\mathrm{abs}}^{\mathrm{cm}}\mathrm{=}\mathrm{-}\frac{\log {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{\left(}\mathrm{5}\hspace{0.17em}\mathrm{mm}\mathrm{\right)}\mathrm{-}\log {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{\left(}\mathrm{10}\hspace{0.17em}\mathrm{mm}\mathrm{\right)}}{\log \mathrm{5}\hspace{0.17em}\mathrm{mm}\mathrm{-}\log \mathrm{10}\hspace{0.17em}\mathrm{mm}}\\ & & \\ & & \end{array}$where ρ_gas is the gas mass density, ρ_dust = ρ_gas·δ the dust mass density, δ the local dust-to-gas mass ratio, dV = 2πr dr dz the volume element, d the distance, ${\mathit{F}}_{\mathit{\lambda }}\mathrm{=}\frac{\mathit{\nu }}{\mathit{\lambda }}{\mathit{F}}_{\mathit{\nu }}$ the spectral flux per wavelength interval, ${\mathit{F}}_{\mathit{\lambda }}^{\mathit{\star }}$ the flux from the naked star and ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}$ the dust absorption opacity [cm²/ g(dust)].

Appendix E: Fluxes of optically thick emission lines

It is noteworthy that all observable emission lines discussed in this paper are optically thick in the reference model, with the only exception being C¹⁸O J = 2 → 1 which is borderline optically thin. For large optical depths we can use the Eddington-Barbier approximation ${\mathit{I}}_{\mathit{\nu }}\mathrm{\approx }{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{(}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{\approx }\{\begin{array}{c}\\ {\mathit{S}}_{\mathit{\nu }}^{\mathrm{cont}}\hspace{0.17em}\mathrm{(}{\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{cont}}\mathrm{=}\mathrm{1}\mathrm{)}& \mathrm{,}\mathrm{if}\mathrm{|}\mathit{\nu }\mathrm{-}{\mathit{\nu }}_{\mathrm{0}}\mathrm{|}\mathrm{\gg }\mathrm{\Delta }\mathit{\nu }\\ {\mathit{S}}_{\mathit{\nu }}^{\mathrm{line}}\hspace{0.17em}\mathrm{(}{\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{line}}\mathrm{=}\mathrm{1}\mathrm{)}& \mathrm{,}\mathrm{if}\mathrm{|}\mathit{\nu }\mathrm{-}{\mathit{\nu }}_{\mathrm{0}}\mathrm{|}\mathrm{\le }\mathrm{\Delta }\mathit{\nu ,}\end{array}$(E.1)where S_ν is the general source function, ${\mathit{S}}_{\mathit{\nu }}^{\mathrm{cont}}$ the continuum source function and ${\mathit{S}}_{\mathit{\nu }}^{\mathrm{line}}$ the line source function. ${\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{cont}}$ and ${\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{line}}$ are the continuum and line centre optical depth, ν₀ and Δν = FWHM/ 2 are the line centre frequency and observed frequency width, respectively. If the continuum is optically thin ${\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{cont}}\mathit{<}\mathrm{1}$, its contribution can be neglected in Eq. (E.1).

Integration over solid angle and frequency, and continuum subtraction, results in the line flux $\begin{array}{ccc}{\mathit{F}}_{\mathrm{line}}& \mathrm{=}& \int \mathrm{\int }\mathrm{(}{\mathit{I}}_{\mathit{\nu }}\mathrm{-}{{\mathit{I}}_{\mathit{\nu }}^{\mathrm{cont}}}^{\mathrm{)}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\hspace{0.17em}\mathrm{d\Omega }\\ & \mathrm{\approx }& \\ & \mathrm{\approx }& \end{array}$where ΔΩ is the solid angle occupied by the part of the disk that is optically thick in the line, and 2Δν = FWHM is the observed frequency full width half maximum. The first approximation (Eq. (E.2)) is valid only if the excitation conditions (${\mathit{S}}_{\mathit{\nu }}^{\mathrm{line}}$, ${\mathit{S}}_{\mathit{\nu }}^{\mathrm{cont}}$) are about constant in the line forming region, otherwise the integration over the solid angle cannot be carried out this way. The second approximation (Eq. (E.3)) is valid for LTE only. Figure 20 shows the typical situation encountered in disks where we look through a warm gas toward the cold, optically thick midplane.

Appendix F: Numerical convergence

Figure F.1 shows three series of ProDiMo models with increasing spatial resolution in the underlying numerical grid. The resulting continuum predictions like near-IR excess, 10 μm amplitude, millimetre and centimetre slopes etc. (see Sect. D) are robust. Even quick 30 × 30 models are sufficient to predict the SED and derived quantities, with an accuracy better than 5% with respect to the results from the big 160 × 150 reference model.

The grid resolution is more critical, however, when studying emission lines. A too coarse spatial grid usually leads to an over-prediction of the emission line fluxes. Most critical are lines which originate in a small portion of the disk volume, like the weak o-H₂ and high-J CO lines, but also [OI] 63.2 μm and [OI] 6300 Å. Here, the radial grid resolution is more important than the vertical grid resolution. The only counter-example to this rule are the ro-vibrational CO lines which are mostly emitted directly from the surface of the inner rim in this model. Here, the vertical grid resolution is more important. In summary, we need a grid resolution of about 100 × 100 to achieve an accuracy better than 10% for all predictions, see grey shaded box in Fig. F.1.

Figure F.2 shows the numerical convergence of MC→ProDiMo chain models, when varying the spatial grid resolution in the Monte-Carlo (MC) programs. It is reassuring to see that the MC→ProDiMo chain models actually produce results that are very similar as compared to the pure ProDiMo models. For sufficient spatial resolution in the MC models (again about 100 × 100 grid points), the deviations in continuum results are smaller than 10%, and line results agree better than 15%, where most critical are the faint mid-far IR o-H₂ and high-J CO lines. Other, e.g. (sub-)mm line results are more robust. There is also no obvious asymmetry in Fig. F.2, e.g. some lines are weaker whereas others are stronger when using the MC→ProDiMo chain models.

For sufficient spatial resolution in the MC models, the MCFOST and MCMax results show a similar pattern with respect to the pure ProDiMo reference results, i.e. the deviations between MCFOST→ProDiMo and MCMax→ProDiMo models are actually smaller than the deviations between those models and pure ProDiMo models.

Appendix G: Impacts of additional model parameters

Figure G.1 shows the impact of various model parameters on the mean CO fundamental line emission strengths and profiles. The results are discussed in Sect. 5.3.7.

Figure G.2 compares the results obtained from a model using anisotropic scattering to the reference model. The largest difference concerns the amplitude of the 10 μm silicate emission feature, the anisotropic model has a stronger amplitude by about 14%, making the feature clearly more visible in the SED plot as compared to the reference model.

Concerning the gas emission lines, we see only little effects, and no clear trend. The disk is mostly UV illuminated from the top, which requires at least one scattering event in the borderline optically thin disk surface layers. Since an

angle-dependent treatment of dust scattering favours forward scattering, one would expect the UV illumination of the disk from above to be reduced, because ≈90° scattering events seem to be required. However, this does not seem to be entirely true. The question whether a scattered UV photon reaches the disk or not is simply determined by whether the photon is scattered upwards or downwards, which is a 50%/50% chance even for anisotropic scattering. With anisotropic scattering, we rather redistribute the entry points where the scattered photons enter the disk, favouring the outer disk regions, and their initial entry angle. Multiple scattering also reduces the effects.

We measure line flux and FWHM differences of order 10%, with no clear trend, the only exception being CO J = 10 → 9 with is actually enhanced by 28%. These results are close to the “noise level” expected from various numerical effects in the models, compare Figs. F.1 and F.2.

All Tables

All Figures