# Gas modelling in the disc of HD 163296

I. Tilling, P. Woitke, G. Meeus, A. Mora, B. Montesinos, P.

Riviere-Marichalar, C. Eiroa, W. -F. Thi, A. Isella, A. Roberge, et al.

► To cite this version:

I. Tilling, P. Woitke, G. Meeus, A. Mora, B. Montesinos, et al. Gas modelling in the disc of HD 163296. Astronomy and Astrophysics - A&A, 2012, 538, 10.1051/0004-6361/201116919. insu-03612453

HAL Id: insu-03612453

https://insu.hal.science/insu-03612453

Submitted on 17 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gas modelling in the disc of HD 163296

I. Tilling¹, P. Woitke²,³,⁴, G. Meeus⁵, A. Mora⁶, B. Montesinos⁷, P. Riviere-Marichalar⁷, C. Eiroa⁵, W.-F. Thi⁸, A. Isella⁹, A. Roberge¹⁰, C. Martin-Zaidi⁸, I. Kamp¹¹, C. Pinte⁸, G. Sandell¹², W. D. Vacca¹², F. Ménard⁸, I. Mendigutía⁷, G. Duchêne⁸,¹³, W. R. F. Dent¹⁴, G. Aresu¹¹, R. Meijerink¹⁵, and M. Spaans¹¹

¹ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
e-mail: it@roe.ac.uk

² University of Vienna, Dept. of Astronomy, Türkenschanzstr. 17, 1180 Vienna, Austria

³ UK Astronomy Technology Centre, Royal Observatory, Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK

⁴ SUPA, School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK

⁵ Dep. de Física Teórica, Fac. de Ciencias, UAM Campus Cantoblanco, 28049 Madrid, Spain

⁶ ESA-ESAC Gaia SOC, PO Box 78. 28691 Villanueva de la Cañada, Madrid, Spain

⁷ Departamento de Astrofísica, Centro de Astrobiología (CAB, INTA-CSIC), ESAC Campus, PO Box 78, 28691 Villanueva de la Cañada, Madrid, Spain

⁸ UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique (IPAG) UMR 5274, 38041 Grenoble, France

⁹ Department of Astronomy, California Institute of Technology, MC 249-17, Pasadena, CA 91125, USA

¹⁰ NASA Goddard Space Flight Center, Exoplanets & Stellar Astrophysics, Code 667, Greenbelt, MD 20771, USA

¹¹ Kapteyn Astronomical Institute, Postbus 800, 9700 AV Groningen, The Netherlands

¹² SOFIA-USRA, NASA Ames Research Center, Mailstop 211-3 Moffett Field CA 94035, USA

¹³ Astronomy Department, University of California, Berkeley, CA 94720-3411, USA

¹⁴ ESO-ALMA, Avda Apoquino 3846, Piso 19, Edificio Alsacia, Las Condes, Santiago, Chile

¹⁵ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands

Received 18 March 2011 / Accepted 2 November 2011

## ABSTRACT

We present detailed model fits to observations of the disc around the Herbig Ae star HD 163296. This well-studied object has an age of ~4 Myr, with evidence of a circumstellar disc extending out to ~540 AU. We use the radiation thermo-chemical disc code ProDiMo to model the gas and dust in the circumstellar disc of HD 163296, and attempt to determine the disc properties by fitting to observational line and continuum data. These include new Herschel/PACS observations obtained as part of the open-time key program GASPS (GAS in Protoplanetary Systems), consisting of a detection of the [OI] 63 µm line and upper limits for several other far infrared lines. We complement this with continuum data and ground-based observations of the ¹²CO 3–2, 2–1 and ¹³CO J = 1–0 line transitions, as well as an upper limit for the H₂ 0–0 S(1) transition. We explore the effects of stellar ultraviolet variability and dust settling on the line emission, and on the derived disc properties. Our fitting efforts lead to derived gas/dust ratios in the range 9–100, depending on the assumptions made. We note that the line fluxes are sensitive in general to the degree of dust settling in the disc, with an increase in line flux for settled models. This is most pronounced in lines which are formed in the warm gas in the inner disc, but the low excitation molecular lines are also affected. This has serious implications for attempts to derive the disc gas mass from line observations. We derive fractional PAH abundances between 0.007 and 0.04 relative to ISM levels. Using a stellar and UV excess input spectrum based on a detailed analysis of observations, we find that the all observations are consistent with the previously assumed disc geometry.

**Key words.** circumstellar matter – astrochemistry – stars: pre-main sequence – protoplanetary disks – stars: fundamental parameters – line: formation

## 1. Introduction

Protoplanetary discs set the initial conditions for planet formation, and the processing of chemical species and dust grains as discs evolve may ultimately provide the seed for life. Indeed, conditions within circumstellar discs are conducive to a rich chemistry, with a large variety of conditions ranging from icy mantles on grain surfaces in the cold disc midplane to ionised gas bathed in stellar ultraviolet radiation at the disc surface. In order to attempt to probe the gas in discs it is necessary to observe gas emission lines from a variety of chemical species (molecules, atomic species, ions), and couple these observations with detailed thermochemical disc modelling.

One important disc property for planet formation models is the amount of gas present as it evolves. The dust characteristics are generally better-understood than the gas, and

the ISM gas/dust mass ratio of 100 is often adopted to give an estimate of the gas mass. However, gas masses derived from millimetre and sub-mm CO emission observations are typically lower than from dust observations with this assumption, sometimes by a factor ~100 (Thi et al. 2001). This is variously attributed to gas depletion due to photoevaporation or planet formation, and CO freeze-out on grains (Zuckerman et al. 1995; Kamp & Bertoldi 2000). These low level rotational CO transitions are also strongly dependent on the disc size, and often optically thick under certain disc conditions, limiting their ability to probe the deeper disc layers (Panić & Hogerheijde 2009; Hughes et al. 2008). In order to better constrain the disc gas mass it is necessary to observe additional gas species, such as those observable in the far-infrared. The FIR range is a crucial observing window which can help to resolve some of these ambiguities. The fine structure lines of atomic oxygen and ionised carbon
<figure><img src="image_2.png" /><figcaption>Fig. 1. Observed [OⅠ]63 µm emission line with Herschel/PACS.</figcaption></figure>

probe the gas in the warm disc surface layers, and as products of CO photodissociation provide a good test for this proposed ex- planation of CO under-abundance. In addition, the high level ro- tational transitions of CO arise from warm gas in the inner disc, probing different disc regions to the low level mm CO lines.

In addition to the discrepancy between gas mass estimates from millimetre continuum and CO observations, there is uncer- tainty regarding the derived disc outer radii from CO and con- tinuum observations. The outer radius suggested by model fits to dust emission is often found to be smaller than that of the gas disc. In the case of the Herbig Ae star AB Aurigae, Piétu et al. (2005) suggested that a change in dust grain properties, leading to reduced opacity in the outer disc, could explain this discrepancy. For HD 163296, Isella et al. (2007) attributed the apparent difference in disc radii to a sharp drop in surface den- sity, opacity or dust-to-gas ratio beyond 200 AU, while Hughes et al. (2008) proposed that the apparently conflicting dust and gas observations could be reconciled using a different treatment of the disc surface density at the outer edge, motivated by sim- ilarity solutions for the time evolution of accreting discs. This study attributed smaller derived disc radii from dust observations to detector sensitivity thresholds, arguing in favour of larger dust discs, albeit with an exponentially-tapered density at the outer edge. A steepening of the density profile in the outer disc has also been observed in the case of DM Tau, LkCa 15 and MWC 480 by Piétu et al. (2007).

The Herschel open-time key program GASPS (GAS in Protoplanetary Systems) aims to characterise the gas in discs at each stage of their evolution, by observing the far-IR lines of a wide range of discs, from young gas-rich accretion discs to gas-poor debris discs (Mathews et al. 2010). This paper rep- resents the second in-depth study of a Herbig disc as part of this program, following after a study focusing on the young Herbig Ae star HD 169142 (Meeus et al. 2010).

HD 163296 is an isolated Herbig Ae star with spectral type A1Ve, mass ~2.5 $M_\odot$ and stellar luminosity ~38 $L_\odot$ (this work), situated at a distance of 118.6$^{+12.7}_{-10.0}$ pc (van Leeuwen 2007). Scattered light and millimetre continuum observations indicate the presence of an inclined circumstellar disc extend- ing out to a radius of 540 AU (Mannings & Sargent 1997; Grady et al. 2000; Wisniewski et al. 2008). In addition, there is evi- dence of an asymmetric outflow perpendicular to the disc, with a chain of six Herbig-Haro knots tracing its mass-loss history (Devine et al. 2000). The derived disc dust mass is in the range (5–17) × 10$^{-4} M_\odot$ (Natta et al. 2004; Tannirkulam et al. 2008a; Mannings & Sargent 1997; Isella et al. 2007).

Recent near-infrared studies of the inner disc of HD 163296 indicate an inner dust rim feature enclosing a bright emission region extending inwards towards the star. This emission in- wards of the dust rim has been attributed to an optically thin inner disc, although there is uncertainty regarding the size and composition of such a disc. Derived radii for the dust rim lie in the range (0.2–0.55) AU (Renard et al. 2010; Eisner et al. 2009; Tannirkulam et al. 2008b; Benisty et al. 2010; Monnier et al. 2006). An increase in opacity at the inner rim due to sublimation of grains could cause the rim to puff up, and it has been sug- gested that this leads to time-variable self-shadowing of the outer disc (Sitko et al. 2008; Wisniewski et al. 2008). Mid-infrared imaging of warm dust in the disc surface layers seems to indicate little flaring (Doucet et al. 2006), and scattered light brightness profiles are consistent with a non-flared outer disc (Wisniewski et al. 2008). Fits to the ISO-SWS spectra of HAeBe discs led Meeus et al. (2001) to classify HD 163296 as one of their group II objects, in which an optically thick inner disc shadows the outer disc from stellar radiation, resulting in a non-flared geom- etry.

In this study we aim to fit a wide variety of observed emis- sion lines and the dust spectral energy distribution (SED) simul- taneously with a disc with this observed geometry.

## 2. Herschel/PACS observations

We have obtained observations of HD 163296 in the far- infrared using Herschel/PACS, consisting of spectroscopic line observations and photometry. These include a detection of the [OⅠ]63 µm fine structure line, and upper limits for a number of other atomic, ionic and molecular far-IR lines. We observed HD 163296 with the PACS spectrometer in line spectroscopy mode (obsid 1342192161, 1252 s, 3 repetitions) and range scan mode (obsid 1342192160, 5141 s, 3 repetitions), while chopping and nodding.

The spectra were reduced using HIPE version 7.0, with the standard pipeline scripts for a chopped line scan of a point source. The reduction steps include division by the spectral re- sponse function and MedianOffsetCorrection, as the new flat fielding task is not yet robust for short range scans such as the present observations (PACS instrument team, priv. comm.). For this paper, we only make use of the data contained in the central spaxel, to which we apply a wavelength-dependent correction factor for the loss of flux due to diffraction losses. With these additional steps, the absolute flux calibration ac- curacy is 30% (see PACS Spectroscopy Performance and Calibration Document, http://herschel.esac.esa.int/ AOTsReleaseStatus.shtml).

The line fluxes and upper limits are given with continuum data in Table 1. The PACS photometry data are given in Table 2. Comparison of the azimuthally averaged radial profiles with the observed point spread function (PSF), using Vesta, show that the disc is unresolved at both 70 µm and 160 µm. The spatial ex- tent of the observed [OⅠ]63 µm emission is illustrated in Fig. 2, where it can be seen that the emission is dominated by the central 9.4 × 9.4 arcsec spaxel. This corresponds to a maximum radial extent of 560 AU for the emission.

The pipeline reduction of the [Cn] 158 µm range shows this line in absorption, which is unexpected. Therefore, we also looked at the unpchopped spectra in the different chop positions. It is clear that the region around HD 163296 is full of [Cn] 158 µm emission, and that this causes the emission feature at the position of the star to cancel out (and even become an absorption feature). This is likely due to cloud material along the line-of-sight, since
<table><caption>**Table 1.** Observed line data with Herschel/PACS.</caption><thead><tr><th>Species</th><th>λ [μm]</th><th>ν [GHz]</th><th>Line flux [10<sup>-18</sup> W m<sup>-2</sup>]</th><th>Continuum (rms) [Jy]</th></tr></thead><tbody><tr><td>Oi</td><td>63.18</td><td>4745.05</td><td>193.1 ± 5.8</td><td>16.46 ± 0.07</td></tr><tr><td>Oi</td><td>145.52</td><td>2060.15</td><td>&lt;8.5</td><td>22.12 ± 0.03</td></tr><tr><td>Cii</td><td>157.74</td><td>1900.55</td><td>–</td><td>23.90 ± 0.05</td></tr><tr><td>p-H<sub>2</sub>O</td><td>89.99</td><td>3331.40</td><td>&lt;9.4</td><td>19.14 ± 0.04</td></tr><tr><td>o-H<sub>2</sub>O</td><td>179.52</td><td>1669.97</td><td>&lt;14.5</td><td>22.28 ± 0.08</td></tr><tr><td>o-H<sub>2</sub>O</td><td>180.42</td><td>1661.64</td><td>&lt;16.2</td><td>21.58 ± 0.09</td></tr><tr><td>o-H<sub>2</sub>O</td><td>78.74</td><td>3810.01</td><td>&lt;15.0</td><td>21.00 ± 0.10</td></tr><tr><td>OH</td><td>79.11</td><td>3792.19</td><td>&lt;17.0</td><td>20.70 ± 0.10</td></tr><tr><td>OH</td><td>79.18</td><td>3788.84</td><td>&lt;17.0</td><td>20.60 ± 0.10</td></tr><tr><td>CO J = 36–35</td><td>72.85</td><td>4115.20</td><td>&lt;11.6</td><td>17.30 ± 0.03</td></tr><tr><td>CO J = 33–32</td><td>79.36</td><td>3777.63</td><td>&lt;22.8</td><td>18.12 ± 0.07</td></tr><tr><td>CO J = 29–28</td><td>90.16</td><td>3325.12</td><td>&lt;11.1</td><td>19.10 ± 0.05</td></tr><tr><td>CO J = 18–17</td><td>144.78</td><td>2070.68</td><td>&lt;13.1</td><td>22.05 ± 0.04</td></tr></tbody></table>

**Notes.** Detections are listed as $F_{i} \pm \sigma$) whereas non-detections quote $<3\sigma$. There is an additional absolute flux calibration error of 30%. In the case of the [Cii] 158 µm line, strong contamination from the surrounding background prevented us from estimating a meaningful upper limit.

<table><caption>**Table 2.** Photometric fluxes observed with Herschel/PACS.</caption><thead><tr><th>λ [μm]</th><th>70</th><th>160</th></tr></thead><tbody><tr><td>Flux [Jy]</td><td>18.5 ± 3.3</td><td>20.2 ± 4.0</td></tr></tbody></table>

**Notes.** These observations (AORS 1342228401 and 1342228401) are done with the standard miniscan.

HD 163296 is located close to the Galactic plane, and is consistent with the presence of CO emission at radial velocities offset from this object, but unresolvable with PACS (Thi et al. 2001). We also inspected all the other lines for their unchopped appearance, but did not find contamination through the offset positions in those cases.

### 3. Properties of HD 163296

The study of the disc of HD 163296 is directly linked to the knowledge of the absolute parameters of the star itself and the behaviour of the SED in the ultraviolet (UV), optical and near infrared (near-IR), which is mostly dominated by stellar radiation; in particular, the energy emitted at wavelengths shorter than 2500 Å is a fundamental piece of information used to model the disc, due to its impact on the photo-chemistry and gas heating (Woitke et al. 2010). This section describes the stellar properties and studies the SED at the range mentioned above.

#### 3.1. Stellar parameters

The determination of the stellar parameters of HD 163296 was done in an extensive work devoted to the study of a sample of Herbig Ae/Be stars (Montesinos et al. 2009). Details of the methodology followed and the observations used can be found in that paper. We outline in this section the relevant steps of the process.

An iterative distance-independent algorithm based on the analysis of the SED and optical high-resolution spectra (for the estimation of the effective temperature and metallicity), and mid-resolution spectra of the Balmer lines Hβ, Hγ and Hδ (for the estimation of the gravity), was applied.

The knowledge of the effective temperature, gravity and metal abundance allows us to place the star in a $\log g_{*} - \log T_{\text{eff}}$ diagram and superimpose the appropriate set of evolutionary tracks and isochrones for that specific metallicity (Yi et al. 2001).

This gives directly – or after a simple interpolation between the tracks and isochrones enclosing the position of the star – the stellar mass and the age. Since there is a one-to-one correspondence between a pair ($T_{\text{eff}}$, $g_{*}$) on a given track and a pair ($T_{\text{eff}}$, $L_{*}/L_{\odot}$), the stellar luminosity can also be estimated.

Table 3 gives the stellar parameters of HD 163296 after a slight refinement of the results presented in Montesinos et al. (2009).

#### 3.2. The SED: UV, optical, and near-IR observations

Simultaneous optical *UBVRI* and near-infrared *JHK* photometry obtained as part of the EXPORT collaboration (Eiroa et al. 2001; Oudmaijer et al. 2001), were used to build the optical and near-IR photospheric SED of HD 163296. Table 4 shows the wavelengths, magnitudes and corresponding fluxes with their uncertainties. The calibration from magnitudes to fluxes has been done using the zero points given by Bessell (1979) for the optical and Cox (1999) for the near-IR magnitudes.

Ultraviolet IUE spectra of HD 163296 were extracted from the INES archive<fnref n="1" />. From the collection of 68 SW and LW (for “short” and “long wavelength”) spectra in low resolution, obtained through the large aperture of the spectrograph, all those that looked clean (few or none bad pixels) and with the correct exposure classification codes were selected. The spectra (SW + LW) cover the range 1250–3000 Å.

A high-resolution, far-UV spectrum of HD 163296 was constructed by combining data taken with HST STIS and FUSE. The STIS spectra, covering 1150 Å to 3000 Å, were obtained from the HST STIS Echelle Spectral Catalog of Stars (StarCat<fnref n="2" />). The FUSE spectrum of HD 163296 is extremely noisy below 970 Å; therefore, we used only the portion between that wavelength and 1150 Å in our high-resolution composite UV spectrum for HD 163296. FUSE spectra contain strong terrestrial airglow emission lines, due to the large size of the observing aperture. We removed these lines from the HD 163296 spectrum before stitching it to the STIS spectrum.

A comparison of the IUE and FUSE+STIS spectra shows that they differ in the level of their continuum. The IUE spectra with the largest flux levels are a factor ~1.6 higher than the FUSE+STIS spectrum, therefore, regarding the modelling, the

<footnote n="1">¹ http://sdc.cab.inta-csic.es/ines/</footnote>

<footnote n="2">² http://casa.colorado.edu/~ayres/StarCAT/</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** [OⅠ]63 µm spectra observed by Herschel/PACS towards HD 163296. Each spectrum corresponds to a 9.4 × 9.4 arcsec pixel centered at the given coordinates. Emission is dominated by the central pixel, corresponding to a maximum outer radius of ~560 AU, assuming a distance of 118.6 parsecs to this object.</figcaption></figure>

<table><caption>**Table 3.** Parameters for HD 163296.</caption><tr><td>Temperature</td><td>9250 ± 150 K</td></tr><tr><td>Gravity (log g<sub>*</sub>)</td><td>4.07 ± 0.09</td></tr><tr><td>Metallicity ([Fe/H])</td><td>+0.20 ± 0.10</td></tr><tr><td>Mass</td><td>2.47 ± 0.10 M<sub>☉</sub></td></tr><tr><td>Luminosity</td><td>37.7 ± 7.0 L<sub>☉</sub></td></tr><tr><td>Age</td><td>4.2 ± 0.4 Myr</td></tr></table>

ultraviolet contribution of the star will be represented by two different UV input spectra, namely a “low UV” state, corresponding to the FUSE+STIS spectrum, and a “high UV” state, corresponding to the average of the highest IUE spectra<fnref n="3" /> (see Sect. 4.1).

In addition, a near-infrared spectrum of HD 163296 was obtained at the NASA InfraRed Telescope Facility (IRTF) on Mauna Kea on 2010 July 8 with SpeX (Rayner et al. 2003). Ten individual exposures of the target, each lasting 32 s, were taken using the short-wavelength cross-dispersed (SXD) mode of SpeX. This mode yields spectra spanning the wavelength range 0.8–2.4 µm divided into six spectral orders. The slit width was set to 0.′3, which yields a nominal resolving power of 2000

<table><caption>**Table 4.** Photometry for HD 163296.</caption><thead><tr><th>Band</th><th>λ (μm)</th><th>Magnitude and error</th><th>Flux and error (Jy)</th></tr></thead><tbody><tr><td>U</td><td>0.360</td><td>6.95 ± 0.03</td><td>3.00 ± 0.08</td></tr><tr><td>B</td><td>0.440</td><td>6.92 ± 0.04</td><td>7.27 ± 0.27</td></tr><tr><td>V</td><td>0.550</td><td>6.86 ± 0.03</td><td>6.56 ± 0.18</td></tr><tr><td>R</td><td>0.640</td><td>6.73 ± 0.05</td><td>6.26 ± 0.29</td></tr><tr><td>I</td><td>0.790</td><td>6.67 ± 0.10</td><td>5.48 ± 0.51</td></tr><tr><td>J</td><td>1.215</td><td>6.17 ± 0.05</td><td>5.55 ± 0.26</td></tr><tr><td>H</td><td>1.654</td><td>5.49 ± 0.05</td><td>6.69 ± 0.31</td></tr><tr><td>K</td><td>2.179</td><td>4.71 ± 0.05</td><td>8.56 ± 0.39</td></tr></tbody></table>

for the SXD spectra. Observations of HD 156717, an A0 V star, used as a “telluric standard” to correct for absorption due to the Earth’s atmosphere and to flux calibrate the target spectra, were obtained immediately preceding the observations of HD 163296. The airmass difference between the observations of the object and the standard was 0.06. The seeing was estimated to be ~0.′3–0.′5 at 2.2 µm during the observations and conditions were clear. The statistical S/N varies across the spectral range but is of the order of several hundred over the entire SXD spectrum. The SpeX data were reduced with Spextool (Cushing et al. 2004) and the telluric correction was performed using the procedures and software described by Vacca et al. (2003).

<footnote n="3">³ The average IUE spectrum was built with the following individual spectra: SWP 28805, SWP 28811, SWP 28813 and LWR 02065.</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Spectral energy distribution for HD 163296. The green dots represent the fluxes corresponding to the *UBVRIJHK* photometry. The size of the dots is of the order or larger than the uncertainties. The IUE average spectrum and the FUSE+STIS spectrum are plotted in blue and red, respectively; the latter has been slightly smoothed to reduce the noise. The SpeX/IRTF spectrum is plotted in magenta. The black solid line is a photospheric model computed for the specific stellar parameters given in Table 3, reddened with $E(B-V) = +0.15$ and normalized at the flux in V. See text for details.

Figure 3 shows the SED of HD 163296. The solid black line is a photospheric model computed for $T_{\text{eff}} = 9250$ K, $\log g_* = 4.07$ and [Fe/H] = +0.20 from the GAIA grid of models created with the PHOENIX code (Brott & Hauschildt 2005). The model has been reddened with $E(B-V) = +0.15$ ($R_V = 3.1$); this will be discussed in Sect. 3.4.

### 3.3. Spectral type from the ultraviolet and near-infrared spectra

Despite the fact that the precise determination of the stellar parameters has been done through a careful analysis of a variety of observations, as we have mentioned in Sect. 3.1, an attempt to determine the spectral type of HD 163296 in a more qualitative way has been done using the ultraviolet and near-infrared spectra.

A comparison of the ultraviolet data with IUE spectra of stars classified as “spectral-type standards” was done. The spectra were taken from the “IUE Ultraviolet Spectral Atlas of Standard Stars”<fnref n="4" /> (Wu et al. 1983, 1992). It is interesting to note that below 2000 Å, for a given spectral type, the ultraviolet spectra of the “spectral-type standards” themselves are very different, in particular for A0 V and A1 V. If we consider only the region above 2000 Å, the comparison suggests that HD 163296 is closer to A3 V (and possibly slightly later), however, if a mild extinction is present, the UV spectra would be also close to that of A1 V stars. A larger extinction would push the spectral type determination towards earlier spectral types.

A similar exercise of spectral typing was done by a comparison of the SpeX spectrum in the *YJ* bands (0.85–1.35 µm) with Kurucz synthetic models (see Montesinos et al. (2009) for references regarding the synthesis codes and models used). The spectrum is complex, the range of $T_{\text{eff}}$’s explored with the Kurucz models (8750–9250 K, corresponding to A3V and A1V) does

not produce significant Paschen line variations. The Paschen line wings of HD 163296 are narrower than in any of the models, which would suggest temperatures of 9250 K or higher. The metallic lines are faint, their typical intensity being lower than 5% that of the continuum. Some lines are deeper than in any photospheric spectrum, which suggests the presence of some circumstellar absorption. Some lines are much fainter, which suggests an spectral type A1V or earlier, or dust emission veiling. Good agreement is rarely found. Although we fail to find a unique combination of $T_{\text{eff}}$ and $\log g$ that provides an adequate fit to both the UV and NIR spectra of HD 163296, all data are consistent with an A1-A3 star, seen through mild foreground extinction.

### 3.4. The extinction

The determination of the extinction is an intricate problem. The classical method of comparing the observed colours with those of standard stars of similar properties does not lead to any conclusive result, given the slight variability of the star and the photospheric (and spectral type) uncertainties. For normal stars, the UV part of the SED can be used to estimate this parameter, but in the case of HD 163296 the ultraviolet spectrum is variable and has contributions both from the photosphere and the accretion shock that are difficult to quantify. The extinction seems to be mostly circumstellar, therefore the exact value of the constant $R_V$ is not known. The observations available do not allow us to disentangle completely this problem.

The strategy followed to assign a value of $E(B-V)$ to this star has four steps: 1) introduce extinctions from $E(B-V) = 0.0$ to 0.2 to the photospheric model; 2) normalize the reddened model photosphere to the flux at V; 3) deredden the normalized model and estimate the integrated stellar flux, $F_*$; note that $F_*$ is the flux that would be observed in the complete absence of extinction. Obviously, different extinctions imply different values of $F_*$. A value of $R_V = 3.1$ has been used.

The fourth step is based on the following argument: since the stellar luminosity $L_*/L_\odot$ is known through the distance-independent method outlined in Sect. 3.1, an estimation of the distance, *d*, can be obtained simply from the expression $L_* = 4\pi d^2 F_*$, where it is assumed that the whole star is visible, i.e. and no part of the stellar flux is completely blocked by the disc of other circumstellar material (see a complete discussion of this issue in Montesinos et al. 2009). The value $E(B-V) = +0.15$ allows us to recover a distance of 128 pc, which matches the Hipparcos distance to within the uncertainties, therefore it has been adopted as the optimum extinction for consistency with the whole set of parameters. We would like to point out that that value is a simple parameterization of a very complex problem and that a deeper study, outside the scope of this paper, would be needed to study all its peculiarities.

In Fig. 3 the photospheric model reddened with $E(B-V) = +0.15$ has been plotted. It goes through the two ultraviolet spectra but falls below both of them at wavelengths shorter than ~1260 Å, and it also underestimates the flux at U. Note that the model is purely photospheric, therefore the excess fluxes inferred from the observed data when compared with the model might be due to emission from the accretion shock, i.e. with a non-photospheric origin. The observed SED seems to depart from the photospheric SED at around ~1 µm, where the contribution from the disc starts to become noticeable.

<footnote n="4">⁴ http://www-int.stsci.edu/~jinger/iweb/proj/project.html</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Spectral energy distribution for the best-fit fully mixed disc model ("preferred" model), obtained from a simultaneous fit to the observed SED, ISO-SWS spectrum and line emission from HD 163296 (solid black line). Black dotted line indicates the SED for the model with power-law density profile, which is unable to fit the mm continuum image for HD 163296. Blue circles indicate (with increasing wavelength) simultaneous *UBVRIJHK* photometry (Eiroa et al. 2001; Oudmaijer et al. 2001), LM photometry (de Winter et al. 2001), IRAS photometry (12–100 mic), sub-mm photometry (Mannings 1994), and millimetre photometry (Isella et al. 2007). Also marked are PACS photometric observations (red circles), PACS continua derived from the spectroscopic observations (black circles), SCUBA photometry (green circles) (Sandell et al. 2011), scaled VLT/UVES spectrum (blue line, Martin-Zaidi, in prep.) and the ISO-SWS spectrum (green line). The red line shows the stellar+UV input spectrum. Downwards arrow denotes upper limit. All fluxes were corrected for interstellar reddening using the Fitzpatrick parameterisation (Fitzpatrick 1999) with $R_V = 3.1$ and $E(B-V) = 0.15$.</figcaption></figure>

### 3.5. Herschel observational constraints

In addition to the new *Herschel* data we have further constrained our disc models using SMA interferometric observations of the *J* = 3–2 transition in $^{12}$CO and PBI interferometry of the $^{12}$CO 2–1 and $^{13}$CO 1–0 transitions (Isella et al. 2007), an upper limit for the S(1) transition in H$_2$ obtained by VLT/VISIR (Martin-Zaidi et al. 2010), the ISO-SWS spectrum in the infrared, and a wealth of photometric data, including SCUBA photometry in the sub-mm (Sandell et al. 2011). See Fig. 4 for a full list of the photometry sources.

There is a wealth of observational data available for HD 163296 which constrains to varying extents our choice of model parameters. As well as fitting to the data mentioned in the previous paragraph, we have taken care not to contradict current understanding of the disc geometry. There is at present some uncertainty regarding the radial extent of the disc. Scattered light imaging indicates an outer edge for the dusty disc of ~540 AU, consistent with that derived for the gas disc from millimetre CO emission by Isella et al. (2007). However, the resolved mm continuum emission seems to indicate a smaller outer disc radius ~200 AU (Isella et al. 2007). Hughes et al. (2008) propose that the millimetre continuum and CO 3–2 emission can be fit simultaneously by a disc with an exponentially-tapered outer edge.

<table><caption>**Table 5.** Fixed disc model parameters.</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>stellar mass</td><td>M<sub>*</sub></td><td>2.47 M<sub>☉</sub></td></tr><tr><td>effective temperature</td><td>T<sub>eff</sub></td><td>9250 K</td></tr><tr><td>stellar luminosity</td><td>L<sub>*</sub></td><td>37.7 L<sub>☉</sub></td></tr><tr><td>inner disc radius</td><td>R<sub>in</sub></td><td>0.45 AU</td></tr><tr><td>dust material mass density</td><td>ρ<sub>gr</sub></td><td>3.36 g cm<sup>-3</sup></td></tr><tr><td>strength of incident ISM UV</td><td>χ<sup>ISM</sup></td><td>1</td></tr><tr><td>cosmic ray H<sub>2</sub> ionisation rate</td><td>ζ<sub>CR</sub></td><td>1.7 × 10<sup>-17</sup> s<sup>-1</sup></td></tr><tr><td>turbulent Doppler width</td><td>ν<sub>turb</sub></td><td>0.15 km s<sup>-1</sup></td></tr><tr><td>α viscosity parameter</td><td>α</td><td>0</td></tr><tr><td>disc inclination</td><td>i</td><td>50°</td></tr><tr><td>distance</td><td>d</td><td>118.6 pc</td></tr></tbody></table>

For the purposes of our modelling, we adopt the surface density profile derived by Hughes et al. (2008) for this object, and attempt to fit simultaneously the remaining observational data (SED, line fluxes, CO line profiles) by varying the remaining disc parameters. As an additional test, we compute a model with a power-law density profile, varying in addition the column density power index to fit the same data (SED, line fluxes, CO line profiles). For the purposes of this power-law model the disc outer edge is fixed at the scattered light radius of 540 AU. These adopted disc outer radii are consistent with our far-infrared photometry measurements, in which the disc is unresolved, corresponding to the emission at these wavelengths being dominated by the inner ~200–300 AU of the disc.

The disc is also resolved at the inner edge, and while there is some uncertainty regarding the precise structure in this region (Renard et al. 2010; Benisty et al. 2010; Eisner et al. 2009; Tannirkulam et al. 2008b), we adopt the median value of 0.45 AU for the inner disc radius. The disc is assumed to have constant flaring with radius, i.e. the assumed inner rim structure does not cause shadowing. The reference scale heights referred to in our results (see Table 6) are defined at the inner rim, i.e. $R = 0.45$ AU. The disc inclination is fairly well-constrained by imaging at various wavelengths (Benisty et al. 2010; Wassell et al. 2006; Isella et al. 2007; Tannirkulam et al. 2008a; Grady et al. 2000), and is fixed at 50° throughout the modelling in this paper. See Table 5 for a full list of the parameters which are fixed during the model-fitting process.

While we have not used the spatial CO emission maps as a constraint for our modelling efforts, we do use the extracted line profiles for the three CO lines (Isella et al. 2007) as a constraint for the radial origin of the emission. The interferometry allows us to isolate the emission from the disc, eliminating contamination from cloud material. The bulk of the modelling uses a surface density profile derived from fitting to the spatial CO 3–2 emission (Hughes et al. 2008).

### 4. Modelling

The disc modelling was carried out using the radiation thermochemical disc code ProDiMo (Woitke et al. 2009; Kamp et al. 2010). The code solves in turn the radiative transfer problem, chemical network, and gas heating and cooling balance. Finally, the level populations are calculated, followed by line transfer calculations to give the predicted line emission.

The frequency-dependent 2D radiative transfer solver calculates the dust temperature structure and internal continuous radiation field for a given disc and stellar spectrum. For the purposes of this paper it was used to determine an SED-fitting dust model
<table><caption>**Table 6.** Model parameters for four well-fitting disc models.</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>(1)<br>"low UV"</th><th>(2)<br>"high UV"</th><th>(3)*<br>Fully mixed</th><th>(4)<br>Power-law density profile</th></tr></thead><tbody><tr><td>disc mass</td><td>M<sub>disc</sub></td><td>2.2 × 10<sup>-2</sup> M<sub>☉</sub></td><td>2.0 × 10<sup>-2</sup> M<sub>☉</sub></td><td>7.1 × 10<sup>-2</sup> M<sub>☉</sub></td><td>1.1 × 10<sup>-2</sup> M<sub>☉</sub></td></tr><tr><td>column density power index</td><td>ε</td><td>n/a</td><td>n/a</td><td>n/a</td><td>0.085</td></tr><tr><td>reference scale height</td><td>H<sub>0</sub></td><td>0.019 AU</td><td>0.019 AU</td><td>0.019 AU</td><td>0.027 AU</td></tr><tr><td>flaring power</td><td>β</td><td>1.068</td><td>1.066</td><td>1.066</td><td>1.019</td></tr><tr><td>gas-to-dust mass ratio</td><td>ρ/ρ<sub>d</sub></td><td>21.8</td><td>20.5</td><td>101.1</td><td>9.1</td></tr><tr><td>minimum dust particle radius</td><td>a<sub>min</sub></td><td>5.95 × 10<sup>-3</sup> µm</td><td>1.25 × 10<sup>-3</sup> µm</td><td>9.63 × 10<sup>-3</sup> µm</td><td>1.24 × 10<sup>-2</sup> µm</td></tr><tr><td>maximum dust particle radius</td><td>a<sub>max</sub></td><td>1134 µm</td><td>1405 µm</td><td>2041 µm</td><td>1015 µm</td></tr><tr><td>dust size dist. power index</td><td>p</td><td>3.61</td><td>3.57</td><td>3.68</td><td>3.75</td></tr><tr><td>dust settling parameter</td><td>s</td><td>0.57</td><td>0.39</td><td>0 (fixed)</td><td>0.56</td></tr><tr><td>minimum dust settling radius</td><td>a<sub>set</sub></td><td>1.09 µm</td><td>0.44 µm</td><td>n/a</td><td>0.81</td></tr><tr><td>PAH abundance relative to ISM</td><td>f<sub>PAH</sub></td><td>4.3 × 10<sup>-2</sup></td><td>3.5 × 10<sup>-2</sup></td><td>6.8 × 10<sup>-3</sup></td><td>2.2 × 10<sup>-2</sup></td></tr><tr><td>fit to observed photometry</td><td>χ<sub>PHOT</sub></td><td>1.723</td><td>2.061</td><td>1.323</td><td>1.279</td></tr><tr><td>fit to ISO-SWS spectrum</td><td>χ<sub>ISO</sub></td><td>1.326</td><td>1.378</td><td>1.780</td><td>0.827</td></tr><tr><td>fit to observed line data</td><td>χ<sub>LINE</sub></td><td>0.661</td><td>0.741</td><td>0.742</td><td>0.607</td></tr></tbody></table>

**Notes.** The first three columns refer to models with exponentially-tapered density profiles, while model 4 is a simple power law. Models 1, 3 and 4 are for a "low UV" input spectrum while model 2 is for a "high UV" flux. Finally, models 1, 2, and 4 allow dust grains to settle towards the midplane while model 3 is fully mixed. Model 3 represents our overall "preferred" model, which is later discussed in detail, and is marked by an asterisk.

for HD 163296, by varying the total dust mass, grain size parameters and the spatial distribution of dust in the disc, including dust settling as a function of grain size (see Sect. 5.1).

Once the dust temperature structure $T_{\text{dust}}(r, z)$ and internal radiation field $J_v(r, z)$ have been computed, PRODiMo solves the chemical network assuming kinetic equilibrium (such that the net concentration of each species does not change over time), and the gas energy balance. This gives the chemical composition of the disc, the gas temperature $T_{\text{gas}}(r, z)$ and the gas emission lines, and depends on the gas-to-dust ratio and PAH abundance. The chemical network includes 973 reactions between 76 gas phase and solid ice species composed of 10 elements (Woodall et al. 2007; Schöier et al. 2005). There is a detailed treatment of UV-photoreactions (see Kamp et al. 2010), as well as $H_2$ formation on grain surfaces, vibrationally excited $H_2^*$ 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).

The level populations of the various atoms, molecules and ions are calculated using an escape probability method, and used to solve the line transfer for selected spectral lines, in this case those of CO (Flower 2001; Jankowski & Szalewicz 2005; Wernli et al. 2006; Yang et al. 2010), H$_2$O (Barber et al. 2006; Dubernet & Grosjean 2002; Faure et al. 2007; Green et al. 1993; Tennyson et al. 2001), O$_1$ (Abrahamsson et al. 2007; Bell et al. 1998; Chambaud et al. 1980; Jaquet et al. 1992; Launay & Roueff 1977) and C$_1$ (Flower & Launay 1977; Launay & Roueff 1977; Wilson & Bell 2002). See Woitke et al. (2011) for a description of the line transfer calculations, as well as recent improvements to the chemistry and gas heating/cooling balance in PRODiMo.

## 4.1. Input spectrum

The stellar input spectrum used for the modelling is a composite of the available observed UV data with a model photosphere spectrum in the wavelength range for which detailed spectral data are not available. We have accounted for the observed ultraviolet variability of this object (see Sect. 3) by running two sets of disc models with different UV input spectra, namely, one representing the "low UV" state and the other a "high UV" state, as described in Sect. 3.2. These were provided by the FUSE+STIS and the IUE average spectra, respectively. At the upper wavelength boundary we switch in each case to the GAIA PHOENIX model photosphere computed with $T_{\text{eff}} = 9250$ K, $\log g_* = 4.07$ and [Fe/H] = +0.20. The observations were de-reddened using the Fitzpatrick parameterization (Fitzpatrick 1999) with $E(B-V) = +0.15$ and $R_V = 3.1$.

The UV spectrum is of central importance to the gas modelling, both in the chemical network and in the gas heating/cooling balance. This requires UV input for wavelengths of 912 Å and above (the wavelength at which atomic hydrogen is ionised). The FUSE spectrum in this case only extends down to 970 Å, and the IUE spectrum down to 1150 Å. In order to give a full spectrum for modelling purposes, the "high UV" input is extended down to 970 Å by scaling the STIS+FUSE (low UV) spectrum accordingly. In both cases a power-law fit to the data is employed for the remaining shortfall. The total photon particle fluxes are 60% greater for the high UV input spectrum than for the low UV state.

## 4.2. Fitting procedure

An evolutionary $\chi^2$-minimisation strategy was employed to simultaneously fit the observed SED, ISO-SWS spectrum and gas line observations (see Sect. 3.5). In the case of the line observations the models were fitted to the observed fluxes, and the line widths in cases where profiles were available. The upper limits were not used as a fitting constraint, but were checked against the predicted best-fit model fluxes for any contradictions.

We have utilised a Monte-Carlo evolutionary method, varying the 11 parameters listed in Table 6 until the minimum $\chi^2$ value was reached (see Woitke et al. 2011 for more information regarding the evolutionary strategy and method of $\chi^2$ calculation). We note that the best-fit parameters depend in general on the weighting scheme employed for the three groups of observations, and in practice it was found that weighting slightly in favour of the line emission and ISO-SWS observations at the expense of the photometry gave the best fit overall. We also note the possibility that the converged parameter values correspond to a local $\chi^2$ minimum, and not a single global minimum, and we
<table><caption>**Table 7.** Dust grain composition.</caption><thead><tr><th>Dust species</th><th>Mass fraction</th><th>Optical constant ref.</th></tr></thead><tbody><tr><td>Amorphous FeMgSiO<sub>4</sub></td><td>0.745</td><td>Dorschner et al. (1995)</td></tr><tr><td>Amorphous carbon</td><td>0.15</td><td>Jager et al. (1998)</td></tr><tr><td>Crystalline Mg<sub>2</sub>SiO<sub>4</sub></td><td>0.035</td><td>Servoin & Piriou (1973)</td></tr><tr><td>Water ice</td><td>0.035</td><td>Warren & Brandt (2008)</td></tr><tr><td>Iron oxide</td><td>0.02</td><td>Henning et al. (1995)</td></tr><tr><td>Iron</td><td>0.015</td><td>Posch et al. (2003)</td></tr></tbody></table>

do not claim a unique fit to the observations. Indeed, it is clear from our modelling efforts that there exist a lot of degeneracies between the various parameters. However, each of the models discussed in the following sections are the result of several hundred generations of $\chi^2$-minimisation.

### 4.3. Dust properties

The rich solid-state spectrum of HD 163296 was observed by van den Ancker et al. (2000), using the ISO-SWS and ISO-LWS spectrometers. This study found the dust in this object to consist of amorphous silicates, iron oxide, water ice and a small fraction of crystalline silicates, with the presence of large millimetre-sized grains indicated by the continuum temperature.

We assume homogeneous and spherical dust grains (Mie theory), with a power-law size-distribution defined by a minimum radius $a_{min}$, maximum radius $a_{max}$, and power-law index $p$. The grain composition is assumed to be constant throughout the disc, and we adopt the dust species mixture determined by Bouwman et al. (2000) for this object, averaging their fractional species abundances over the disc (see Table 7). All grains have the same species composition regardless of their size, and the grain composition is not treated consistently with the physics and chemistry of grains. For instance, the assumed grain water ice fraction is constant throughout the disc, independent of the water ice abundances calculated in the solution of the chemical network. For the dust material mass density, we take the average value for this mix of $3.36 \text{ g cm}^{-3}$. Optical constants for the various dust species were taken from measurements made by the studies listed in Table 7.

Bouwman et al. (2000) carried out a detailed study of this object, fitting to spectral data over a large wavelength range. We note however that the dust opacity law represents a potentially large source of uncertainty when deriving the disc parameters. For instance, the two olivine species FeMgSiO$_4$ and Mg$_2$SiO$_4$, whilst having broadly similar spectra in the wavelength range observable by ISO-SWS, have a factor >30 difference in absorption opacity at ~1 micron. This will inevitably affect the derived disc parameters, although the difference in opacities at millimetre wavelengths is negligible. The models referred to in this study have extinction opacities in the range $(10.4-12.8) \text{ cm}^2 \text{ g}^{-1}$ (dust) at 1 mm.

### 4.4. Gas/dust ratio

The models considered have constant dust properties with radius. However, in general the models do allow for differential dust grain scale heights as a function of grain size, to represent the major effects of dust settling. This leads to a local gas/dust ratio which varies with height in the disc, yet at any given radius the ratio of the vertical gas and dust total column densities is constant, and equal to the global gas/dust ratio. The degree of dust settling is determined by a simple parameterisation (see

Sect. 5.1), and models in which no dust settling is present are referred to as “fully mixed”.

## 5. Results

Table 6 gives the derived parameter values for four runs of the evolutionary scheme outlined in the previous section. The first three columns refer to models with a fixed exponentially-tapered density profile. The fourth column refers to a model with power-law density profile, $\Sigma \propto R^{-\epsilon}$, in which the outer radius is fixed at 540 AU but the power index, $\epsilon$, is allowed to vary.

The surface density profile in the models with an exponential outer edge (Col s. 1–3) takes the form

$$ \Sigma(R) \propto R^{-\gamma} \exp \left[ - \left( \frac{R}{R_0} \right)^{2-\gamma} \right], \quad (1) $$

where $\gamma$ is analagous to the power law index, and $R_0$ is the scale length over which the disc surface density tapers exponentially. We adopt the previously derived values of $\gamma = 0.9$ and $R_0 = 125$ AU (Hughes et al. 2008). For the purposes of our modelling, the disc chemistry etc. was computed out to 850 AU, at which point the column density is negligible.

The first, “low UV” column uses the FUSE+STIS UV spectrum as input for the modelling, and the second “high UV” column uses the average IUE spectrum, as described in Sect. 4.1. Both runs produce settled discs, with variable scale heights for dust grains of different sizes. The third column in Table 6 gives the results for a run with well-mixed dust grains, i.e. no settling is introduced. This run uses the FUSE+STIS (low UV) spectrum, as does the run with power-law surface density referred to in Col. 4.

None of the models described in Table 6 represent a perfect fit to all the available data for HD 163296, but in all cases the models are able to fit almost all the observations. Considering the large number of data used for the modelling, and the range of parameters covered by the well-fitting models, is is clear that there is some degree of parameter degeneracy present.

### 5.1. Continuum emission and spectral energy distribution

The SED for the best-fit fully mixed disc model (Col. 3 in Table 6) is shown in Fig. 4. With the constraints outlined in previous sections, our aim was to fit the observed SED with a simple dust model, i.e. a continuous disc with constant flaring, and dust species composition constant throughout the disc. All of the models described in Table 6 require the presence of mm-sized grains in order to fit the observed SED. It was not possible to fit both the 10 micron silicate feature and the millimetre tail with a well-mixed disc. This can be seen in Table 6, where the fully mixed model gives a worse fit to the ISO-SWS spectrum than for the models in which dust settling is present. Figure 4 shows that the best-fit fully mixed model gives a smaller than observed silicate emission feature. Sitko et al. (2008) observe a spectral variability of ~10% in the wavelength range covered by the ISO-SWS data, and while we are unable to fit the observations to within this range with even a settled model (simultaneously with the overall SED and line data), the fit is considerably better than is possible with a well-mixed disc. This result is similar to that found for IM Lupi by Pinte et al. (2008), in which a settled disc was needed to fit the 10 µm silicate feature and SED millimetre tail simultaneously.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** 1.3 mm continuum emission maps for HD 163296. *From left:* observed map; emission map for preferred model (Col. 3 in Table 6); residuals for preferred model; emission map for power-law model (Col. 4 in Table 6); residuals for power-law model. Residuals computed by subtracting the model intensity from the observed intensity. Contours spaced at 12 mJy intervals, corresponding to [3,6,9,12,15,18,21,24,27,30,33,36,39,42,45]×σ.</figcaption></figure>

On the other hand, the constraint of fitting to the ISO-SWS spectrum leads to models which produce too much emission in the near-infrared (J, H and K bands), giving an overall worse fit to the observed photometry. This effect is less pronounced in the best-fit fully mixed model (Col. 3 cf. Cols. 1 and 2 in Table 6). The flux overprediction in the settled models is likely mainly due to the smaller average grain size in the strongly illuminated disc surface layer. The condition of radiative equilibrium between grains means that this gives a higher average grain temperature at the disc surface, and the re-emitted light peaks at shorter wavelength than for the fully mixed models (in the near-IR as opposed to the mid-IR). This flux overprediction is still present to a smaller extent in the fully mixed model (see Fig. 4). This could be a consequence of our fixing the disc inner radius for modelling purposes, and a better fit in the near-IR would be possible if we had allowed the inner radius to vary. It is also probable that the structure of the inner rim is more complex than allowed by our parameterised model, with evidence for puffing-up of material, which would certainly affect the emission at these wavelengths. There is also significant variability observed for this object in this wavelength region (de Winter et al. 2001; Sitko et al. 2008).

Dust settling enhances the silicate emission feature for an optically thick disc since it removes the large grains from the surface layers, and places them at smaller heights in the disc where they cannot be observed at mid-infrared wavelengths. The flat blackbody opacity of these larger grains is overwhelmed by the characteristic spectrally-varying opacity of smaller micron-size grains which remain at lower optical depth, dominating the observed emission.

In cases such as this where a simple parameterised disc structure is assumed, PRODiMo implements a simple recipe to account for the major effects of vertical dust settling. We assume that the dust grains are distributed vertically with a scale height which decreases for large dust particles:

$$H'(a, r) = H(r)\max\{1, a/a_s\}^{-s/2} \quad (2)$$

where $H(r)$ is the gas scale height, and $s$ and $a_s$ are two free parameters (see Woitke et al. 2010 for details).

Introducing two additional parameters to the modelling inevitably improves the fit to the observed data, and indeed the combined fit to the photometry, ISO-SWS spectrum and line emission does improve overall as dust settling is introduced. The fact remains, however, that it is possible to obtain a good fit to almost all the available data with a well-mixed disc, with a gas/dust ratio almost equal to the canonical value of 100, and so we adopt this model, described in Col. 3 of Table 6, as our

“preferred” model. However, the inability of a well-mixed disc such as this to fit the observed emission in the mid-infrared means that there remains compelling evidence for a disc exhibiting dust-settling, able to fit the line emission with a depleted gas/dust ratio (see Col. 1 in Table 6). The effect of dust settling on the line emission is further explored in Sect. 5.2.2.

The model with power-law density profile (Col. 4 in Table 6) gives the best overall fit to the SED, ISO-SWS spectrum and line emission (see Fig. 4), with quite different disc parameters to the three models with exponential outer edges (Cols. 1–3). The disc requires a very flat density profile, $\epsilon \sim 0.085$, in order to fit these observations (not including spatial data). The predicted millimetre continuum emission from this model is plotted with the observed maps in Fig. 5. It is clear that the flat density profile leads to an emission deficit in the inner disc in comparison to observations, and so this cannot be considered to be a good model for the disc of HD 163296. A better power-law model could be found by using the maps as a further constraint for the modelling, but that is beyond the scope of this paper. The fully mixed “preferred” model gives a better match to the observed millimetre emission, as should be expected since its density profile is derived from fitting to this data (Hughes et al. 2008). The factors driving the power-law model to such a flat density profile will be discussed in Sect. 5.2.

The power-law density profile model is, however, able to fit the non-spatially-resolved data with a smaller flaring index (~1.02) than the exponentially-tapered discs (~1.07), more in keeping with the observational evidence for a disc that is not strongly flaring (Meeus et al. 2001; Doucet et al. 2006; Wisniewski et al. 2008). Meeus et al. (2001) classify HD 163296 as a group II object, in which the inner disc shadows the outer disc, as opposed to group I objects which have a flared geometry. Meijer et al. (2008) find HD 163296 to lie on the transition between flared and non-flared geometry, i.e. a flaring index close to 1. It is clear that our “preferred” model is unable to account for every observed property of the disc of HD 163296, and it is probable that a more complex model – with non-constant flaring and variable dust properties with radius – would be required to simultaneously fit the full set of observational data available for this object.

The models described in Table 6 have disc dust masses in the range $(7-12) \times 10^{-4} M_{\odot}$, which is within the range of masses $(5-17) \times 10^{-4} M_{\odot}$ found in the literature for HD 163296. This small spread in dust mass from our fitting efforts is unsurprising given the fixed disc size and constant grain composition.
<table><caption>**Table 8.** Integrated line fluxes for the transitions observed by PACS, and additional CO and H₂ lines, as predicted by the models listed in Table 6.</caption><thead><tr><th>Species</th><th>λ [µm]</th><th>ν [GHz]</th><th>"low UV" model</th><th>"high UV" model</th><th>Fully mixed model</th><th>Power-law</th><th>Observed flux</th></tr></thead><tbody><tr><td>Oi</td><td>63.18</td><td>4745.05</td><td>200.6</td><td>222.0</td><td>191.4</td><td>170.3</td><td>193.1 (58.2)</td></tr><tr><td>Ot</td><td>145.52</td><td>2060.15</td><td>6.29</td><td>7.33</td><td>5.39</td><td>4.59</td><td>&lt;8.5</td></tr><tr><td>Cii</td><td>157.74</td><td>1900.55</td><td>8.32</td><td>8.40</td><td>11.2</td><td>9.98</td><td>—</td></tr><tr><td>p-H₂O 3₂₂ → 2₁₁</td><td>89.99</td><td>3331.40</td><td>1.22</td><td>1.68</td><td>1.20</td><td>0.353</td><td>&lt;9.4</td></tr><tr><td>o-H₂O 2₁₂ → 1₀₁</td><td>179.52</td><td>1669.97</td><td>4.40</td><td>4.19</td><td>2.43</td><td>3.94</td><td>&lt;14.5</td></tr><tr><td>o-H₂O 2₂₁ → 2₁₂</td><td>180.49</td><td>1661.64</td><td>0.847</td><td>0.993</td><td>0.618</td><td>0.243</td><td>&lt;16.2</td></tr><tr><td>o-H₂O 4₃₂ → 3₁₂</td><td>78.74</td><td>3810.01</td><td>2.33</td><td>2.94</td><td>2.20</td><td>0.840</td><td>&lt;15.0</td></tr><tr><td>o-H₂ S(1)</td><td>17.03</td><td>17603.78</td><td>1.15</td><td>1.59</td><td>1.15</td><td>0.572</td><td>&lt;28</td></tr><tr><td>OH</td><td>79.11</td><td>3792.19</td><td>4.59</td><td>6.31</td><td>5.40</td><td>2.71</td><td>&lt;17.0</td></tr><tr><td>OH</td><td>79.18</td><td>3788.84</td><td>4.78</td><td>6.54</td><td>5.59</td><td>2.82</td><td>&lt;17.0</td></tr><tr><td>CO J = 36–35</td><td>72.85</td><td>4115.20</td><td>0.335</td><td>0.489</td><td>0.423</td><td>0.132</td><td>&lt;11.6</td></tr><tr><td>CO J = 33–32</td><td>79.36</td><td>3777.63</td><td>0.497</td><td>0.704</td><td>0.578</td><td>0.197</td><td>&lt;22.8</td></tr><tr><td>CO J = 29–28</td><td>90.16</td><td>3325.12</td><td>0.763</td><td>1.03</td><td>0.732</td><td>0.311</td><td>&lt;11.1</td></tr><tr><td>CO J = 18–17</td><td>144.78</td><td>2070.68</td><td>1.91</td><td>2.31</td><td>1.49</td><td>0.851</td><td>&lt;13.1</td></tr><tr><td>CO J = 3–2</td><td>866.96</td><td>345.80</td><td>1.26</td><td>1.27</td><td>1.22</td><td>1.49</td><td>1.65 (0.39)</td></tr><tr><td>CO J = 2–1</td><td>1300.40</td><td>230.54</td><td>0.401</td><td>0.401</td><td>0.397</td><td>0.464</td><td>0.379 (0.118)</td></tr><tr><td><sup>13</sup>CO J = 1–0</td><td>2720.41</td><td>110.20</td><td>0.0122</td><td>0.0115</td><td>0.0120</td><td>0.0127</td><td>0.0124 (0.007)</td></tr></tbody></table>

**Notes.** Observed fluxes for comparison. Fluxes in $[10^{-18} \text{ W/m}^2]$. Errors are a quadratic sum of the calibration error (20% for the interferometric CO observations and 30% for the PACS observations) and the rms continuum noise.

## 5.2. Gas properties

The predicted line fluxes for the various models are given in Table 8. We have ignored the [Cii] 158 µm result during the model-fitting process, due to systematic uncertainties arising from strong emission at the offset positions.

The observed line profiles for the ¹²CO J = 3–2, ¹²CO J = 2–1 and ¹³CO J = 1–0 transitions (Isella et al. 2007) are plotted with the preferred model profiles in Fig. 6. The observed profiles are double-peaked, consistent with a disc in Keplerian rotation. The CO 3–2 line is observed to be brighter than predicted by the model. This behaviour is repeated across the three exponentially-tapered models, but is less pronounced in the power-law model, whose flat surface density profile gives a larger CO 3–2 flux. Discs with flatter density profiles allow stellar radiation to penetrate deeper into the inner disc, resulting in a slight increase in temperature throughout the disc, including its outer regions. This leads to brighter CO 3–2 lines since this transition is optically thick in all of our models, and is largely dependent on the disc outer radius (which is fixed) and the outer disc temperature at intermediate height (see Fig. 8). The CO 2–1 line also follows this behaviour, with the line flux ratio in our models staying roughly constant at a value of CO 3–2/1–0 ~ 3, as expected for lines formed under optically thick LTE conditions (Kamp et al. 2010). This is slightly lower than the observed ratio of CO 3–2/2–1 = 4.3. We note that recent observations of this object give a peak intensity of 25 Jy for the 3–2 line (Hughes et al. 2011)<fnref n="5" />, compared with 43 Jy in the Isella et al. (2007) data and the 30 Jy predicted by our preferred model. Our predicted flux is within the calibration error margin for both observations, but fitting to the lower value would in general tend to lead to power-law models with steeper density profiles, in contrast to the flat power-law profile we obtain from fitting to the brighter observation. Another factor driving the power-law models to flat density profiles is the overprediction in near-IR continuum emission by our models. The inner radius is fixed in accordance with high

resolution imaging, and so the power-law models tend to reduce their near-IR emission by removing material from the inner disc, giving a better fit to the photometry.

PAHs are one of the main sources of gas heating in the disc, as UV photons cause them to emit excited electrons via the photoelectric effect, which thermalise in the gas. For the purposes of our modelling we consider a typical size of PAH molecule (circumcoronene, $N_C$ = 54 carbon atoms and $N_H$ = 18 hydrogen atoms) and include PAH⁻, PAH, PAH⁺, PAH²⁺ and PAH³⁺ as additional species in the chemical network (see Woitke et al. 2010 for further details). The fractional PAH abundance, $f_{\text{PAH}}$, is defined relative to the standard ISM particle abundance with respect to hydrogen nuclei, $X_{\text{PAH}}^{\text{ISM}} = 3 \times 10^{-7}$ (Tielens 2008), so that the PAH abundance, $\epsilon(\text{PAH}) = f_{\text{PAH}} X_{\text{PAH}}^{\text{ISM}} / N_C$. Small $f_{\text{PAH}}$ values (0.007–0.04 relative to ISM abundances) are required in our models in order to fit the observed line data. This is consistent with the findings of Geers et al. (2006), who require PAH abundances $f_{\text{PAH}} < 0.1$ in order to fit the observed PAH emission from Herbig discs. PAH emission was not observed in HD 163296 by Acke et al. (2010), and we have used the Monte Carlo radiative transfer code McFost (Pinte et al. 2006) to check that our results are consistent with this non-detection. A disc model with parameters and PAH abundance identical to our preferred model gives an infrared spectrum in which no PAH emission is seen. This is due in part to the low flaring in the disc, $\beta = 1.066$, and the large fraction of small dust grains in the disc, with significant opacity in the optical-UV range over which PAH molecules absorb. This results in the PAHs being “hidden” from direct stellar illumination, absorbing a fraction ~$10^{-4}$ of the energy absorbed by the disc. In our preferred model 0.3% of the total carbon mass is tied up in PAH molecules.

The total hydrogen number density and gas temperature within the disc are plotted in Fig. 7, and the spatial origin of the various emission lines is visualised in Fig. 8. These plots all refer to the “preferred” disc model, i.e. one in which no dust settling is present. The CO J = 3–2 and [Oi] 63 µm lines are optically thick throughout the disc, and cannot be used alone to trace the gas mass. Even the ¹³CO line is optically thick throughout much of the disc, only becoming optically thin outside of ~400 AU. Also, there is evidence that this line (and most others) can be

<footnote n="5">⁵ This CO 3–2 data also seems to indicate the presence of turbulence of order 0.3 km s⁻¹ in the disc of HD 163296 (Hughes et al. 2011). We have not yet explored the effect of this possible stronger turbulence on our results.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Plotted in black are the line profiles observed by Isella et al. (2007) for the ¹²CO J = 3–2 (upper panel), ¹²CO J = 2–1 (middle panel) and ¹³CO J = 1–0 (lower panel) transitions, with the corresponding profiles from the preferred disc model in red. This refers to a single simultaneous fit to the observed continuum and line data (Col. 3 in Table 6). Observed profiles are obtained by integrating over the whole disc.</figcaption></figure>

affected by the degree of dust settling in the disc (see Sect. 5.2.2). The [CII] 157.74 µm line is optically thin, but traces only the warm ionised gas in the disc surface. Clearly care must be taken when attempting to use individual emission lines as a tracer of

gas mass. Note the contrast in disc radii at which line and continuum become optically thin in the case of CO 3–2, reflecting the conflicting derived radii for the gas and dust discs in HD 163296. The line optical depths in Fig. 8 were computed using an escape probability formalism, and do not take into account the effects of Keplerian shear in the disc.

As well as optical depth effects, disc gas mass estimates from CO emission alone can be affected by the extent of CO freeze-out on grain surfaces. PRODiMo calculates the grain adsorption and desorption as part of the solution of the chemical network, and the spatial abundances of gas-phase and ice CO in our best-fit model are plotted in Fig. 9. The total mass of gas-phase CO in the disc is $1.51 \times 10^{-5} M_{\odot}$, and the CO ice mass is $4.32 \times 10^{-5} M_{\odot}$. We would expect a smaller fraction of CO to be frozen out on grain surfaces in discs with flatter density profiles, due to weaker dust shielding in the inner disc leading to warmer temperatures throughout the disc, and less CO freeze-out.

The derived chemical abundances and gas properties in our preferred model have been checked by re-computing the chemistry using a time-dependant solver, and comparing the results to those obtained with the assumption of kinetic equilibrium. We assume molecular cloud initial abundances, and after running the solver for 4 Myr there is no major departure from the equilibrium chemistry. The biggest change in line flux is a reduction of 20% in the 180.4 µm water line, while the [OII] 63 µm line decreases by 4% and the CO lines by <1%. The assumption of a constant ¹³CO/¹²CO ratio is also valid since dust shielding dominates over CO self-shielding in the models, with negligible change to the results when CO self-shielding is swtiched off. We would therefore expect there to be no fractionation effect present in this disc.

All of the disc models considered are passive discs, i.e. the viscous “$\alpha$” heating parameter is set to zero. This is not strictly consistent with the presence of accreting gas in this object, since this implies some form of viscosity to transport angular momentum through the disc. However, viscous heating is likely to be unimportant in the case of Herbig Ae discs such as this, which are thought to be dominated by radiative heating by the central star (D’Alessio et al. 1998).

### 5.2.1. Effect of UV variability

We have investigated the effects of UV variability on the gas in the disc in two ways. Firstly, by computing a test model with parameters identical to our “preferred” model, but with the high UV state spectrum used as input. This allows us to isolate the effect of increasing the UV intensity on the line emission, since all other model parameters remain the same. The results of this test are summarised in Table 9. The second approach consists of a separate evolutionary run of models using the high UV state as input, and has been covered in Sect. 5. This method allows us to assess how the disc model parameters change in order to fit the same line data with an increased UV intensity (see Table 6). By examining firstly the physical effect of the UV on the gas in the disc, and secondly the effect of the UV on the model parameter fit, we can estimate the degree of uncertainty introduced by UV variability in this object.

An increased UV strength affects the gas chemistry by promoting photochemical reactions, and the photodissociation of H₂ and CO. It also leads to increased desorption of ice species from grain surfaces in regions where the UV is able to penetrate. It also strongly affects the gas heating, via the photoelectric effect in dust grains and PAHs.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** The total hydrogen number density (**left panel**) and gas temperature (**right panel**) are plotted for a vertical cross-section through the preferred disc model (Col. 3 in Table 6). Dashed lines indicate contours of gas temperature and visual extinction $A_V$, as marked.</figcaption></figure>

The “high UV” input spectrum represents a 60% increase in UV intensity over the low UV input. Table 9 summarises the effects of this increase on the gas in the preferred model. There is a slight overall increase in the disc gas temperature, and an increase in the atomic hydrogen and $C^+$ masses, due to photoelectric effects and PAH heating, and increased photodissociation of gas molecules. The total gas-phase CO mass actually increases slightly, since the increased rate of photodissociation is balanced by an increase in the rate of ice desorption. The $J = 36–35$ line in particular sees the biggest increase in flux, ~50%. By considering only the gas in the warm layers of the inner disc, where the high $J$ CO lines form, it can be seen that this increase in line flux is caused by a substantial increase in the gas temperature in this region. The [OI] 63 µm flux increases by ~30%, roughly equal to the calibration uncertainty margin, and the low $J$ CO lines increase by ~5%. All non-detected lines stay within the observed upper limits. In summary, it would not be possible to distinguish between the high and low UV states from these disc observations, since all fluxes are within the observational error margins. The general influence of the UV intensity on the disc gas is discussed by Woitke et al. (2010), in the context of a large grid of models.

The derived parameters for the evolutionary run using the “high UV” spectrum as input are listed in Table 6. They are largely similar to those derived from the low UV run (Col. 2 cf. 1), as would be expected from the small fractional change in UV intensity. The slight decrease in gas mass and flaring will tend to counteract any increase in line emission. The PAH fractional abundance decreases slightly as might be expected. The main change lies in the grain size distribution, with a reduction in the minimum grain size. This makes it harder for the UV to penetrate the disc, balancing the increase in UV intensity. This also produces even more continuum emission in the near-IR, worsening the SED fit still further in comparison to the “low UV” model.

### 5.2.2. Effects of dust settling

The effect of dust settling on various line fluxes can be seen in Fig. 10. This shows the increase in line flux across the various

transitions as smaller and smaller grains are allowed to settle towards the midplane. The models are otherwise identical to the preferred (fully mixed) model, but with a settling parameter of 0.5. The models run from strongly settled discs through to entirely well-mixed, i.e. identical to the preferred model. For the purposes of this exercise we have computed some extra line fluxes in addition to those observed in HD 163296, namely the $^{12}$CO 6–5 and 1–0 transitions, the $^{13}$CO 3–2 transition, and the $^{1}$H$_2$O 538.3 µm transition.

Kamp et al. (2011) computed a large grid of disc models using PRODiMo, in tandem with McFost (Pinte et al. 2006). This study noted a trend of increasing disc dust temperature with settling, due to a combination of reduced emissivity at long wavelengths in the disc surface, and increased illumination of the dust as the stellar radiation is able to penetrate further into the disc. Our models of HD 163296 follow this same behaviour, albeit with an accompanying increase in gas temperature. The effect of dust settling on the gas and dust temperature contours in our models is demonstrated in Fig. 11, where the chemical abundances of various chemical species are plotted. The chemical structure is seen to “follow” the settled dust grains towards the midplane, since the UV penetrates further into the disc, causing the characteristic layering of the various gas species to move closer to the midplane. This is despite the gas scale heights being identical in both models. In the settled models the warm gas in the disc resides at lower heights where the gas density is higher, leading to a general brightening of the emission lines. This seems to be a firm result for this object, with serious implications for the derivation of a precise disc gas mass from the emission lines.

The brightening of the line fluxes seems to be a general result for all the transitions considered. The exact behaviour of the various line fluxes with variable dust settling, and the extent to which they are affected, depends on the height and radial position in the disc from which they originate. This is well-illustrated by the CO lines in Fig. 10. As the largest grains begin to settle, the $^{13}$CO 1–0 line is the first to show an increase in flux. This is due to an increase in thermal desorption of CO ice from grain surfaces in the outer disc midplane, which this line traces (see Fig. 8). This behaviour is echoed in the other $^{13}$CO line,
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Spatial origin of the various gas emission lines in the preferred model (Col. 3 in Table 6). From top-left clockwise: CO J = 3–2, CO J = 18–17, [Cu] 157.74 µm, o-H₂ S(1) (with red Tg contours), [Oi] 63.18 µm, ¹³CO J = 1–0. In each panel, the *upper plot* shows the line optical depth as a function of radius (blue line) and the continuum optical depth at the corresponding wavelength (black line). The *middle plot* shows the cumulative contribution to the total line flux with increasing radius. The *lower plot* shows the gas species number density, and the black line marks the cells that contribute the most to the line flux in their vertical column. The two vertical dashed lines indicate 15% and 85% of the radially cumulative face-on line flux respectively, i.e. 70% of the line flux originates from within the two dashed lines.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** CO abundances in the preferred model (Col. 3 in Table 6). Left panel: gas-phase CO abundances. White and blue dashed lines plot the gas and dust temperature contours respectively. Right panel: CO ice abundance.</figcaption></figure>

<table><caption>**Table 9.** Effects of UV variability. Comparison between the preferred model, and an identical model illuminated by an input spectrum typical of a “high UV” state for this object.</caption><thead><tr><th></th><th>Best-fit model</th><th>Test model (high UV)</th></tr></thead><tbody><tr><td>L<sub>UV</sub>/L<sub>*</sub></td><td>0.097</td><td>0.155</td></tr><tr><td>M<sub>H</sub></td><td>2.67 × 10<sup>-5</sup> M<sub>☉</sub></td><td>2.82 × 10<sup>-5</sup> M<sub>☉</sub></td></tr><tr><td>M<sub>C+</sub></td><td>8.67 × 10<sup>-8</sup> M<sub>☉</sub></td><td>9.42 × 10<sup>-8</sup> M<sub>☉</sub></td></tr><tr><td>M<sub>CO</sub></td><td>1.54 × 10<sup>-5</sup> M<sub>☉</sub></td><td>1.66 × 10<sup>-5</sup> M<sub>☉</sub></td></tr><tr><td>⟨T<sub>CO</sub>⟩</td><td>55.3 K</td><td>55.0 K</td></tr><tr><td>⟨T<sub>CO</sub>⟩ (z/R > 0.2, R < 20 AU)</td><td>262.8 K</td><td>332.6 K</td></tr><tr><td>F<sub>[Oi] 63 µm</sub></td><td>1.91 × 10<sup>-16</sup> W m<sup>-2</sup></td><td>2.49 × 10<sup>-16</sup> W m<sup>-2</sup></td></tr><tr><td>F<sub>[Cii] 58 µm</sub></td><td>1.12 × 10<sup>-17</sup> W m<sup>-2</sup></td><td>1.41 × 10<sup>-17</sup> W m<sup>-2</sup></td></tr><tr><td>F<sub>CO 3-2</sub></td><td>1.22 × 10<sup>-18</sup> W m<sup>-2</sup></td><td>1.31 × 10<sup>-18</sup> W m<sup>-2</sup></td></tr><tr><td>F<sub>CO 36-35</sub></td><td>4.23 × 10<sup>-19</sup> W m<sup>-2</sup></td><td>6.22 × 10<sup>-19</sup> W m<sup>-2</sup></td></tr></tbody></table>

**Notes.** Shown are total species masses, average overall CO temperature, average CO temperature for the warm gas regions in the inner disc, and line fluxes. $L_{UV}$ refers to the energy emitted in the wavelength range 912–2500 Å.

the $J = 3-2$ transition. In contrast, the $^{12}$CO 3–2 line, which is also formed in the outer disc, shows a smaller increase in flux. This is because this line is optically thick, and formed higher in the disc where it remains unaffected by the thermal desorption of ice in the midplane. The $^{12}$CO 2–1 and 1–0 lines follow this same behaviour (not plotted). As smaller grains begin to settle, the $^{13}$CO lines increase less rapidly than for the lines formed in the warm gas higher in the disc. This is most pronounced in the high $J$ CO lines, which form in the inner disc where the temperature and density gradients are more extreme, and the lines more sensitive to the downwards shift in chemical structure. The $^{12}$CO 6–5 line forms at intermediate distance in the disc, and this is reflected in the level of flux enhancement (Fig. 10; top panel).

The behaviour of the water lines (middle panel in Fig. 10) is less straightforward. The $\mathrm{O-H}_2\mathrm{O}$ 179.5 µm line shows a stronger increase in flux than both the 180.5 µm line, which is at higher excitation, and the 538.3 µm line, which is at lower excitation. The $\mathrm{p-H}_2\mathrm{O}$ 89.99 µm and $\mathrm{o-H}_2\mathrm{O}$ 78.7 µm lines are at higher excitation, forming in the inner disc, and these also show a large increase in flux as the dust grains settle. It is likely that the relative behaviour of the various water lines is a complex function of the changing temperature structure and $\mathrm{H}_2\mathrm{O}$ chemical

abundance structure in the disc, affecting the excitation of the various levels.

The largest increase in line flux occurs in the $\mathrm{H}_2\mathrm{S(1)}$ transition (Fig. 10; lower panel). The emission region for this line is seen to follow the 300 K gas temperature contour (see Fig. 8), which intersects a greater molecular hydrogen mass in the settled models. This leads to a dramatic enhancement in flux in the strongly settled models, a factor of ~30 increase on the fully mixed model. The $[\mathrm{Oi}]$ 63 µm line also starts to increase as the smallest grains are allowed to settle, coinciding with an increase in atomic gas, but the enhancement is less extreme than since this line is formed further out in the disc (see Fig. 8). The $\mathrm{OH}$ 79.11 µm line is formed at intermediate distance between $\mathrm{H}_2\mathrm{S(1)}$ and $[\mathrm{Oi}]$ 63 µm, and its behaviour with increasing dust settling reflects this. The $[\mathrm{Cu}]$ 158 µm line is less strongly affected by settling. This line is formed in a thin layer at the disc surface, with emission dominated by the outer disc, and so it is less sensitive to changes in the internal disc temperature structure.

The settling of dust grains might be expected to drive line formation conditions closer to LTE, with the various species
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** The fractional enhancement in line flux is plotted against the minimum grain size affected by dust settling, relative to the preferred (fully mixed) model ($a_s = a_{max}$). All models have a settling parameter of 0.5. Top panel: CO lines; middle panel: water lines; lower panel: [OⅠ]63 µm, [CuⅠ]158 µm, OH 79.11 µm and H₂ S(1) lines.</figcaption></figure>

residing lower in the disc, in regions of higher density, with more frequent collisions between particles. However, there is no evidence for this in our models. The majority of lines show only minor departures from LTE, in both settled and well-mixed models. The largest departure from the line fluxes assuming LTE level populations occurs in the water lines, but the LTE line fluxes increase at roughly the same rate as those computed using escape probability, so the ratio $F_{\text{NLTE}}/F_{\text{LTE}}$ remains roughly constant. The critical density for the water lines, $n_{\text{crit}} \sim 10^8-10^{10} \text{ cm}^{-3}$. These densities are reached only in the inner disc in our models, with each of the lines arising in regions with $n < n_{\text{crit}}$, even in the settled models. The water line which is closest to LTE in our models is the 538.3 µm line, which has the lowest level of excitation and the smallest critical density. In general, the critical densities decrease rapidly at large optical depths, driving most of the transitions towards LTE.

These results are in contrast to the findings of Jonkheid et al. (2007), who find a trend of decreasing gas temperature and line emission as the dust grains are allowed to settle towards the mid-plane. However, a direct comparison is not appropriate, since their models assumed a simultaneous reduction in the dust/gas ratio and PAH abundances with dust settling. The reduction in gas temperature can largely be attributed then to the drop in photoelectric heating from dust grains and PAH molecules. Jonkheid et al. (2007) also note that the dust temperatures in their models are high enough to prevent CO freeze-out, whereas considerable freeze-out is present here. Meijerink et al. (2009) also predict an increase in line emission with dust settling, from models in which the effect of settling is simulated by simply increasing the global gas/dust ratio, as opposed to considering the vertical distribution of grains in the disc.

Our models indicate that the ¹³CO 1–0 line flux can change by a factor ~2 in strongly settled models. This is in a high mass disc for which the line is optically thick throughout most of the disc, and so this settling flux enhancement might be expected to be even stronger in a lower mass disc. In any case, variable grain freeze-out would seem to limit the ability of this line to trace the total disc gas mass. In general, the effect of dust settling on the vertical thermal structure of the gas in discs is seen to introduce new degeneracies when attempting to fit the disc parameters to the observed line emission. While our preferred model fits the observed data with a well-mixed disc with the canonical gas/dust ratio, we cannot rule out the possibility of a gas-depleted disc in which dust settling gives an enhancement in the various line fluxes. Indeed, such a disc with gas/dust ~20 allows a better fit to the line data and the observed 10 micron silicate emission (Col. 1 in Table 6).

### 5.3. Effect of X-rays

It is currently unclear how important a role X-rays play in determining the gas chemistry and temperature structure in discs. Aresu et al. (2011) studied the effects of stellar X-ray emission on models of T Tauri discs computed with PRODiMo. This study found that while Coulomb heating by X-rays introduced an extended hot gas surface layer to the discs, the [OⅠ] and [CuⅠ] fine structure lines were only affected for X-ray luminosities $L_X > 10^{30} \text{ erg s}^{-1}$. One would expect this to hold for the warmer gas in Herbig discs, and Günther & Schmitt (2009) derive an X-ray luminosity of $L_X = 10^{29.6} \text{ erg s}^{-1}$ for HD 163296, lower than but close to the threshold value noted by Aresu et al. (2011) for T Tauri discs. HD 163296 represents an interesting test case for the influence of X-rays in Herbig discs.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Effect of dust settling on the disc chemical abundances. Left column indicates abundances in a fully mixed model, and right column represents a strongly-settled model with otherwise identical parameters. Top row: CO abundance with dust (blue dashed lines) and gas (white dashed line) temperature contours. Middle row: H₂ abundance with 300 K gas temperature contour (black dashed line) and visual extinction contour (white dashed line). Bottom row: H₂O abundance. White dashed lines indicate T<sub>gas</sub> contours enclosing “hot water” region, red dashed lines indicate contours of UV field strength per hydrogen nucleus, log((χ/n<sub>H</sub>)) as defined in Draine & Bertoldi (1996), enclosing the cool water belt.</figcaption></figure>
We have computed a disc model with parameters identical to our preferred model, but with an additional X-ray component in the input spectrum. The X-ray luminosity was set equal to the value of $L_X = 10^{29.6}$ erg s$^{-1}$ observed by Günther & Schmitt (2009) for HD 163296. The effect of this on the gas temperature is illustrated by Fig. 12 (cf. Fig. 7). The X-ray heating processes produce an extended hot surface gas layer, as noted by Aresu et al. (2011), with gas temperatures >5000 K extending out to the outer disc.

There is little discernible effect on the disc chemistry, with none of the chemical species masses or predicted line fluxes changing by more than a few percent. In all cases the effect of the additional X-rays is less than produced by switching to the “high UV” input spectrum. We would therefore expect the UV to dominate the gas chemistry in a Herbig disc of the sort considered here.

## 6. Summary and conclusions

This paper presents new observations of the far-IR lines of the Herbig Ae star HD 163296, obtained using Herschel/PACS as part of the GASPS open time key program. These consist of a detection of the [OI] 63.18 µm line, as well as upper limits for eight additional lines. We have supplemented these observations with additional line and continuum data, as well as high resolution observations of the variable optical spectra obtained as part of the EXPORT collaboration.

We have computed radiation thermo-chemical disc models using the disc code PRODiMo, and employed an evolutionary $\chi^2$-minimisation strategy to find the best simultaneous model fit to the continuum and line data. The stellar parameters and UV input spectrum for the modelling were determined through detailed analysis in the UV, optical and near-infrared.

We obtain reasonable fits to the observed photometry, ISO-SWS spectrum, far-IR line fluxes and millimetre CO profiles for a variety of discs, and note that parameter degeneracies preclude the precise derivation of the disc properties. In particular, the effects of dust settling on the vertical thermal structure of the gas in discs strongly influence the line emission, placing further limits on the derivation of a disc gas mass from the line data. We are able to fit the photometry and line data with a well-mixed disc with gas/dust ~100, but an equally good fit is possible with a disc with gas/dust ~22, in which dust settling is present. The main advantage of the settled model is that the fit to the observed ISO-SWS spectrum is greatly improved, and it seems probable that some degree of dust settling and gas-depletion is present in the disc. These discs have a radial surface density profile as derived by Hughes et al. (2008), with an exponentially-tapered outer edge.

A model fit to the same data with power-law radial density profile gives a gas/dust ratio of ~9. This power-law model gives the best fit to the observed SED data, although it does very badly when it comes to fitting the radial intensity distribution of the resolved millimetre continuum. This is due to the flat density power index required to fit the SED and the observed CO 3–2 emission, and presents a challenge for future modelling efforts.

The emitted line fluxes are in general sensitive to the degree of dust settling in the disc. This is a firm result in our models, and has serious implications for attempts to derive the disc gas mass and other properties from line observations. This settling flux enhancement arises from changes to the vertical temperature and chemical structure, as settled dust grains allow stellar UV to penetrate deeper into the disc. The effect is strongest in lines which are formed in the warm gas in the inner disc (e.g. a

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Gas temperature structure for the preferred model (Col. 3 in Table 6), irradiated by an additional X-ray spectral component with $L_X = 10^{29.6}$ erg s$^{-1}$, as observed by Günther & Schmitt (2009) in HD 163296. This is in comparison to the right hand panel of Fig. 7, in which no X-ray spectral component is present.</figcaption></figure>

factor ~30 increase in the H$_2$ S(1) line), but the low excitation molecular lines are also affected, e.g. a factor ~2 increase in the $^{13}$CO 1–0 line.

We explore the effects of the observed UV variability in this object on the gas chemistry in our models, and conclude that this effect is not large enough to affect the observable line fluxes beyond the current range of instrumental uncertainty.

We also examine the effect of X-rays on the gas chemistry of our models, and find that while X-rays present a significant source of gas-heating in the disc surface layers, the observed X-ray luminosity of $L_X = 10^{29.6}$ erg s$^{-1}$ does not significantly alter the gas chemistry or line emission. Any effects are smaller than those expected as a result of the observed UV variability in this object.

It is difficult to reach any firm conclusions regarding the evolutionary state of the disc of HD 163296. There is some evidence to suggest that the disc is gas-depleted. This is in contrast to the result found for the Herbig Ae star HD 169142 by Meeus et al. (2010), where despite indications that the disc is transitional, it was found to be gas-rich. We note that there are uncertainties associated with our gas/dust values arising from uncertainties in the disc dust composition, and the associated opacity law. This is reflected by the factor of ~3 spread in the range of derived dust masses found in the literature for this object (Natta et al. 2004; Tannirkulam et al. 2008a; Mannings & Sargent 1997; Isella et al. 2007). All of our derived dust masses are within the range (5–17) $\times 10^{-4}$ M$_\odot$ from the literature. We also note the restrictions placed on our conclusions by the assumption of constant dust grain properties throughout the disc, and it has been suggested that the dust properties in discs should in general be variable with radius (Birnstiel et al. 2010; Guilloteau et al. 2011).

It would appear that the modelling of HD 163296 is a far from straightforward task, and it is difficult to fit it into the standard evolutionary picture developed over the past decade. We are unable to find a single disc model which fits perfectly the entire wealth of observational data for this object. As well as
<figure><img src="image_14.png" /></figure>

**Fig. 13.** Transient absorption components in 16 May 1998. The high resolution observed spectrum of several lines (solid lines) is plotted along with the synthetic photospheric spectrum (dashed lines; Kurucz 1993). The short vertical dashes mark the approximate location of the metallic-only transient absorptions at +50 km s⁻¹. A transient red-shifted component of velocity *v* ~ 50 km s⁻¹ is clearly seen in the metallic but not the hydrogen lines. In addition a blue-shifted absorption is identified in all lines for *v* ~ -100 km s⁻¹.

the gas/dust ratio, there is uncertainty regarding the disc flar-
ing, which itself is strongly tied in to the disc gas heating and
line emission, as well as being indicative of the disc’s evolution-
ary state. Evidence of in-falling material and the presence of a
bipolar outflow in HD 163296 seem to indicate that the star is
actively accreting material, typical of a young object with a gas-
rich disc. However, it has been suggested that the observed disc
geometry and flaring could indicate a star at a later stage of its
evolution. This would be consistent with the evidence for grain
growth from our modelling, with all the models requiring large
grains to fit the observations. There is also possible further evi-
dence for particle growth from high resolution optical spectra of
this object (see Appendix). This is also consistent with possible
evidence for dust settling and gas-depletion from this study, but
is hard to reconcile with evidence for substantial ongoing activ-
ity in the inner disc. It is clear that this object is at a fascinating
stage in its evolution.

**Appendix: Échelle spectra analysis**

Section 5 provides possible evidence of a low gas-to-dust ratio
in the HD 163296 circumstellar disc as compared to interstellar
abundances. This result is also supported in this Appendix from
the analysis of visible high resolution échelle spectra.

The EXPORT collaboration (Eiroa et al. 2001) obtained high
resolution échelle spectra (R = λ/Δλ = 48 000) in the range
3800–5900 Å during five nights in May and July 1998 using the

<figure><img src="image_15.png" /></figure>

**Fig. 14.** Transient absorption components in 16 July 1998. The high resolution observed spectrum of several lines (solid lines) is plotted along with the photospheric spectrum (dashed lines). No clear transient absorption is detected in the metallic lines. On the other hand a significant blue-shifted absorption is detected in all Balmer lines at *v* ~ -100 km s⁻¹.

Utrecht Echelle Spectrograph (UES) on the William Herschel Telescope, La Palma.

Significant spectral variability was found in the optical spec-
trum. The EXPORT data monitored those changes on time scales
of hours, days and months (Mendigutía et al. 2011). Almost any
strong photospheric line is variable to some extent. In particular,
transient absorption components (TACs) are observed in many
hydrogen (Balmer series) and metallic lines, such as Ca II K &
H, Na I D2 and D1 and the Fe II 4924, 5018 and 5169 Å triplet.
The Appendix summarizes the main results on the dynamics of
the circumstellar gas obtained from the TACS.

Many TACs observed are similar to those studied by Natta
et al. (2000) and Mora et al. (2002, 2004) in similar age
Herbig Ae stars. These authors propose magnetospheric accre-
tion of primordial gas as the most likely origin. That is, the TACs
in Herbig Ae stars typically trace the kinematics of material with
gas to dust ratios similar to that of the interstellar medium.

However, a distinctive feature of HD 163296 is the presence
of low velocity (v ~ 50 km s⁻¹) red-shifted TACs only observed
in refractory elements: Fe I, Fe II, Ti II and Sc II. These “metal-
lic” TACs have no clear counterpart neither in H nor Ca II and
Na I. They were observed during all nights except 28 July 1998.
These transient absorptions could be better explained in terms
of the model developed by Lagrange-Henri et al. (1988) for the
evaporation of solid bodies in the β Pic debris disc and subse-
quently postulated by Grinin et al. (1991, 1994) for accretion
discs around Herbig Ae/Be stars.

Figure 13 shows selected line profiles observed in 16 May
1998. Wavelengths have been converted to radial velocities in
the stellar reference system. The approximate location of the
<figure><img src="image_16.png" /></figure>

**Fig. 15.** Kinematics of the circumstellar gas. Features with similar radial velocity observed in different lines and different dates are assumed to share a common origin and grouped into events. Each event is represented by a box enclosing a number of points. Each point represents the weighted average of the radial velocity simultaneously measured for several lines. The error bar is the standard deviation. The weighted number of lines is given next to each point. Red boxes represent in-falling material, blue boxes out-flowing material and green boxes in-falling material only detected in metallic lines.

metallic-only transient absorptions at +50 km s⁻¹ is displayed with a short vertical dash. The underlying photosphere is represented by a synthetic Kurucz model computed with $T_{\text{eff}} = 9250$ K, $\log g_* = 4.07$ and [Fe/H] = +0.20, using the ATLAS9 and SYNTHE codes by Kurucz (1993). In this work we have used the GNU Linux version of the codes available online<fnref n="6" /> (Sbordone et al. 2004). It is clear that no Balmer line absorption counterpart can be identified, even for the higher $H\delta$ and $H\zeta$ lines in the series for which no significant core line emission is apparent. On the other hand, blue-shifted absorption with $v \sim -100$ km s⁻¹ is clearly seen both in metallic and hydrogen lines. The metallic only transient events are not an artifact generated by poor photosphere characterisation, as revealed by Fig. 14, which shows the spectra obtained in 28 July 1998, the only night where no metallic only events where identified.

Similar results were found for WW Vul by Mora et al. (2004), but not for the other four stars studied. In addition, both hydrogen-rich and hydrogen-less TACs are simultaneously detected in HD 163296 and WW Vul. If particle growth indeed explains the metallic-only events, that would suggest the circumstellar disc around HD 163296 (and incidentally WW Vul) has an evolutionary status between a pure primordial gas rich accretion disc and a gas-depleted second generation debris disc. That finding is fully consistent with the results obtained in Sect. 5 from the analysis of Herschel/PACS far-IR spectroscopy and ground-based mm CO observations.

## Kinematics of the circumstellar gas

A detailed study of the dynamics of the circumstellar gas around HD 163296 using the method presented in Mora et al. (2004) is out the scope of this paper. The main results are, however, summarized here.

The following lines were studied: the Balmer series ($H\beta$, $H\gamma$, $H\delta$, $H\epsilon$ and $H\zeta$), The Na I D and Ca II HK doublets, the Fe II 4924, 5018 and 5169 Å triplet, Fe I 4046 Å, Ti II 4444 and 4572 Å and Sc II 4247 Å. Three types of variability were always simultaneously observed: core line emission (most prominent in the Balmer, Na I and Ca II lines), transient red-shifted absorptions (associated to in-falling material, observed in all lines) and transient blue-shifted absorptions (associated to outflows, observed in all lines). HD 163296 was more active than the stars studied by Mora et al. (2002, 2004). In particular, the larger line core emission made the analysis more difficult.

The in-falling material always span a large range of velocities, from near zero to +300 km s⁻¹ for the red end of the higher velocity components. The stellar winds maximum velocities can also be high. A maximum of -400 km s⁻¹ was observed in 17 May 1998.

Figure 15 shows a summary of all the transient absorptions identified. Features with similar radial velocity observed in different lines are assumed to share a common origin. The average radial velocity and the standard deviation are displayed as points and error bars in the figure. The number of lines used to compute the average is given next to each point. The fainter lines are given half weight in the average (see Mora et al. 2004 for additional details). Points with similar velocity on nearby nights are also assumed to be produced by the same clump of circumstellar material and are grouped in “events”. Events are displayed in the figure with boxes enclosing a certain number of average points. Red boxes represent in-falling material, blue boxes outflows and green boxes in-falling material only detected in metallic lines. Metallic-only events have smaller radial velocities than hydrogen-rich in-falling material.

One of the most interesting events observed is #4, a strong wind observed in plenty of Balmer and metallic lines during July 1998. The most relevant feature is an almost constant

<footnote n="6">⁶ http://wwwuser.oat.ts.astro.it/atmos/</footnote>
acceleration rate of ~0.4 m/s². The average radial velocity starts at -120 km s⁻¹ on 28 July 1998 and grows up to -235 km s⁻¹ in 31 July 1998. The width of the absorption features has an almost constant value of 120 km s⁻¹, supporting the hypothesis of a common origin of all the transient absorptions in the event.

## References

Abrahamsson, E., Krems, R. V., & Dalgarno, A. 2007, ApJ, 654, 1171

Acke, B., Bouwman, J., Juhász, A., et al. 2010, ApJ, 718, 558

Aresu, G., Kamp, I., Meijerink, R., et al. 2011, A&A, 526, A163

Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087

Bell, K. L., Berrington, K. A., & Thomas, M. R. J. 1998, MNRAS, 293, L83

Benisty, M., Natta, A., Isella, A., et al. 2010, A&A, 511, A74

Bessell, M. S. 1979. PASP, 91, 589

Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79

Bouwman, J., de Koter, A., van den Ancker, M. E., & Waters, L. B. F. M. 2000, A&A, 360, 213

Brott, I., & Hauschildt, P. H. 2005, in The Three-Dimensional Universe with Gaia, ed. C. Turon, K. S. O'Flaherty, & M. A. C. Perryman, ESA SP, 576, 565

Chamaud, G., Levy, B., Millie, P., Roueff, E., & Tran Minh, F. 1980, J. Phys. B Atomic Mol. Phys., 13, 4205

Cox, A. N., ed. 1999, Allen's Astrophysical Quantities, 4th ed. (New York: Springer-Verlag)

Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362

D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411

de Winter, D., van den Ancker, M. E., Maira, A., et al. 2001, A&A, 380, 609

Devine, D., Grady, C. A., Kimble, R. A., et al. 2000, ApJ, 542, L115

Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503

Doucet, C., Pantin, E., Lagage, P. O., & Dullemond, C. P. 2006, A&A, 460, 117

Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269

Dubernet, M.-L., & Grosjean, A. 2002, A&A, 390, 793

Eiroa, C., Garzón, F., Alberdi, A., et al. 2001, A&A, 365, 110

Eisner, J. A., Graham, J. R., Akeson, R. L., & Najita, J. 2009, ApJ, 692, 309

Faure, A., Crimier, N., Ceccarelli, C., et al. 2007, A&A, 472, 1029

Fitzpatrick, E. L. 1999, PASP, 111, 63

Flower, D. R. 2001, J. Phys. B Atomic Mol. Phys., 34, 2731

Flower, D. R. & Launay, J. M. 1977, J. Phys. B Atomic Mol. Phys., 10, 3673

Geers, V. C., Augereau, J., Pontoppidan, K. M., et al. 2006, A&A, 459, 545

Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ, 544, 895

Green, S., Maluendes, S., & McLean, A. D. 1993, ApJS, 85, 181

Grinin, V. P., Kiselev, N. N., Chernova, G. P., Minikulov, N. K., & Voshchinnikov, N. V. 1991, Ap&SS, 186, 283

Grinin, V. P., The, P. S., de Winter, D., et al. 1994, A&A, 292, 165

Guilloteau, S., Dutrey, A., Piétu, V., & Boehler, Y. 2011, A&A, 529, A105

Günther, H. M., & Schmitt, J. H. M. M. 2009, A&A, 494, 1041

Henning, T., Begemann, B., Mutschke, H., & Dorschner, J. 1995, A&AS, 112, 143

Hughes, A. M., Wilner, D. J., Qi, C., & Hogerheijde, M. R. 2008, ApJ, 678, 1119

Hughes, A. M., Wilner, D. J., Andrews, S. M., Qi, C., & Hogerheijde, M. R. 2011, ApJ, 727, 85

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Jager, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291

Jankowski, P., & Szalewicz, K. 2005, J. Chem. Phys., 123, 104301

Jaquet, R., Staemmler, V., Smith, M. D., & Flower, D. R. 1992, J. Phys. B Atomic Mol. Phys., 25, 285

Jonkheid, B., Dullemond, C. P., Hogerheijde, M. R., & van Dishoeck, E. F. 2007, A&A, 463, 203

Kamp, I., & Bertoldi, F. 2000, A&A, 353, 276

Kamp, I., Tilling, I., Woitke, P., Thi, W., & Hogerheijde, M. 2010, A&A, 510, A18

Kamp, I., Woitke, P., Pinte, C., et al. 2011, A&A, 532, A85

Kurucz, R. L. 1993, VizieR Online Data Catalog, 6039, 0

Lagrange-Henri, A. M., Vidal-Madjar, A., & Ferlet, R. 1988, A&A, 190, 275

Launay, J.-M., & Roueff, E. 1977, J. Phys. B Atomic Mol. Phys., 10, 879

Mannings, V. 1994, MNRAS, 271, 587

Mannings, V., & Sargent, A. I. 1997, ApJ, 490, 792

Martin-Zaïdi, C., Augereau, J., Ménard, F., et al. 2010, A&A, 516, A110

Mathews, G. S., Dent, W. R. F., Williams, J. P., et al. 2010, A&A, 518, L127

Meeus, G., Waters, L. B. F. M., Bouwman, J., et al. 2001, A&A, 365, 476

Meeus, G., Pinte, C., Woitke, P., et al. 2010, A&A, 518, L124

Meijer, J., Dominik, C., de Koter, A., et al. 2008, A&A, 492, 451

Meijerink, R., Pontoppidan, K. M., Blake, G. A., Poelman, D. R., & Dullemond, C. P. 2009, ApJ, 704, 1471

Mendigutía, I., Eiroa, C., Montesinos, B., et al. 2011, A&A, 529, A34

Monnier, J. D., Berger, J., Millan-Gabet, R., et al. 2006, ApJ, 647, 444

Montesinos, B., Eiroa, C., Mora, A., & Merín, B. 2009, A&A, 495, 901

Mora, A., Natta, A., Eiroa, C., et al. 2002, A&A, 393, 259

Mora, A., Eiroa, C., Natta, A., et al. 2004, A&A, 419, 225

Natta, A., Grinin, V. P., & Tambovtsseva, L. V. 2000, ApJ, 542, 421

Natta, A., Testi, L., Neri, R., Shepherd, D. S., & Wilner, D. J. 2004, A&A, 416, 179

Oudmaijer, R. D., Palacios, J., Eiroa, C., et al. 2001, A&A, 379, 564

Panić, O., & Hogerheijde, M. R. 2009, A&A, 508, 707

Piétu, V., Guilloteau, S., & Dutrey, A. 2005, A&A, 443, 945

Piétu, V., Dutrey, A., & Guilloteau, S. 2007, A&A, 467, 163

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pinte, C., Padgett, D. L., Ménard, F., et al. 2008, A&A, 489, 633

Posch, T., Kerschbaum, F., Fabian, D., et al. 2003, ApJS, 149, 437

Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362

Renard, S., Malbet, F., Benisty, M., Thiébaut, E., & Berger, J. 2010, A&A, 519, A26

Sandell, G., Weintraub, D. A., & Hamidouche, M. 2011, ApJ, 727, 26

Sbordone, L., Bonifacio, P., Castelli, F., & Kurucz, R. L. 2004, Mem. Soc. Astron. Ital. Suppl., 5, 93

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Servoin, J. L., & Piriou, B. 1973, Physica Status Solidi (b), 55, 677

Sitko, M. L., Carpenter, W. J., Kimes, R. L., et al. 2008, ApJ, 678, 1070

Tannirkulam, A., Monnier, J. D., Harries, T. J., et al. 2008a, ApJ, 689, 513

Tannirkulam, A., Monnier, J. D., Millan-Gabet, R., et al. 2008b, ApJ, 677, L51

Tennyson, J., Zobov, N. F., Williamson, R., Polyansky, O. L., & Bernath, P. F. 2001, J. Phys. Chem. Ref. Data, 30, 735

Thi, W. F., van Dishoeck, E. F., Blake, G. A., et al. 2001, ApJ, 561, 1074

Tielens, A. G. G. M. 2008, ARA&A, 46, 289

Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389

van den Ancker, M. E., Bouwman, J., Wesselius, P. R., et al. 2000, A&A, 357, 325

van Leeuwen, F., ed. 2007, Hipparcos, the New Reduction of the Raw Data, Astrophys. Space Sci. Lib., 350

Warren, S. G., & Brandt, R. E. 2008, J. Geophys. Res., 113, D14220

Wassell, E. J., Grady, C. A., Woodgate, B., Kimble, R. A., & Bruhweiler, F. C. 2006, ApJ, 650, 985

Wernli, M., Valiron, P., Faure, A., et al. 2006, A&A, 446, 367

Wilson, N. J., & Bell, K. L. 2002, MNRAS, 337, 1027

Wisniewski, J. P., Clampin, M., Grady, C. A., et al. 2008, ApJ, 682, 548

Woitke, P., Kamp, I., & Thi, W. 2009, A&A, 501, 383

Woitke, P., Pinte, C., Tilling, I., et al. 2010, MNRAS, 405, L26

Woitke, P., Riaz, B., Duchêne, G., et al. 2011, A&A, 534, A44

Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197

Wu, C., Ake, T. B., Boggess, A., et al. 1983, NASA IUE Newsl., No. 22, 2+324 pp., 22

Wu, C., Reichert, G. A., Ake, T. B., et al. 1992, NASA Reference Publication, 1285

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062

Yi, S., Demarque, P., Kim, Y.-C., et al. 2001, ApJS, 136, 417

Zuckerman, B., Forveille, T., & Kastner, J. H. 1995, Nature, 373, 494