© ESO, 2011

1. Introduction

Gas-rich dust disks around young stars (hereafter, protoplanetary disks) provide the raw material to build up new planets. The physical, thermal, and chemical conditions in the disk, the timescale over which the gas disperses, and the physical mechanisms contributing to the gas dispersal are keys to understanding what type of planets can form and on what timescales.

Significant progress has been made in the past few years in measuring the dispersal timescale of the dust component of the disk. Infrared surveys of nearby star-forming regions and associations have established that the frequency of optically thick dust disks decreases exponentially with time (e.g. Mamajek 2009). By an age of 10 Myr only a few percent of pre-main sequence Sun-like stars (T Tauri stars) still retain an optically thick dust disk (see e.g. Hernández et al. 2008; Pascucci & Tachibana 2010, for reviews). Since the near-mid infrared excess (λ ≲ 30 μm) is sensitive to the presence of small dust grains, not larger than a few microns in size, these observations effectively trace the dispersal of small grains within about 10 AU from T Tauri stars. Millimetre observations, tracing colder dust at hundreds of AU from the central star, indicate a similarly fast clearing for the outer disk, within about 10−30 Myr for T Tauri stars (Carpenter et al. 2005). There is growing observational evidence that the dust disk lifetime depends on stellar mass. Disks around intermediate-mass (≳1.5   M_⊙) stars disperse in less than 10 Myr, whereas disks around low-mass stars (M dwarfs and brown dwarfs) persist for longer times (Carpenter et al. 2006; Currie et al. 2007; Riaz & Gizis 2008).

Due to observational challenges in detecting gas lines from disks and difficulties in interpreting them, much less is known about the evolution of the gas component of the disks. Three observables point to a dispersal timescale similar to (or possibly shorter than) the dust dispersal timescale: the exponential decrease with time in the frequency of accreting stars (Fedele et al. 2010); the non-detections of infrared gas lines from abundant molecules and atoms in tenuous dust disks (Hollenbach et al. 2005; Pascucci et al. 2006); upper limits on the H₂/dust mass ratio of less than 10 in two  ~12 Myr old edge-on disks (Lecavelier des Etangs et al. 2001; Roberge et al. 2005).

The aim of the Herschel open time key program “Gas in Protoplanetary Systems” (GASPS, Dent et al., in prep.) is to provide new insights into the chemical and gas temperature structure of protoplanetary disks, the gas/dust ratio, the gas dispersal timescale, and disk evolution. GASPS will acquire a large sample of sensitive far-infrared Herschel/PACS spectra for 240 disks in nearby star-forming regions and associations that span the critical 1−30 Myr age range over which disks are known to disperse. The primary signatures of the gas in the disk are expected to be the forbidden [OI] 63.2 μm, [OI] 145.5 μm, and [CII] 157.7 μm lines, as well as some CO and H₂O lines. The first GASPS papers have shown that (1) the [OI] 63 μm line can be used as primary gas indicator and is often detected toward protoplanetary disks (Mathews et al. 2010), (2) a combination of far-IR and (sub-)millimetre gas lines provides a promising tool to estimate the total gas mass of protoplanetary disks (Pinte et al. 2010), and (3) detailed models of individual sources allow for a better characterisation of the disk structure and shape, and the dust and gas components of protoplanetary disks (Meeus et al. 2010; Thi et al. 2010).

In this paper, we present an analysis of the circumstellar disk of ET Cha, an approximately 8 Myr old late-type T Tauri star, with the goals of characterising in detail its dust and gas content. ET Cha is one of the few nearby relatively old stars still possessing an optically thick dusk disk (Sicilia-Aguilar et al. 2009) and still accreting disk gas (Lawson et al. 2004). TW Hya and PDS 66 are two other well-known old stars with properties similar to ET Cha. Both disks have been studied in detail and show evidence of evolution with respect to 1−2 Myr old T Tauri disks in Taurus, for example, depleted inner disk in TW Hya (Calvet et al. 2002) and flatter disk structure for PDS 66 (Cortes et al. 2009). Both disks are likely to have too low disk masses to form giant planets at this evolutionary stage. ET Cha would be the third such old disk system where observational data allows for an in-depth-study of its dust and disk properties.

The paper is structured as follows. Section 2 provides an overview of the prior knowledge of the source. We then present new multi-wavelength observations of ET Cha in Sect. 3. Section 4 presents a detailed dust and gas disk model for ET Cha. We describe the main results of our models in Sect. 5. Finally, we discuss some critical aspects of the modelling, and the implications of both models, in Sect. 6, before we finish the paper with our conclusions in Sect. 7.

2. ET Cha: an old T Tauri star with active accretion

ET Cha (2 MASS J08431857-7905181, ECHA J0843.3-7905,also referred to as RECX 15¹) is a low-mass T Tauri star that was identified by (Lawson et al. 2002) as a member of the nearby, 8 Myr old η Chamaeleontis moving group (Mamajek et al. 2000)². The association is located only 97 pc away from the Sun (Mamajek et al. 1999) and is virtually unaffected by extinction (Luhman & Steeghs 2004), an ideal set of conditions to study circumstellar disks in detail. A slighly lower distance to ET Cha of 94.3 pc was reported by van Leeuwen (2007), but we have used the earlier and better known value of 97 pc for the modelling in this paper. Brandeker et al. (2006) obtained high-angular resolution images of ET Cha and concluded that it has no companions (brown dwarfs) outside of 10 AU (30 AU).

ET Cha is one of the few association members that possess a circumstellar disk, as indicated by a series of Spitzer observations that revealed the presence of dust including strong mid-infrared silicate features (Megeath et al. 2005; Bouwman et al. 2006; Gautier et al. 2008; Sicilia-Aguilar et al. 2009). Despite the age of the system, the infrared colours of the source are reminiscent of those of much younger (1–2 Myr) circumstellar disks. Furthermore, optical spectroscopy of ET Cha has shown that it is undoubtedly accreting, with a very strong and broad H_α emission line (Lawson et al. 2002; Lyo et al. 2004; Luhman & Steeghs 2004). Based on the observed H_α line, the mass accretion rate of ET Cha has been estimated to be 10^-9   M_⊙/yr (Lawson et al. 2004), and the disk inclination (by modelling the H_α line profile) to be about 60° as measured from face-on. The spectra also reveal a series of forbidden optical emission lines ([OI], [SII], [NII]) that unambiguously indicate the presence of a jet/outflow. In addition, the stellar absorption lines in these spectra allowed for precise spectral typing of the central star; all estimates agree with a M3−M3.5 spectral classification.

Among all members of the η Cha association, ET Cha shows the strongest variability in the visible of order 0.3−0.4 mag, which places it at the high end of the variability distribution of WTTS (see Fig. 1 in Grankin et al. 2008). The main feature in the observed lightcurve is a  ~12 day period which has been found in two consecutive years. Lawson et al. (2002) also noted a flare lasting for about 1.7 days. The more regular 12 day variations are most easily interpreted in terms of an accretion hotspot co-rotating with the star. However, our re-analysis of optical absorption line profiles (see Sect. 3.6) suggests that the stellar rotational period is much shorter, around 2 days. Therefore, the physical origin of the  ~12 day period remains uncertain. The stellar variability casts some doubt on the derivation of stellar parameters, because the photometry reported by Lawson et al. (2002) was taken “near maximum light”, when e.g. an accretion hot spot could contribute significantly to the observed flux. Furthermore, it prevents using optical and near-infrared colours in case of non-simultaneous data. For this reason, we have obtained new, simultaneous photometry, which is presented in Sect. 3.2.

Ramsay Howat & Greaves (2007) detected a clear gas signature from the disk of ET Cha in form of the H₂ v = 1 → 0 S(1) line in emission at 2.122 μm with Gemini/Phoenix, with an integrated line flux of (2.5 ± 0.1) × 10^-18 W/m². These high angular resolution observations showed a narrow line (FWHM = 18 ± 1.2 km s^-1) centred on the stellar velocity to within  ~1 km s^-1. No angular offset between the line and the star was detected at the level of 4 AU. Therefore, Ramsay Howat & Greaves (2007) argue that this line is emitted by H₂ gas in Keplerian rotation at  ~2 AU. Ramsay Howat & Greaves (2007) observed 3 other disk-bearing members of the η Cha association, but ET Cha was the only one with detectable H₂ emission. Bary et al. (2008) and Martin-Zaïdi et al. (2009, 2010) reported on several detections of H₂-lines toward other sources where the emission is also likely originating from the disk rather than from an outflow.

3. Observations and data reduction

3.1. Herschel/PACS

Herschel/PACS observations were obtained for ET Cha during the science demonstration phase. The photometric observations were obtained in the scan map mode in the blue (70 μm) and the red (160 μm) filters. Two different scan map angles were used, 45° (obsid 1342187338, 133 s) and 63° (obsid 1342189366, 220 s). Both scans were obtained at a medium scan speed of 20″/s, with a cross scan step of 5″ and scan leg length of 3′. The number of scan legs was 8 for the 63° scan, and 4 for the 45° scan. For spectroscopic observations, a 1669 s PacsLineSpec (obsid 1342186314) and a 5150 s PacsRangeSpec (obsid 1342187019) were obtained. The PacsLineSpec provides two simultaneous spectra at wavelengths 62.93−64.43 μm, and 180.76−190.29 μm. PacsRangeSpec covers six spectral ranges of 71.81−73.28 μm, 78.37−79.76 μm, 89.28−90.48 μm, 143.59−146.53 μm, 156.70−159.47μm, and 178.51−180.96 μm. Observations were taken in the chop-nod mode, with a narrow 2″ dither. The target was centred at the central spaxel of the 9.4″ × 9.4″ grid of the PACS integral field unit. The data was reduced using the Herschel interactive processing environment (HIPE; Ott 2010) developer build version 3.0.1212, and the data reduction scripts provided at the Herschel data reduction workshop held in January 2010.

For the photometric data, a mosaic was created from the two scan maps. Aperture photometry was performed using an aperture radius of 16″ in the blue, and 19.2″ in the red. An aperture correction of 0.922 (blue), and 0.885 (red) was applied to the photometry. The aperture corrections were obtained from the PACS PhotChopNod release note (Feb. 22, 2010). The flux calibration uncertainty is estimated to be 5% in the blue and 10% in the red. For the spectroscopic data, we extracted the spectrum from the central spaxel, and then applied an aperture correction in order to minimise the flux loss in the neighbouring spaxels. Spectra from the central spaxel were extracted for both the A and the B nods. We then applied wavelength binning to each nod spectrum, using a bin size that is half the width of the instrumental resolution. The final spectrum is the mean of the wavelength-binned spectra from the two nods. The absolute flux calibration uncertainty is estimated to be 40%. We have detected the [OI] 63.2 μm emission line for ET Cha, while all other lines are undetected (see Table 1). We used the IDL routine MPFITPEAK to fit an error-weighted Gaussian to the observed [OI] 63.2 μm line, and measured the integrated flux of the Gaussian line fit. The 1-σ error to the line flux was calculated by setting the height of the Gaussian equal to the continuum rms value, and the width equal to the instrumental resolution. The continuum emission at the rest wavelength of 63.18 μm was estimated by fitting a first-order polynomial to the spectral region.

3.2. CTIO/ANDICAM photometry

We obtained new simultaneous optical/NIR photometry using the dual-channel ANDICAM instrument on the CTIO 1.3 m telescope on March 6, 2009. Photometric calibration was performed using the PG 1047 Landolt field. Both the optical and near-infrared data were reduced using standard procedures (flat-fielding, cosmetic cleaning, shift-and add). The photometry on ET Cha and the photometric standard was extracted within a 3″ aperture and airmass and colour corrections were applying using coefficients from the ANDICAM website³ in the optical and from Frogel (1998) in the near-infrared. The resulting photometry is listed in Table 2 along with the mid- and far-infrared fluxes adopted in our analysis. The new photometric fluxes are substantially lower than previously published fluxes, see Table 2.

3.3. APEX/LABOCA photometry

An upper limit at 870 μm has been obtained from new continuum maps of the η Chamaeleontis association taken with the LABOCA bolometer array on APEX (Siringo et al. 2009). The data was reduced using the Bolometer array data Analysis (BoA) software, with a pipeline optimised for faint compact sources. Fluxes were extracted within BoA by fitting the amplitude of a beam-sized Gaussian at specified positions within a map with a pixel scale of 4.6″/px. For each target the flux was extracted at points in a 5 × 5 rectangular grid centred on the target with a spacing of 37″ (twice the beam FWHM). The sample standard deviation of the 24 off-source measurements is the 1-σ uncertainty quoted here. Using aperture photometry instead of beam fitting, with a variety of aperture sizes, and aperture corrections computed from the beam, yields very similar results (9 ± 11 mJy/beam for apertures with r = HWHM).

3.4. UV observations

ET Cha was observed in the far-ultraviolet with the Cosmic Origins Spectrograph (COS) and in the near-ultraviolet and blue with the Space Telescope Imaging Spectrograph (STIS) on board of the Hubble Space Telescope (HST) as part of the programme “Disks, Accretion and Outflows of T Tau Stars” (P.I. G. Herczeg) on 5 Feb. 2010. The far-ultraviolet observations were obtained with the G130M and G160M gratings with 3891s and 4501s integration times, respectively, to cover the 1140−1790 Å spectral range with a resolution of  ~18 000. The data were reduced with the COS calibration pipeline CALCOS and individual segments were combined with the IDL coaddition procedure described by Danforth et al. (2010). The near-ultraviolet and blue observations were obtained with the G230L and G430L gratings with integrations of 45 s and 680 s, respectively, to cover the 1700−5700 Å spectral range with a resolution of  ~2000. The data were reduced with custom-built software in IDL.

A FUSE spectrum of ET Cha covers the spectral region shortward of 1140 Å but has poor signal-to-noise. Instead, the 912–1140 Å flux was estimated from a high-quality FUSE spectrum of TW Hya (Johns-Krull & Herczeg 2007), scaled to the flux in the 1250–1700 Å bandpass. All different data have been smoothed to a common low resolution, see Fig. 2. Integrated fluxes over three different UV bands as listed in Table 3.

3.5. Re-analysis of high-resolution optical spectroscopy

A high-resolution optical spectrum of ET Cha was acquired on June 26, 2002 with the 3.9 m Anglo-Australian Telescope AAT and University College London coudé Echelle Spectrograph (UCLES) as first published by Lawson et al. (2004). A 1.5″ slit was used delivering a resolving power of  ~30 000 (10 km s^-1 at 6300 Å) and covering the wavelength range between 4980 and 9220 Å. The spectrum was calibrated using dome-flats, bias frames and ThAr arc frames, making use of standard library routines such as ccdproc within IRAF (see Lawson et al. 2004, for a detailed description). We have paid particular attention to the removal of telluric contributions to these lines. We removed the OI 6300 Å telluric contribution using the RECX 10 spectrum, a star of similar spectral type but without any evidences of an accretion disk or outflow (Lyo et al. 2003), meaning that the measured OI 6300 Å emission from RECX 10 is thoroughly due to the atmosphere of the Earth.

The average air-masses during observations of ET Cha and RECX 10 are 2.6 and 2.7, respectively. We made a [OI] 6300 Å-map by first removing the [OI] feature from the RECX 10 spectrum, and then subtracting this edited spectrum from the original RECX 10 spectrum. The ET Cha spectrum was divided by the [OI] 6300 Å-map to remove the OI telluric line. We used the photospheric Li I absorption lines at 6707.76 Å and 6707.91 Å to measure a stellar radial velocity of  ~22 km s^-1 (~34.6 km s^-1 after heliocentric correction).

Figure 3 shows the observed profiles of ET Cha in the oxygen line [OI] 6300 Å, with 3 other optical forbidden emission lines overplotted that trace outflows, plus a Li I absorption line to determine the systematic stellar velocity. The [OI] 6300 Å line shows a broad component centred around the stellar systematic velocity (low-velocity component LVC), and a blue component shifted by about 42 km s^-1 (high-velocity component HVC). We have fitted the HVC and LVC by two Gaussian profiles. Measured equivalent widths are listed in Table 4. Equivalent widths for these lines are converted to line luminosities using the procedure outlined in Hartigan et al. (1995), assuming a distance of 97 pc and zero visual extinction. As expected for outflows, the [NII] 6583 Å, [SII] 6716 Å and [SII] 6731 Å emission lines emphasise the HVC, and are all blue-shifted by about the same margin.

3.6. Projected stellar rotational velocity

The high resolution spectrum presented in Sect. 3.5 was additionally used to determine the projected rotational velocity of the star v_rotsini. This quantity is typically derived from empirical relations between the full-width at half-maximum of stellar absorption lines and v_rotsini (see e.g. Martínez-Arnáiz et al. 2010). However, in this paper we use the Fourier transform of the line profile (Gray 1992). The full power of this method is revealed in case of very high spectral resolution and signal to noise (Reiners & Schmitt 2003), but has also been applied successfully to moderate and low quality data (Mora et al. 2001), using the following simplifications.

The location of the first minimum of the Fourier transform of any isolated photospheric line profile provides a measurement of vsini. If many lines are available, the average value and standard deviation can be estimated. The advantage of this method is that the measurements are independent of empirical calibrations and synthetic spectra, provided that enough lines can be identified and the dominant line broadening mechanism is stellar rotation. In terms of Fourier analysis, the latter requirement is equivalent to the following condition: the first minimum of the Fourier transform of photospheric line profiles must be at a lower frequency than the Fourier transform of lines not affected by rotation velocity (e.g. arc or telluric lines). In addition, the algorithm is robust against the introduction of noise, and capable of working with low signal to noise data.

The method has been applied to the visible UCLES spectrum of ET Cha as described in Sect. 3.5. Few lines could be used for the analysis, due to the low signal to noise ratio (~10 for unbinned data). A synthetic Kurucz spectrum with physical parameters similar to those listed in Sect. 4.1 was used to identify isolated atomic photospheric lines. After clipping bad lines, only 16 lines remained in the final average used below. Figure 4 shows a sample Fourier transform of both a photospheric and a telluric line to illustrate the method. We obtain ${\mathit{v}}_{\mathrm{rot}}\sin \mathit{i}\mathrm{=}\mathrm{(}\mathrm{15.9}\mathrm{\pm }\mathrm{2.0}\mathrm{)}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}$(1)

3.7. APEX ¹²CO J = 3 → 2

The source was studied in the J = 3 → 2 rotational transition of ¹²CO at 345.796 GHz using the APEX telescope and heterodyne receiver system twice, in April and July 2010, for a total on-source time of 100 minutes. A standard on-off switching scheme was used, with reference position 180″ off in Azimuth. With an aperture efficiency of η_mb = 0.6, the effective area of the telescope is 67.9 m², and the beamsize Ω_A at 345.796 GHz is 13.7″. The CO J = 3 → 2 line is commonly detected from disks of radii 100 AU or more; however, in the case of ET Cha, the line was not detected (Fig. 6), down to an rms of 17 mK (antenna temperature) for 0.425 km s^-1 wide velocity bins. Assuming the underlying linewidth of the emitting region is 7.5 km s^-1, as predicted by the best fitting model (Sect. 5) and typical of CO lines from disks, the upper limit to the line flux is 3σ = 5.5 × 10^-20 W/m².

4. Disk modelling

The strong near–mid IR continuum excess of ET Cha (see Fig. 7) as well as the narrow H₂ v = 1 → 0 S(1) emission line at 2.122 μm (Ramsay Howat & Greaves 2007) indicate the presence of a gas-rich protoplanetary disk. The optical high-resolution spectrum with blue-shifted emission lines shows, however, that in addition to the disk, there also exists an outflow that might contribute to the gas emission lines. In the following, we present a combined gas+dust disk model to explore in how far a disk model alone can explain all photometric and spectroscopic observations of ET Cha. The possible contribution of an outflow is discussed in Sect. 6.4.

The disk modelling procedure in this paper consists of three phases. In phase 1, we make a Bayesian analysis, based on the new BVRIJ photometric data, to obtain best fitting values for the stellar parameters. In phase 2, we fix a few more parameters like the inner disk radius based on physical assumptions, e.g., according to stellar luminosity and dust sublimation temperature, and in phase 3, we use the thermo-chemical code ProDiMo (Woitke et al. 2009; Kamp et al. 2010) to calculate the gas- and dust temperature structure in the disk, the chemical and ice composition, and to fit all remaining disk, dust, and gas parameters of the model to the observations as good as possible.

4.1. Stellar parameters

While our new JHK photometry is consistent with the 2 MASS photometry, the VRI photometry is about 0.5−0.7 mag fainter than the data published in (Lawson et al. 2002). Indeed, our simultaneous photometry appears to represent the “faint” state for ET Cha which, in the accretion hotspot scenario, is purely photospheric.

To estimate the stellar properties, we therefore adopt our new simultaneous photometry. We only fit the BVRIJ fluxes to ensure that our estimates are not biased by the near-IR excess from the disk. We computed a grid of photospheric models, varying T_eff (using the log g = 4.5 NextGen models from (Allard et al. 1997), A_V (using the R_V = 3.1-law from Cardelli et al. 1989), R_⋆ and the distance d to ET Cha within generous ranges. The optical extinction is defined by A_V = 2.5log ₁₀(e) τ_V where τ_V is the optical depth in the visual due to interstellar dust extinction. We note that the last two parameters are degenerate.

The results are shown in Fig. 5. A very good fit is obtained with T_eff = 3300−3400 K, d = 97−100 pc, R_⋆ = 0.84−0.89 R_⊙ (corresponding to L_⋆ = 0.085 L_⊙) and A_V ≤ 0.2 mag. This is in excellent agreement with the findings of Luhman & Steeghs (2004), including the effective temperature which they derived exclusively from their spectrum and not from photometry. For the disk modelling, we adopt T_eff = 3400 K and R_⋆ = 0.84 R_⊙, which is both the most probable combination of parameters and the closest to the spectroscopically derived stellar properties.

Finally, we note that conducting the same analysis on the previous photometric dataset (obtained in the “bright” state) yields both a hotter central star (T_eff ≈ 4000 K) and a substantial foreground extinction (A_V = 0.75 mag), both of which are inconsistent with our prior knowledge of the source (see Sect. 2).

In the following, we use the nearest NextGen stellar input spectrum for T_eff = 3400 K, L_⋆ = 0.085 L_⊙ and solar metallicities (resulting in log  g = 3.89 and R_⋆ = 0.84 R_⊙ for M_⋆ = 0.2 M_⊙). According to the stellar evolutionary models of Siess et al. (2000), our choice of T_eff and L_⋆ is consistent with a 6−8 Myr old star of mass (0.2−0.3) M_⊙ and solar metallicities.

4.2. Stellar UV-excess

In the UV, stellar activity and mass accretion create an excess emission with respect to classical stellar atmosphere models. Therefore we use our composite observed UV spectrum of ET Cha (Fig. 2) at λ < 450 nm as input instead. Although strong Lyα emission from ET Cha was detected with HST/COS, no emission was detected near line centre because of neutral hydrogen in our line of sight to the star. Based on an analysis of a similar Lyα emission profile from TW Hya (Herczeg et al. 2004), we roughly estimate that the total flux in this line is twice the detected value. Therefore, we have doubled the detected flux between 1210 Å and 1220 Å for the construction of our UV input spectrum “as seen by the disk”. This yields an integrated flux of 2.22 × 10^-12 erg/cm²/s in the 1110−1450 Å band, in contrast to Table 3. We have converted these integrated UV fluxes into additional photometry points, see Table 2, and plotted these points with large errorbars in Fig. 7. These data result in a fractional UV luminosity f_UV = L_UV/L_⋆ = 0.017 with the UV luminosity L_UV being integrated from 912 Å to 2500 Å as introduced by Woitke et al. (2010). The “extra UV” is not important for the spectral energy distribution (SED) modelling, but is essential for the gas modelling, because the UV irradiation is a decisive factor for the gas heating and chemistry in the disk.

4.3. Disk inclination and inner radius

For the disk modelling in Sect. 4.5 we have fixed two parameters to reduce the dimension of parameter space for the fitting problem, namely the disk inclination and inner radius.

First, the disk inclination is assumed to be i = 40° as measured from face-on orientation. As long as the disk is not intersecting the line of sight to the star, this parameter does not have a strong impact on the model results concerning both SED-analysis and line flux predictions, and we do not have clear observations, e.g., images, that would allow us to determine this quantity unambiguously. High values for the disk inclination are supported by the H_α line analysis (60°, Lawson et al. 2004), the low outflow velocities observed (Sect. 3.5), and by the strongly rotation broadened stellar absorption lines (Sect. 3.6), favouring an edge-on rather than a face-on disk orientation.

However, our disk modelling suggests big disk scale heights in the inner disk regions, needed to reproduce the strong near-IR excess (see Fig. C.4). For inclinations of 60° and higher, the observer’s line of sight to the central star would therefore intercept the disk. This would cause a dramatic reduction of the observed fluxes at optical and UV wavelengths, as well as a strong reddening of the optical colours. This is inconsistent with our initial assumption that the star is seen directly, used to estimate the stellar properties. In principle, one could increase both the intrinsic stellar luminosity and effective temperature of the star to compensate exactly for this reddening, but to match the observed SED, the required T_eff is well beyond values consistent with the M3–M3.5 spectroscopic classification of ET Cha. For inclinations  ≳60° we furthermore observe that the 10 μm and 20 μm silicate features are no longer in emission, which clearly contradicts of observed SED of ET Cha.

Second, concerning the disk inner radius, the strong and continuous near-IR fluxes of ET Cha suggest a disk which extends inward close to the star, with no dust holes or gaps. T Tauri disks are generally assumed to be truncated by the stellar magnetosphere near the co-rotation radius, with material accreting along magnetic field lines onto high-latitude regions of the star (Königl 1991; Shu et al. 1994). From the analysis of optical absorption lines of ET Cha (Sect. 3.6) we have inferred a projected stellar rotation velocity of v_rotsini = (15.9 ± 2.0) km s^-1. At i = 40°, this translates into v_rot = (25 ± 3) km s^-1 (about 12% of the break-up velocity) and, by assuming M_⋆ = 0.2   M_⊙ and R_⋆ = 0.84 R_⊙ (Sect. 4.1), into a stellar rotation period of P = (1.7 ± 0.2) days, and a co-rotation radius of (0.015 ± 0.002) AU.

However, the dust temperature at the co-rotation radius (0.015 AU) is higher than 1500 K, associated with the sublimation temperature of silicates. Therefore, we have assumed a slightly larger inner disk radius, R_in = 0.022 AU, where the dust temperature is about 1500 K.

4.4. A single disk model

For a single disk model, we use the radiation thermo-chemical disk code ProDiMo (Woitke et al. 2009; Kamp et al. 2010) to calculate the dust continuum radiative transfer, and the gas thermal balance and chemistry throughout the disk. We use 10 elements, 76 gas phase and solid ice species, and 992 reactions including a detailed treatment of UV-photorates (see Kamp et al. 2010), H₂ formation of grain surfaces, vibrationally excited H$\begin{array}{}\mathit{\star }\\ \mathrm{2}\end{array}$ chemistry, and ice formation (adsorption, thermal desorption, photo-desorption, and cosmic-ray desorption) for a limited number of ice species (see Woitke et al. 2009, for details). We also use ProDiMo in this paper to compute all observables including the SED, images, and gas emission line fluxes and profiles. Latest improvements to the ProDiMo model include X-rays chemistry and heating, a parametric prescription for dust settling, UV fluorescent pumping, PAH ionisation and heating/cooling, [OI] 6300 Å pumping by OH-photodissociation, H₂-pumping by its formation on grain surfaces, formal solutions of the line transfer problem, and chemical heating. These improvements are explained in Appendices A.2 to A.8.

Our modelling of ET Cha is based on a prescribed disk density structure, using power-laws for the surface density distribution and disk flaring, see Appendix A.1. This approach allows for a rapid model computation (avoiding the structure iteration loop, see Fig. 1 in Woitke et al. 2009), and is hence more suitable for a deep search for fitting values in parameter space. In this mode, the disk code uses altogether 25 physical parameters, most of which are considered as fixed for the modelling of ET Cha (for instance the stellar properties, see Table 5).

4.5. Parameter fitting procedure

The determination of the remaining free disk, gas and dust parameters, like the total disk gas mass for instance, is a key objective of the modelling of individual targets. These parameters are “determined” in this paper by considering the values (better the ranges of values) that lead to a satisfying match between predictions and observations in the frame of the model, henceforth called “solutions”. Practically all determinations of properties of astrophysical objects work that way, even if it’s just a simple one-dimensional dependence between property and observable, because the desired quantities are rarely accessible via direct observations. This is the well-known problem of model inversion in astronomy (Lucy 1994), with all the usual shortcomings and concerns which can be subdivided into four families:

• 1.

  Concerns about the quality of the model itself (missing physics,poorly determined input quantities like cross-sections and ratecoefficients, numerical issues like grid resolution,etc.).

• 2.

  Concerns about the completeness of solutions found in parameter space, for instance, how can we be sure of having found the best solution?

• 3.

  Under-determination, i.e. a weak dependence of observables on model parameters in combination with considerable uncertainties in the observables.

• 4.

  Non-uniqueness of inversion procedure (multiple solutions).

Given the complexity of the disk model at hand, with 8−14 free parameters for the modelling of the disk of ET Cha, we have to take into account the possibility that different solutions exist which practically achieve the same degree of agreement between model results and observations. Whereas in a low-dimension parameter space this seems to a be a weird, seldom special case, it occurs frequently in high dimensions, with numerous local minima. The manifold of solutions in a multi-dimensional parameter space can be (and usually is) astonishingly rich and complicated in structure.

Furthermore, it is important to note that, in a high-dimensional parameter space, an exhaustive search is practically impossible. This is even more so for our disk modelling, since one complete model takes about 1 CPU hour on a single 3 GHz processor machine. For 10 parameters, with 20 well-selected values around a main solution each, one would need to run 20¹⁰ ≈ 10¹³ models which would take about 1.2 × 10⁹ CPU yrs.

Therefore, we are not able to fully devitalise the concerns (1) to (4), but have to face the fact the any parameter determination by model inversion in high-dimensional parameter space is an intrinsically uncertain business. Our strategy in this paper is as follows. We use an evolutionary strategy to find a handful of well-fitting disk models, corresponding to different local minima in parameter space. Among these, we select a “best-fitting” model, the properties of which are described in detail in Appendix C. The best-fitting model is then re-run with different choices of input physics in Appendix D, to explore the principal effects of e.g. X-rays and the treatment of H₂-formation.

We then conclude about the confidence intervals of derived parameter values in Sect. 5 by considering (i) small deviations of single parameters around the best fitting solution, (ii) overall experience from fitting by hand and variance of different solutions found by the evolutionary strategy, and (iii) direct constraints from physical arguments like optical depth and column densities. We admit that this modelling procedure is not entirely satisfactory. A more complete discussion of our modelling strategy will be the topic of a forthcoming paper.

4.6. Fit quality and evolutionary strategy

Mathematically, all model inversion techniques can be formulated in terms of a certain strategy to minimise χ², i.e. to find a minimum of the deviations between model predictions and observations in parameter space. In this paper, we use the following logarithmic measure of these deviations as $\begin{array}{ccc}& & {\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Delta }}_{\mathit{i}}^{\mathrm{2}} \mathrm{with}\\ & & {\mathrm{\Delta }}_{\mathit{i}}\mathrm{=}\{\begin{array}{c}\\ & & \\ & & \end{array}\end{array}$(2)${\mathit{F}}_{\mathrm{obs}}^{\mathit{i}}$ and ${\mathit{\sigma }}_{\mathit{F}}^{\mathit{i}}$ are an observed flux and its uncertainty, and ${\mathit{F}}_{\mathrm{mod}}^{\mathit{i}}$ is the predicted flux by the model. The logarithmic nature of our χ² is motivated by the need to assign an equally high (bad) number to χ if a model flux is a factor 10 too strong as if it is a factor 10 too weak, analog to considering deviations in magnitudes⁴.

Equation (2) is applied separately to all photometric data, all λ-points in the Spitzer spectrum, and to the line fluxes of CO J = 3 → 2, [OI] 63 μm, [OI] 6300 Å (LVC) and o-H₂ 2.122 μm, resulting in ${\mathit{\chi }}_{\mathrm{Phot}}^{\mathrm{2}}$, ${\mathit{\chi }}_{\mathrm{Spit}}^{\mathrm{2}}$, and ${\mathit{\chi }}_{\mathrm{Line}}^{\mathrm{2}}$, respectively. The total fit quality of a model is then calculated as ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}{\mathit{w}}^{\mathrm{Phot}}{\mathit{\chi }}_{\mathrm{Phot}}^{\mathrm{2}}\mathrm{+}{\mathit{w}}^{\mathrm{Spit}}{\mathit{\chi }}_{\mathrm{Spit}}^{\mathrm{2}}\mathrm{+}{\mathit{w}}^{\mathrm{Line}}{\mathit{\chi }}_{\mathrm{Line}}^{\mathrm{2}}$(3)with adjustable weights w^Phot, w^Spit and w^Line, normalised as w^Phot + w^Spit + w^Line = 1.

We have applied the (1,11) -evolutionary strategy with adaptive step-size control of Rechenberg (2000) to minimise χ, i.e. to find best-fitting values of our remaining free gas, dust and disk parameters. The strategy uses 1 parent producing 11 offsprings with slightly modified parameters, the best of which will become the parent of the next generation (the parent always dies). The step-size δ is transmitted to the children and treated as additional parameter to be optimised. After some experiments, we have chosen w^Phot = 0.35, w^Spit = 0.45 and w^Line = 0.2 for optimum performance of the evolutionary strategy.

Figure 8 visualises one exemplary run of the evolutionary strategy, showing the changing model parameters and results over 300 generations. We started from a generic disk setup with M_gas = 10^-3   M_⊙, gas/dust ratio of 100, and an outer disk radius of R_out = 100 AU. This model fails badly to explain the SED and, in particular, the line fluxes. CO J = 3 → 2 is about 10 times too strong and both [OI] 6300 Å and o-H₂ 2.122 μm are more than a factor of 100 too weak. However, after about 300 generations, the model has achieved a good fit of all available continuum and line observations. The final parameter values are far from these initial guesses, yielding a disk that is less than 10 AU in radius and much less massive. All models from about generation 50 onwards fit the dust observations about equally well, the SEDs from these models are almost indistinguishable by eye when plotted as in Fig. 7. A good fit to the line observations was achieved only from generation 200 onward, after the disk radius has shrunk to about R_out ≲ 10 AU while continuing to fit the dust observations about equally well. Thus, the line observations can help to break the degeneracy of SED-fitting.

5. Results

We identify the result of the evolutionary run depicted in Fig. 8 as our main “best-fitting” model. The resulting parameters are listed in Table 5, the computed line fluxes are summarised in Table 6 and the resulting SED is plotted in Fig. 7. More details about the internal physical and chemical structure of this disk model are shown in Appendix C.

However, the identification of a best-fitting model was not at all straightforward. Altogether 17 runs of the evolutionary strategy have been executed with about 50 to 300 generations each, choosing different parameters to be varied, different initial guesses of the model parameters, or using different setups of the evolutionary strategy. Not all of these runs succeeded to find well-fitting solutions. Some almost equally well-fitting solutions are listed in Table 7. Based on these diverse model inversion results, we have to be very careful with conclusions about the disk properties of ET Cha. The following section summarises our confidence intervals for the various disk shape, dust and gas parameters and discusses which observations are key for their determination.

5.1. Dust mass

ET Cha’s spectral energy distribution (SED) is characterised by a strong near and mid IR excess relative to the star, similar to other much more massive T Tauri disks, see Fig. 7. However, in the far-IR, the fluxes are very faint. 69 mJy at 160 μm was too weak to be detected by Spitzer, and we have not detected the object at 870 μm with APEX. Thus, the new Herschel photometry points at 70 μm and 160 μm allow for the most complete SED analysis of the source to date.

In all computed disk models for ET Cha, the disk is optically thin at 160 μm. The observed flux at distance d is hence given by ${\mathit{F}}_{\mathit{\nu }}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}{\mathit{B}}_{\mathit{\nu }}(\mathrm{⟨}{\mathit{T}}_{\mathrm{dust}}\mathrm{⟩})\hspace{0.17em}{\mathit{M}}_{\mathrm{dust}}}{{\mathit{d}}^{\mathrm{2}}}$(4)where  ⟨ T_dust ⟩  is the dust mass averaged dust temperature, ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}$ the dust absorption coefficient in  [cm²/g(dust)]  and B_ν the Planck function. In the best-fitting model we measure  ⟨ T_dust ⟩  ≈ 56 K, and ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{=}\mathrm{34}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\mathit{/}\mathrm{g}\mathrm{\left(}\mathrm{dust}\mathrm{\right)}$ at 160 μm, which, according to Eq. (4), results in a 160 μm flux of 48 mJy which is in good agreement with both, the computed flux from the full disk model (51 mJy) and the observed value of 68 mJy.

Therefore, the observed 160 μm-flux leaves no doubt that the mass of the dust in the disk of ET Cha must be low, M_dust ≈ 2.6 × 10^-8 M_⊙ according to our best fitting model. More massive disks produce too strong 160 μm continuum fluxes. Although there are certainly some ambiguities in the dust mass determination, e.g. if other dust size parameters are used, leading to a different ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}$, but among all parameters, M_dust certainly belongs to those which are only little influenced by other parameters and can be determined with some confidence. Across all SED-fitting disk models, with varying dust composition and size parameters, the value for M_dust ranges in (2−5) × 10^-8 M_⊙.

5.2. Dust disk characteristics and particle sizes

According to the low dust mass we derive, the disk of ET Cha has only low optical depths across the disk. Our best-fitting model has vertical dust absorption optical depths  <1 for λ ≳ 4 μm, except for the optically thick 10 μm and 20 μm silicate features. This is valid for all radii since the column density exponent ϵ ≈ 0. However, the radial optical depths are much larger. The midplane radial dust extinction optical depth at 1 μm is about 150. Therefore, the disk of ET Cha is on the borderline between optically thin and thick. It has a complex dust temperature structure due to radial shielding effects (see Fig. C.2), but the vertical dust emission could be treated in the optically thin limit at most wavelengths.

Some conclusions about the dust particle sizes can be drawn from the shape of the 10 μm and 20 μm silicate features, which are clearly seen in emission for ET Cha (Fig. 7). We need a strong amplitude of opacity variation across these features (see Fig. C.1) to model the SED of ET Cha, which favours small, (sub-) micron sized dust particles. Since the peaks are optically thick, warm and small grains must be located in front of cooler dust along the line of sight. We notice that it is easier to fit the observed silicate features with low values of the column density power-law index ϵ, i.e. a roughly flat surface density distribution, and small disk inclinations. For larger ϵ or higher disk inclinations (closer to edge-on), the silicate features weaken and eventually vanish.

If dust and gas are well-mixed (no dust settling, e.g. models 3 and 4 in Table 7), the models clearly favour very small particles. We can fit the entire SED with a uniform dust population that is either truncated at about 1 μm, or, alternatively, with a continuous size distribution ranging from 0.05 μm to 1 mm, but with an unusually large power-law index of p = 4.1. The second approach allows for slightly better SED-fits. These findings with ProDiMo have been carefully checked against the Monte Carlo radiative transfer code MCFOST (Pinte et al. 2006, 2009), showing a very good agreement in calculated dust temperatures and continuum fluxes.

However, if dust settling is taken into account (in the approximate way explained in Appendix A.2), Table 7 demonstrates that lower p = 3.9 or even p = 3.2 are also possible, in which case the volume-integrated dust size distribution in ET Cha would not be unusual at all, close to the default value of p = 3.5. Our honest conclusion about ET Cha is hence that it’s dust must predominantly be made of (sub-) micron particles at the disk surface, where the silicate emission features form. We have therefore decided to refrain from quoting any errorbars to our results concerning the dust size parameters.

5.3. Dust composition

In our best fitting model, we assume an effective mix of about 33% amorphous fosterite Mg₂SiO₄ (Jäger et al. 2003), 24% amorphous carbon (Zubko et al. 2004), 23% amorphous olivine MgFeSiO₄ (Dorschner et al. 1995), 9% amorphous silica SiO₂ (Posch et al. 2003), 8% amorphous enstatite MgSiO₃ (Dorschner et al. 1995), and 3% crystalline fosterite Mg₂SiO₄ (Servoin & Piriou 1973). The citations indicate our sources for the complex refractory indexes of the various pure materials. The dust absorption and scattering opacities are calculated by applying effective mixing theory (Bruggeman 1935) and Mie theory, based on these optical data and volume fractions, and assuming that the chemical dust composition and unsettled dust size distribution is unique throughout the disk. The inclusion of crystalline fosterite is motivated by the fine-structure of the observed second silicate feature at 20 μm (Bouwman et al. 2006), see (Sicilia-Aguilar et al. 2009), which shows several narrow peaks close to the fosterite peak positions in the data of Servoin & Piriou (1973). The resulting volume fractions are a by-product of our automated fitting procedure from additional runs of the evolutionary strategy not shown in Table 7. We do not claim, however, to have determined the dust composition of ET Cha, as our method is focused on fitting the overall shape of SED rather than individual dust features.

The inclusion of 24% amorphous carbon, however, was an important step to understand the SED of ET Cha. A comparison model without amorphous carbon (green line in Fig. 10) demonstrates its impact on the SED. First, amorphous carbon reduces the dust albedo at UV to near-IR wavelengths (see Fig. C.1). Pure laboratory silicates have an albedo of about 80−99% around 1 μm. This leads to a substantial starlight amplification via scattering by the disk at UV to near-IR wavelengths, which is inconsistent with the photometric data. Second, amorphous carbon enhances the absorption and thermal emission in the near-mid IR. As a consequence, the disk is more effectively heated by star light, and produces more thermal emission shortward of the 10 μm silicate feature, just where ET Cha is very bright, resulting in a much better fit of the SED if we include amorphous carbon. However, any other kind of impurities or inclusions, metallic iron for instance, would cause a similar increase of the dust absorption opacities at optical to near-IR wavelengths, amorphous carbon is just one of the options. The formation of small alien inclusions in the dust material has been demonstrated by Helling & Woitke (2006) and Helling et al. (2006) to be a natural consequence of the refractory dust formation process in somewhat different oxygen-rich environments, namely in brown dwarf atmospheres. Concerning ET Cha we conclude that strong near-IR dust absorption opacities are needed to fit the SED, substantially stronger than those of pure amorphous or crystalline silicates which have a “glassy” character.

5.4. Outer disk radius

Any outer disk radius, at least up to 200 AU, works fine to fit the SED alone. However, the non-detection of CO J = 3 → 2 and the modest [OI] 63 μm line flux put severe constraints on R_out. As Figs. 8 and 9 demonstrate, the CO line flux depends very strongly and robustly on R_out. All calculated models with R_out > 25 AU would violate the 3σ upper limit of CO J = 3 → 2. Only models with R_out < 10 AU are consistent with the 1σ CO line flux upper limit. The CO lines are extremely optically thick, τ_line ≈ 10⁵ (see Sect. 5.5). Therefore, even if the CO abundance was reduced by a factor of 10³ throughout the disk, for example by very efficient CO ice formation, the conclusion that the disk of ET Cha must be small would still be valid. The relatively robust determination of R_out ≈ (6−9) AU is demonstrated by Fig. 9 and Table 7. Lower values of R_out are actually inconsistent with the measured [OI] 63 μm line flux (unless interpreted as originating form an outflow, see Sect. 6.4), because the outer disk is responsible for the [OI] 63 μm line emission (see Sect. 5.5).

5.5. Spatial origin and characteristics of gas emission lines

Before we continue to conclude about the determination of the disk gas mass of ET Cha, we first have to clarify where the observed spectral lines come from and what they tell us. Figure 11 shows the calculated line fluxes and profiles of the three detected lines [OI] 63 μm, [OI] 6300 Å (LVC) and o-H₂ 2.122 μm. We also plot here the vertical line and continuum optical depths and the cumulative line flux as function of radius, to facilitate the discussion of the lines’ spatial origin. These plots are based on the best-fitting disk model, but the drawn conclusions about the spatial origin of the spectral lines are quite general and valid for all calculated disk models that fit the observations of ET Cha.

The far-IR [OI] 63 μm line probes the outer disk layers, about 2−7 AU in this tiny model, and is optically thick with vertical optical depths τ_line ≈ 100−1000 in this radial region⁵. The continuum is optically thin. Since the line is collisionally excited, with an excitation energy of about 228 K, the line flux probes first and foremost the existence of warm gas (T_gas ≳ 50 K) in the disk surface, here at relative heights z/r = 0.4−0.5. In this region, PAH heating is usually the most important heating process. Therefore, this line shows a strong correlation with the assumed PAH abundance, f_PAH, UV irradiation, and disk flaring.

The optical [OI] 6300 Å line originates from hot gas in the inner disk, about 0.06−0.7 AU in this model. As a forbidden line, with very low Einstein coefficient A₄₁ ≈ 6.5 × 10^-3 s^-1, it is mostly optically thin τ_line ≈ 0.1−1. With an excitation energy of  ~15 900 K, it probes the existence of hot and dense gas T_gas ≳ 3000 K in front of an optically thick continuum. To our surprise, this line is not significantly excited by OH photo-dissociation in the model. Neglecting this pumping effect (see Appendix A.5) weakens the [OI] 6300 Å line flux by only 14%. Thus, the [OI] 6300 Å line comes from the bottom of the hot atomic layer, from about z/r = 0.3−0.4 in this model. As soon as molecules form, for example OH and CO, the temperature drops significantly by molecular line cooling, and the [OI] 6300 Å line cannot be excited any longer. Since the line is optically thin, its flux reacts quite sensitively on disk mass. We note that the predicted line FWHM = 35.7 km s^-1 is in excellent agreement with the observations, FWHM = (38 ± 15) km s^-1, meaning that the spatial origin of the [OI] 6300 Å line is about correct in the model.

The near-IR o-H₂ v = 1 → 0 S(1) 2.122 μm line also probes the inner disk regions, 0.05−0.6 AU in this model. The line is quite optically thick, τ_line ≈ 5−10 in this region, where the dust is also optically thick. Since its Einstein coefficient is extremely small, A₆₁ ≈ 7 × 10^-7 s^-1, the line needs high H₂ column densities of order 10²³ cm^-2 to become visible over the continuum level. Such column densities are only reached in quite deep layers, z/r ≈ 0.15, so this line forms in much deeper layers than the [OI] 6300 Å line. The temperature contrast between gas and dust is already quite weak in these layers (Fig. C.2), which limits the o-H₂ 2.122 μm line flux. The excitation energy of the line is about 7000 K, i.e. T_gas ≳ 1000 K is required to collisionally excite it. In the best-fitting model such conditions, T_gas ≳ 1000 K ≫ T_dust, are provided by the effect of exothermic chemical reactions (see Appendix A.8) which is active in the warm and dense gas close to the inner rim, deep but not too deep, just where this line forms. Therefore, we see a clear correlation between the o-H₂ 2.122 μm line flux and the assumed heating efficiency of exothermic reactions γ^chem. The line is also substantially pumped by H₂-formation on grains (see Appendix A.6). When neglecting this effect, the line attains a flux of only 35% of the value from the full model. We conclude that the H₂ line is sensitive to a temperature contrast between gas and dust, and H₂-formation, in quite deep disk layers close to the inner rim. However, the predicted line is too broad (FWHM = 37 km s^-1) as compared to the observations (FWHM = 18 ± 1.2 km s^-1). We were unable to find any disk model, among the  ~20 000 models computed, that shows such a narrow o-H₂ 2.122 μm line. Thus, the line forms too close to the star in the model, and we must be very careful when drawing our conclusions about the nature of ET Cha from the observed o-H₂ 2.122 μm line.

The sub-mm CO J = 3 → 2 line (not depicted in Fig. 11) probes the size of the disk and the gas temperature in the outermost disk layers. It is massively optically thick, τ_line ≳ 10⁵ in the line forming region 2−7 AU, where the dust is optically thin. Its flux is roughly proportional to the projected disk area times the gas temperature at relative heights z/r ≈ 0.4−0.5 in these outermost parts of the disk.

5.6. Disk gas mass and disk shape

The determination of the total disk gas mass of ET Cha in this paper is based on three detected gas emission lines, namely [OI] 63 μm, [OI] 6300 Å (LVC), and o-H₂ v = 1 → 0 S(1) 2.122 μm, which probe complementary radial and vertical disk regions. However, the first line, [OI] 63 μm, is massively optically thick, and the latter two lines mainly probe the hot gas in the inner disk regions. CO J = 3 → 2 was not detected. These circumstances immediately suggest that a robust gas mass determination of ET Cha is difficult.

Table 7 shows that various solutions in form of well-fitting disk models can be found where the total gas mass differs by up to two orders of magnitudes (factor 300). A careless reading of Table 7 would suggest that the particular values for M_gas, ϵ and β do not matter much, but this is not true. Careful inspection of the solutions in Table 7 by systematic variation of selected parameters (see e.g. Fig. 9 for the best-fitting model) shows that these solutions all represent well-defined local minima, where slight changes of any parameter leads to a considerable deterioration of the combined line + continuum fit quality χ. If one would measure the uncertainty of gas mass determination from these dependencies alone (for instance the deviation ΔM_gas where χ doubles) one would arrive at relatively small errors, about 40%. But such an error estimate would be misleading as it would not take into account the complicated manifold of local minima in parameter space.

All solutions in Table 7 fit our entire set of line and continuum observations about equally well. The different values for the gas mass, however, come in certain combinations with other disk shape parameters, like the column density power-law exponent ϵ and the flaring power β. Certain fine-tuned combinations of M_gas, ϵ and β do apparently all provide the proper density and temperature conditions in the disk that result in almost exactly the same observables, for example low M_gas in combination with high values for ϵ and β, or vice versa.

These relations can be understood by minimum column densities of warm H₂ and atomic oxygen in the inner disk parts that are inevitably required to make the [OI] 6300 Å (LVC) and o-H₂ 2.122 μm lines visible over the strong continuum at optical and near IR wavelengths. For example, in the least massive model 1 in Table 7, there is so little gas in the disk that a relatively large value of ϵ = 1.16 is required to concentrate the mass in the inner disk parts and so to provide the necessary column densities there. All line fluxes are in good agreement with the observations then, but the model has problems to find a good fit of the 10 μm and 20 μm silicate features. Also the flaring index of β = 1.33 is quite extreme, leading to a relative disk height of about z/r = 1.2 in the outer disk parts – such a “disk” would be taller than wide. From these arguments, we will discard model 1 from our selection of valid solutions in the following, and claim a minimum gas mass of about 5 × 10^-5   M_⊙ to achieve the necessary minimum gas column densities in the inner disk regions and to fit the silicate dust emission features simultaneously.

Concerning the other direction, we find it hard to provide any solid argument why the disk gas mass must not exceed a certain maximum value. In combination with our quite robust dust mass determination of M_dust ≈ 3 × 10^-8 M_⊙, however, even our minimum gas mass of 5 × 10^-5   M_⊙ already implies a very high gas/dust ratio of 2000. An even higher gas mass would translate into an accordingly higher gas/dust ratio. The highest value we found with our evolutionary strategy is 2.5 × 10^-3   M_⊙.

Summarising our modelling efforts, after having tried about 20 000 disk models and various choices of input physics, we conclude that the disk gas mass of ET Cha is only little constrained, our confidence interval is about M_gas ≈ (5 × 10^-5−3 × 10^-3) M_⊙. The surface density exponent ϵ and the flaring parameter β are likewise poorly constrained, although the models favour small ϵ ≈ 0 to fit the silicate dust emission features.

In view of the outflow discussion (Sect. 6.4), we note that the unclear physical origin of [OI] 63 μm obviously increases our uncertainty in gas mass determination. We have set up a similar evolutionary optimisation run as depicted in Fig. 8, but now treating the [OI] 63 μm line flux, as emitted by the disk, as an upper limit with  <2 × 10^-18 W/m². This run resulted in a disk of similar mass as compared to the best-fitting model, M_gas = 5.5 × 10^-4 M_⊙, but of even smaller size R_out = 2.1 AU, with an [OI] 63 μm line flux of  ~4 × 10^-18 W/m². We conclude that the [OI] 63 μm emission line flux is not crucial for the disk gas mass determination of ET Cha, because this disk is tiny and the hot gas responsible for our detected [OI] 6300 Å (LVC) and o-H₂ 2.122 μm lines is located in the inner regions. For other, more extended protoplanetary disks, the [OI] 63 μm line will be less optically thick and hence more useful for the purpose of gas mass determination (Kamp et al. 2010; Pinte et al. 2010; Woitke et al. 2010).

6. Discussion

6.1. Scale heights and unidentified heating

From the SED fitting, we get surprisingly robust results concerning the assumed vertical disk scale height, H₀ = (0.008−0.011) AU at reference radius r₀ = 0.1 AU, throughout all SED-fitting models. In the model, the reference scale height in combination with the flaring power β determines how much star light is captured by the disk. This light (unless scattered away) is absorbed by the dust in the disk surface and then thermally re-emitted, heating also the inner disk parts. Thus, H₀ regulates the dust temperature in the disk. Since ET Cha is very bright in the (3−8) μm region, we require large relative scale heights of the order of 10% close to the star to create enough dust emission from warm and hot grains to fit the near and mid IR.

Appendix D.4 (see lower left plot in Fig. C.4) demonstrates, however, that our prescribed scale heights are significantly larger, by a factor of 2−3, than those derived from self-consistent models where vertical hydrostatic equilibrium is assumed. This mismatch can be interpreted in three ways: (1) the close midplane regions of T Tauri disks are not in hydrostatic equilibrium, (2) the midplane temperatures are actually 4-9 times higher than assumed in the model (assuming $\mathit{H}\mathrm{\propto }{\mathit{c}}_{\mathrm{T}}\mathrm{\propto }\sqrt{\mathit{T}}$), or (3) there is an additional dust heating process active in the close disk midplane regions that leads to an additional energy flux through the dust component resulting in more observable near-mid IR photons without changing the temperatures much. Possibility 1 cannot be excluded per se. Possibility 2 seems unrealistic because it would cause the dust to evaporate. We favour possibility 3. Appendix D.5 shows that the inclusion of non-radiative dust heating via inelastic gas-dust collisions (thermal accommodation, driven by gas-dust temperature differences created through exothermic chemical reactions) leads to similar effects than increasing the scale height. Because of the ρ²-scaling of chemical reactions, this additional heating affects the dust primarily in high-density regions, increasing the production of near IR photons needed to fit the SED of ET Cha with moderate scale heights. By means of an extra run of the evolutionary strategy, we found out that we can reduce the scale height by about 35% to fit the SED, if we include this effect.

More unidentified dust heating processes may be active in the close midplanes of T Tauri disks. However, viscous heating (according to the formulation by Frank et al. 1992) is not doing a particularly good job in explaining the scale height inconsistencies of ET Cha, see Appendix D.3. Viscous heating dumps additional energy  ∝ ρ, i.e., after volume-integration, preferentially into the outer cold regions, leading to the production of more far IR photons. Moreover, viscous heating simply has little effect on the resulting disk temperatures and observable continuum flux in case of low-mass disks as for ET Cha. In fact, the viscous heating according to the formulation by (Frank et al. 1992) produces artifacts in the uppermost tenuous disk layers, where the viscous heating  ∝ ρ cannot be balanced by any cooling  ∝ ρ².

6.2. Effects of other physical processes

Appendix D discusses a number of further physical processes and their influence on the model that have not been included or discussed so far. To summarise, we find that

• (i)

  X-rays have very little effect (see Appendix D.2) and the neglecting of X-rays in our main model for ET Cha is fully justified. An X-ray luminosity of L_X = 6 × 10²⁸ erg/s, as observed for ET Cha (López-Santiago et al. 2010), turns out to be much less important for the disk heating and ionisation as compared to the strong FUV of ET Cha, L_UV = 6.5 × 10³¹ erg/s as measured between 912 and 2500 Å, based on our HST/COS and HST/STIS observations;

• (ii)

  the treatment of H₂-formation on grain surfaces is one of the most important yet quite uncertain processes for astrochemical modelling. Appendix D.6 shows that different approaches to calculate the H₂-formation rate lead to a systematic uncertainty in the model for the computation of the o-H₂ 2.122 μm about a factor of 5. Other spectral lines are less effected.

6.3. Disk inclination and outflow velocity

The low blue-shifts (− 42 km s^-1) seen in the optical emission lines (Fig. 3) are quite unusual for the jets/outflows of T Tauri stars. Typical outflow velocities are  ≳100 km s^-1 (Hartigan et al. 1995). A possible explanation could be the projection effect in case of a strongly inclined disk (i.e., close to edge-on). Adopting 100 km s^-1 as a lower limit for the outflow velocity of ET Cha would suggest a disk inclination of  >65°.

However, as argued in Sect. 4.3, we find that inclinations in the range 0°−50° are in agreement with the observations, but higher inclination angles would result in a partial obscuration of the star by the disk, with dramatic effects on the SED. The SED analysis therefore leads to the conclusion that the inclination of the system, as measured from face-on, is 50° or less which, in turn, translates to an outflow velocity of 65 km s^-1 at most. If this interpretation is correct, this would make the outflow from ET Cha one of the slowest known outflows from a T Tauri star.

6.4. Outflow and disk lifetime

Our analysis of the optical [OI] 6300 Å and [SII] 6731 Å emission line profiles and fluxes (see Appendix B) results in an estimate of the outflow mass-loss rate of ET Cha of Ṁ_outflow ≈ 10^-9   M_⊙/yr. Such an outflow can contribute to the [OI] 63 μm emission, which would render the [OI] 63 μm line flux, as emitted from the disk, weaker as assumed in Sect. 4.5. In summary, Appendix B shows that both approaches, the simple energetic outflow analysis by Hollenbach (1985), as well as the shock models computed by (Hartigan et al. 2004), suggest that the outflow from ET Cha does contribute a substantial, if not dominant, fraction of the observed [OI] 63.2 μm emission line flux. However, without spatially or velocity resolved observations of the outflow, and its subsequent detailed modelling (rather than using some “template” shock models), it is presently impossible to assess the exact contribution of the outflow to the line emission.

Our estimate of the outflow rate of ET Cha is high compared to the total disk mass we derive, M_disk ≈ 6 × 10^-4   M_⊙, as it suggests a disk lifetime of only M_gas/Ṁ_outflow ≈ 0.7 Myr, inconsistent with the age of the Chamaeleontis cluster of  ~8 Myr. If we would assume a generic gas/dust ratio of 100 and take our dust mass determination for granted, M_dust ≈ 3 × 10^-8   M_⊙, it is even worse. The disk lifetime would be even shorter in this case, only 3000 yrs, way too short to be feasible, unless we are just observing a temporary but short-lived 100 ×  peak in outflow rate.

We also notice that the mass accretion rate of ET Cha was estimated by Lawson et al. (2004) to be equally high, Ṁ_acc ≈ 10^-9   M_⊙/yr as the outflow mass loss rate, leading to similar lifetime inconsistencies. Furthermore, a branching ratio of Ṁ_outflow  /  Ṁ_acc ~ 1 is highly unusual, a few percent seems to be a well-established value for T Tauri stars (see e.g. Hartigan et al. 1995).

Murphy et al. (2011) reported on highly variable H_α equivalent widths and, accordingly, mass accretion rate, for the old T Tauri stars in the η Cha cluster. A factor of 100 variation in the accretion rate is observed in one newly-identified halo member of η Cha. But only a simultaneous reduction of both Ṁ_acc and Ṁ_outflow would help to resolve the aforementioned lifetime inconsistency, which doesn’t seem very likely.

One possible explanation would be a massive but short-lived outflow due to a flare from a former epoch, when the mass accretion rate was at least 10 times higher. An outflow of 100 km s^-1 would need  ~50 yrs to reach a distance of 1000 AU. This would be in agreement with the somewhat slow outflow velocity we derive, 42−60 km s^-1, because the outflow might have slowed down ever since. However, the optical emission lines do not show any evidence for a red-shifted HVC, as one would expect for a symmetric bi-polar outflow. The only plausible explanation of the missing red-shifted HVCs is that these components from the far side are attenuated by the dust in the disk. According to our disk model, the disk is optically thick at 6300 Å up to the outer radius  ~10 AU (see Fig. 11). Therefore, the line emitting region of the blue-shifted optical emission lines must be quite small, less than 5 AU if seen under 60° disk inclination, which corresponds to an age of only 0.25 yr.

7. Summary and conclusions

This paper has reported on new observations of ET Cha with several instruments: Herschel/PACS, CTIO/ANDICAM, HST/COS/STIS, and APEX. In combination with published data from Spitzer, Gemini/Phoenix and AAT/UCLES, we have collected an unprecedented observational data set about this object, including photometry, UV spectra, high-resolution optical spectrum, near and mid IR spectra, and far IR and sub-mm line fluxes.

We have calculated united gas and dust models for the disk of ET Cha that can simultaneously fit all line and continuum observations, except for a too broad o-H₂ 2.122 μm emission line profile. The observations also show some blue-shifted components of optical emission lines that point to an outflow and are not included in the models.

This paper has explored the parameter space of the disk models by using an evolutionary strategy to minimise the discrepancies between model predictions and observations. The paper has also introduced a number of basic improvements to the ProDiMo disk modelling code concerning the treatment of PAH ionisation balance and heating, heating by exothermic chemical reactions, several non-thermal pumping mechanisms for selected gas emission lines, and formal solutions of the line transfer problem at given inclination (Appendix A).

From the disk modelling we find a rich variety of fitting disk models that can explain our observations about equally well. Some of the model parameters (like the dust mass and the outer radius) can be determined with some confidence whereas other parameters (like the disk gas mass) are poorly constrained:

• The new Herschel/PACS photometric fluxes at 70 μm and 160 μm constrain the disk dust mass of ET Cha to be about M_dust = (2−5) × 10^-8 M_⊙, putting the object at the borderline between optically thin and optically thick.

• Then strong near IR excess of ET Cha can be fitted with a disk that is truncated at R_in ≈ 0.022 AU (where T_dust ≈ 1500 K) which is located slightly outside of the co-rotation radius of 0.015 AU. The latter is calculated according to the assumption that the star rotates with a period P = 1.7 days as suggested by our v_rotsin(i) analysis of rotationally broadened stellar absorption lines.

• From the APEX CO J = 3 → 2 non-detection, we can infer, with confidence, that the disk of ET Cha must be tiny in radius. The models favour an outer disk radius as small as R_out ≈ (6−9) AU. All disk models with R_out ≳ 25 AU would violate the 3σ CO J = 3 → 2 non-detection limit, independent of chemical details.

• The SED-fitting suggests that the dust grains in the surface of the protoplanetary disk of ET Cha (where the near-mid IR continuum forms) must be small in radius (sub-micron sized) and opaque in the optical and near-mid IR, i.e. absorption must dominate over scattering opacities. In the models, these spectral properties are provided by the inclusion of about 25% amorphous carbon, but other options like metallic iron are also possible.

• The disk gas mass of ET Cha is poorly constrained by our line observations CO J = 3 → 2 (non-detection), [OI] 63 μm, [OI] 6300 Å (LVC), and o-H₂v = 1 → 0 S(1) 2.122 μm. We find a variety of about equally well fitting disk models with total gas masses M_gas = (5 × 10^-5−3 × 10^-3) M_⊙. The forbidden lines of [OI] 6300 Å (LVC), and o-H₂ 2.122 μm are close to optically thin and hence quite useful for gas mass determination, but originate from hot gas in the inner disk regions only. [OI] 63 μm is massively optically thick and, hence, does not discriminate much between the various disk models for the case of ET Cha. We would need to observe additional far-IR or sub-mm spectral lines that originate from the outer disk parts to determine the gas mass with more confidence. However, in order to explain the [OI] 6300 Å (LVC) and o-H₂ 2.122 μm lines with a disk model, we require line optical depths τ_line ≳ 1 which translates into total vertical column densities in the inner disk regions of about N_O ≳ 10²¹   cm^-2 and N_H2 ≳ 10²⁴ cm^-2 for atomic oxygen and molecular hydrogen, respectively. Such column densities are incompatible with M_gas < 5 × 10^-5 M_⊙.

• The wide range of fitting gas masses is related to particular values of the disk shape parameters, namely the column density exponent ϵ and the disk flaring power β. High gas masses need to be combined with low values for ϵ and β, and vice versa.

• From our SED modelling of ET Cha we derive disk scale heights of the order of 10% relative to radius close to the star, which is about a factor of 2−3 larger than the scale heights inferred form self-consistent models that assume hydrostatic equilibrium. This discrepancy can partly be explained by an additional non-radiative heating of the dust close to the star, for example via exothermic chemical reactions.

These results suggest a surprisingly high value for the overall gas/dust mass ratio of ET Cha of at least 2000, or even 20 000. Whether or not the gas responsible for the [OI] 6300 Å and o-H₂ 2.122 μm emissions is still physically connected to the disk is not entirely clear, but the observations show that this gas is at least not moving much with respect to the disk. Possibly, the overwhelming majority of dust particles responsible for the near to far-IR emission of the star has already been transformed into bigger pebbles or solid bodies, which have negligible opacities.

The fluxes in blue-shifted emission line components like [OI] 6300 Å (HVC) suggest an outflow with mass-loss rate Ṁ_outflow ≈ 10^-9   M_⊙/yr, on a similar level as the reported mass accretion rate of ET Cha. The low velocities (≈−42 km s^-1) suggests a high inclination angle of at least 60° (rather more). However, inclinations in excess of 50° are inconsistent with our SED-modelling, which favours i ≲ 40°, from which we determine an outflow velocity of 65 km s^-1 at most (typical values are in excess of 100 km s^-1 for T Tauri outflows Hartigan et al. 1995). If this interpretation is correct, ET Cha would possess one of the slowest known outflows from a T Tauri star.

An outflow with mass loss rate Ṁ_outflow ≈ 10^-9   M_⊙/yr is likely to contribute significantly to the [OI] 63μm line flux as observed with Herschel/PACS. Given the observational data we have collected in this paper, we cannot distinguish between outflow or disk origin of [OI] 63 μm.

To conclude, ET Cha seems to be an extraordinary and puzzling object concerning the evolution of protoplanetary disks. Despite its age of about (6−8) Myr, there is evidence of active accretion onto the central star, on a similar level as the derived outflow mass-loss rate. According to the low disk masses we derive, the disk lifetime M_gas/Ṁ_acc ≈ (0.05−3) Myr is inconsistent with cluster age. If a generic gas/dust ratio of 100 was assumed, the disk lifetime, based on our relatively robust determination of the disk dust mass, would be even shorter, only  ~3000 yrs. Either the object is actually much younger than the age of the Chamaeleontis cluster or the object is going through a phase of unusually high mass accretion rate and outflow.

Online material

Appendix A: New features in ProDiMo

A.1. Fixed disk structure

In contrast to earlier publications, we use a parametrised description for the shape and distribution of gas and dust in the disk in this paper $\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{-}\mathit{ϵ}}\frac{{\mathit{H}}_{\mathrm{0}}}{\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\hspace{0.17em}\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathit{H}\mathrm{(}\mathit{r}{\mathrm{)}}^{\hspace{0.17em}\mathrm{2}}}\right)$(A.1)between an inner and outer disk radius, R_in and R_out, respectively, with sharp edges. ρ(r,z) is the local gas mass density. H(r) is the vertical scale height of the disk, assuming to vary with radius as $\mathit{H}\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 }}\mathrm{·}$(A.2)H₀ is the reference scale height at reference radius r₀. ϵ is the column density power-law index and β the flaring power. The constant ρ₀ is adjusted such that the integrated disk mass 2π ^! ρ(x,z) dzr dr equals M_disk.

A.2. Dust size distribution and dust settling

The dust grains are assumed to have a power-law size distribution in the unsettled case as $\mathit{f}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }{\mathit{a}}^{\hspace{0.17em}\mathrm{-}\mathit{p}}$(A.3)between minimum and maximum grain radius, a_min and a_max, respectively. The free constant in Eq. (A.3) is adjusted to result in the prescribed unsettled dust/gas mass ratio ρ_d/ρ.

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

A.3. PAH ionisation equilibrium and PAH-heating

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

Concerning the photo-ionisation reactions, we consider the PAH absorption cross sections ${\mathit{\sigma }}_{\mathrm{PAH}}^{\mathit{k}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\hspace{0.17em}\mathrm{\left[}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\mathrm{\right]}$ from Li & Draine (2001, see their Eqs. (6) −  with recent updates for the resonance parameters from Draine & Li (2007). We use their neutral PAH cross section for charge k = 0, and the charged cross section otherwise. The photo-ionisation rates  [s^-1]  for PAH molecules with charge k are calculated as ${\mathit{R}}_{\mathrm{ph}}^{\mathit{k}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathit{hc}}{\mathrm{\int }}_{\mathrm{912}\hspace{0.17em}\AA }^{\mathit{\lambda }\mathrm{thr}}{\mathit{\sigma }}_{\mathrm{PAH}}^{\mathit{k}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{\nu }{\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\hspace{0.17em}{\mathit{Y}}_{\mathit{\nu }}^{\mathit{k}}\hspace{0.17em}{\mathit{s}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{d}\mathit{\lambda }$(A.10)where ν [Hz] is the frequency, J_ν(r,z)  [erg/cm²/s/Hz/sr]  the mean intensity as computed by the dust continuum radiative transfer, and s_ν(r,z) the following self-shielding factor $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}& \mathrm{\approx }& \mathit{ϵ}\mathrm{\left(}\mathrm{PAH}\mathrm{\right)}\mathrm{⟨}\mathit{\sigma }{\mathrm{⟩}}_{\mathrm{PAH}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{Min}\mathrm{\{}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}^{\mathrm{ver}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathit{,}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}^{\mathrm{rad}}\mathrm{(}\mathit{r,z}{\mathrm{)}}^{\mathrm{\}}}\\ {\mathit{s}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}& \mathrm{=}& \exp \mathrm{(}{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}}^{\mathrm{)}}\mathit{.}\end{array}$${\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}^{\mathrm{ver}}$ and ${\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}^{\mathrm{rad}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}$ are the radial inward and vertical upward hydrogen nuclei column densities [cm^-2] in the disk, as measured from point (r,z). For simplicity, we put $\mathrm{⟨}\mathit{\sigma }{\mathrm{⟩}}_{\mathrm{PAH}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathrm{0.5}{}^{\mathrm{(}}\mathit{\sigma }_{\mathrm{PAH}}^{\hspace{0.17em}\mathrm{0}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{PAH}}^{\mathrm{+}\mathrm{1}}\mathrm{(}\mathit{\nu }{\mathrm{\left)}}^{\mathrm{\right)}}$. The photo-electron yield is taken from (Jochims et al. 1996) ${\mathit{Y}}_{\mathit{\nu }}^{\mathit{k}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1}\mathit{,}& \mathit{h\nu }\mathit{>}\mathit{I}{\mathit{P}}^{\hspace{0.17em}\mathit{k}}\mathrm{+}\mathrm{9.2}\hspace{0.17em}\mathrm{eV}\\ \frac{\mathit{h\nu }\mathrm{-}\mathit{I}{\mathit{P}}^{\hspace{0.17em}\mathit{k}}}{\mathrm{9.2}\hspace{0.17em}\mathrm{eV}}\mathit{,}& \mathit{I}{\mathit{P}}^{\hspace{0.17em}\mathit{k}}\mathrm{+}\mathrm{9.2}\hspace{0.17em}\mathrm{eV}\mathit{>}\mathit{h\nu }\mathit{>}\mathit{I}{\mathit{P}}^{\hspace{0.17em}\mathit{k}}\\ \mathrm{0}\mathit{,}& \mathit{h\nu }\mathit{<}\mathit{I}{\mathit{P}}^{\hspace{0.17em}\mathit{k}}\mathit{.}\end{array}$(A.13)According to Weingartner & Draine (2001), the ionisation potentials IP^k depend on PAH-size and charge k. Using their Eqs. (1) and (2) for N_C = 54, the results are IP⁻¹ = 3.10 eV, IP⁰ = 6.24 eV, IP^+ 1 = 9.38 eV and IP^+ 2 = 12.5 eV. The formula is an approximation that reproduces the ionisation potential of benzene as the smallest PAH and graphite as infinitely large PAH. The threshold wavelengths are given by ${\mathit{\lambda }}_{\mathrm{thr}}^{\mathit{k}}\mathrm{=}\mathit{hc}\mathit{/}\mathit{I}{\mathit{P}}^{\mathit{k}}$.

With these photo-ionisation and recombination rates, the local particle densities of PAH⁻, PAH, PAH⁺, PAH²⁺ and PAH³⁺ are calculated consistently with the gas-phase and ice chemistry, with a considerable influence on the resultant local electron density. Once these particle densities ${\mathit{n}}_{\mathrm{PAH}}^{\mathit{k}}$ have been determined, the total PAH heating rate  [erg/cm³/s]  can be calculated as ${\mathrm{\Gamma }}_{\mathrm{PAH}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathit{hc}}\sum _{\mathit{k}}{\mathit{n}}_{\mathrm{PAH}}^{\mathit{k}}\underset{\mathrm{912}\hspace{0.17em}\AA }{\overset{{\mathit{\lambda }}_{\mathrm{thr}}^{\mathit{k}}}{\mathrm{\int }}}{\mathit{\sigma }}_{\mathrm{PAH}}^{\mathit{k}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{\nu }{\mathit{J}}_{\mathit{\nu }}\hspace{0.17em}{\mathit{Y}}_{\mathit{\nu }}^{\mathit{k}}\hspace{0.17em}{\mathit{s}}_{\mathit{\nu }}\hspace{0.17em}(\mathit{h\nu }\mathrm{-}\mathit{I}{\mathit{P}}^{\hspace{0.17em}\mathit{k}})\mathrm{d}\mathit{\lambda }\mathit{.}$(A.14)The PAH recombination cooling rate is calculated as ${\mathrm{\Lambda }}_{\mathrm{PAH}}\mathrm{=}\sum _{\mathit{k}}{\mathit{n}}_{\mathrm{PAH}}^{\mathit{k}}\hspace{0.17em}{\mathit{n}}_{\mathrm{e}}\hspace{0.17em}{\mathit{k}}_{\mathrm{PAH}}^{\mathit{k}}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathrm{\right)}\hspace{0.17em}(\mathrm{1.5}\hspace{0.17em}\mathit{k}{\mathit{T}}_{\mathrm{g}})$(A.15)where ${\mathit{k}}_{\mathrm{PAH}}^{\mathit{k}}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathrm{\right)}\hspace{0.17em}\mathrm{[}\mathrm{c}{\mathrm{m}}^{\mathrm{3}}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{]}$ is the PAH recombination rate coefficient, n_e is the electron density and ${\mathit{T}}_{\mathrm{g}}$ the gas temperature. As a general result, we obtain a top-down layered PAH charge structure, with PAH²⁺ in the uppermost layers with reduced PAH-heating, but PAH⁻ and intensified PAH-heating close to the midplane. For massive disks (unlike ET Cha), where the electron density virtually vanishes in the deepest layers, the PAH charge becomes neutral again in the midplane.

A.4. Fluorescent UV-pumping

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

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

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

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

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

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

Our formulation of the H₂-pumping, therefore, leads to a constant de-population of all states l < u with a certain timescale ${\mathit{C}}_{\mathit{l}\hspace{0.17em}\mathit{u}}^{\mathrm{chem}}\hspace{0.17em}\mathrm{\left[}{\mathrm{s}}^{-1}\mathrm{\right]}$ (simulating the disappearance of H₂ in all states due to other chemical channels that are feeding the H₂ formation) and an equally strong population  [cm^-3 s^-1]  of the state u due to the formation of H₂ on grains.

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

A.7. Line transfer

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

The continuum  +  line radiative transfer equation is given by $\frac{\mathrm{d}{\mathit{I}}_{\mathit{\nu }}}{\mathrm{d}\mathit{s}}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{ext}}({\mathrm{S}}_{\mathit{\nu }}\mathrm{-}{\mathit{I}}_{\mathit{\nu }})\mathit{,}$(A.20)where I_ν is the spectral intensity and s the distance on a ray. The source function and extinction coefficient are given by $\begin{array}{ccc}{\mathrm{S}}_{\mathit{\nu }}& \mathrm{=}& \frac{{\mathit{ϵ}}_{\mathit{\nu }}^{\mathrm{D}}\mathrm{+}{\mathit{\phi }}_{\mathit{\nu }}{\mathit{ϵ}}_{\mathit{\nu }}^{\mathrm{L}}}{{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{ext}}}\\ {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{ext}}& \mathrm{=}& {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{D}}\mathrm{+}{\mathit{\phi }}_{\mathit{\nu }}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{L}}\mathit{.}\end{array}$Assuming isotropic dust scattering, the continuum and line transfer coefficients are given by $\begin{array}{ccc}{\mathit{ϵ}}_{\mathit{\nu }}^{\mathrm{D}}& \mathrm{=}& {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\mathrm{+}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{sca}}{\mathit{J}}_{\mathit{\nu }}\\ {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{D}}& \mathrm{=}& {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}\mathrm{+}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{sca}}\\ {\mathit{ϵ}}_{\mathit{\nu }}^{\mathrm{L}}& \mathrm{=}& \frac{\mathit{h\nu }}{\mathrm{4}\mathit{\pi }}{\mathit{n}}_{\mathit{u}}{\mathit{A}}_{\mathit{ul}}\\ {\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{L}}& \mathrm{=}& \frac{\mathit{h\nu }}{\mathrm{4}\mathit{\pi }}({\mathit{n}}_{\mathit{l}}{\mathit{B}}_{\mathit{lu}}\mathrm{-}{\mathit{n}}_{\mathit{u}}{\mathit{B}}_{\mathit{ul}})\mathit{,}\end{array}$where ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{abs}}$ and ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{sca}}\hspace{0.17em}\mathrm{\left[}\mathrm{c}{\mathrm{m}}^{-1}\mathrm{\right]}$ are the local dust absorption and scattering coefficients, B_ν is the Planck function, J_ν is the local mean intensity (taken from the results of the continuum radiative transfer), ν is the line centre frequency, n_u and n_l are the level populations  [cm^-3]  in the upper and lower state, respectively, and A_ul and B_ul,B_lu are the Einstein coefficients.

The profile function φ_ν  [Hz^-1]  is assumed to be given by a Gaussian with thermal  +  turbulent broadening $\begin{array}{ccc}\mathit{x}& \mathrm{=}& {\mathit{v}}^{\mathrm{\prime }}\mathrm{+}{n}\mathrm{·}{v}\\ {\mathit{\phi }}_{\mathit{\nu }}& \mathrm{=}& \frac{\mathit{c}}{\mathit{\nu }\sqrt{\mathit{\pi }}\mathrm{\Delta }\mathit{v}}\exp \left(\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{\Delta }{\mathit{v}}^{\mathrm{2}}}\right)\mathit{,}\end{array}$where v′ is the observers velocity with respect to the star along backward ray direction n, v is the 3D-velocity of the emitting gas with respect to the star (assumed to be given by Keplerian orbits), and x the local velocity shift. The velocity width is $\mathrm{\Delta }{\mathit{v}}^{\mathrm{2}}\mathrm{=}{\mathit{v}}_{\mathrm{th}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{turb}}^{\mathrm{2}}$ and the thermal velocity ${\mathit{v}}_{\mathrm{th}}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathit{k}{\mathit{T}}_{\mathrm{g}}\mathit{/}\mathit{m}$ where m is the mass of the line emitting species.

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

A.8. Chemical heating

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

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

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

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

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

Appendix B: The outflow model

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

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

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

Appendix C: Details of the best-fitting disk model

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

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

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

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

Appendix D: Variation of input physics

D.1. Difference between dust and gas temperature

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

D.2. Influence of X-rays

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

D.3. Influence of viscous heating

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

D.4. Self-consistent disk structure

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

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

D.5. Influence of non-radiative dust heating

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

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

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

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

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

Concerning the formulation of Sternberg & Dalgarno (1995), which is valid for standard ISM size distribution and dust/gas ratio only, we add a scaling factor to account for deviations of the total dust surface per hydrogen nucleus in disks as ${\mathit{R}}_{\mathrm{H}\mathrm{2}}^{\hspace{0.17em}\mathrm{form}}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-18}{\mathit{n}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}\sqrt{{\mathit{T}}_{\mathrm{gas}}}\frac{\mathrm{⟨}{\mathit{a}}^{\mathrm{2}}\mathrm{⟩}\hspace{0.17em}{\mathit{n}}_{\mathrm{dust}}\mathit{/}{\mathit{n}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}}{\mathrm{5.899}\mathrm{\times }{\mathrm{10}}^{-22}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}}\mathrm{·}$(D.3)

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

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

---

Acknowledgments

We thank Catherine Dougados for fruitful discussions about the properties of disk outflows and the interpretation of optical emission lines. W.-F. Thi acknowledges a SUPA astrobiology fellowship. I. Pascucci and B. Riaz acknowledge NASA/JPL for funding support. The LAOG group acknowledges PNPS, CNES and ANR (contract ANR-07-BLAN-0221) for financial support. I. de Gregorio-Monsalvo is partially supported by Ministerio de Ciencia e Innovación (Spain), grant AYA 2008-06189-C03 (including FEDER funds), and by Consejería de Innovación y Ciencia y Empresa of Junta de Andalucía, (Spain).

References

All Tables

All Figures