# Nature of the gas and dust around 51 Ophiuchi. Modelling continuum and Herschel line observations

W. F. Thi, F. Ménard, G. Meeus, A. Carmona, P. Riviere-Marichalar, J. -C.

Augereau, I. Kamp, P. Woitke, C. Pinte, I. Mendigutía, et al.

► To cite this version:

W. F. Thi, F. Ménard, G. Meeus, A. Carmona, P. Riviere-Marichalar, et al. Nature of the gas and dust around 51 Ophiuchi. Modelling continuum and Herschel line observations. Astronomy and Astrophysics - A&A, 2013, 557, 10.1051/0004-6361/201221002. insu-03620913

HAL Id: insu-03620913

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

Submitted on 27 Mar 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Nature of the gas and dust around 51 Ophiuchi*

## Modelling continuum and Herschel line observations

W. F. Thi<fnref n="1" />, F. Ménard<fnref n="2" />,<fnref n="1" />, G. Meeus<fnref n="3" />, A. Carmona<fnref n="1" />, P. Riviere-Marichalar<fnref n="4" />,<fnref n="3" />, J.-C. Augereau<fnref n="1" />, I. Kamp<fnref n="5" />, P. Woitke<fnref n="6" />, C. Pinte<fnref n="1" />, I. Mendigutía<fnref n="7" />, C. Eiroa<fnref n="3" />, B. Montesinos<fnref n="3" />, S. Britain<fnref n="7" />, and W. Dent<fnref n="8" />

Received 24 December 2012 / Accepted 16 April 2013

## ABSTRACT

**Context.** Circumstellar disc evolution is paramount for the understanding of planet formation. The gas in protoplanetary discs large program (GASPS) aims at determining the circumstellar gas and solid mass around ~250 pre-main-sequence Herbig Ae and T Tauri stars.

**Aims.** We aim to understand the origin and nature of the circumstellar matter orbiting 51 Oph, a young (<1 Myr) luminous B9.5 star.

**Methods.** We obtained continuum and line observations with the PACS instrument on board the Herschel Space Observatory and continuum data at 1.2 mm with the IRAM 30 m telescope. The spectral energy distribution and line fluxes were modelled using the physico-chemo radiative transfer code ProDiMo to constrain the gas and solid mass of the disc around 51 Oph. The disc vertical hydrostatic structure was computed self-consistently together with the gas thermal balance.

**Results.** We detected a strong emission by atomic oxygen [O I] at 63 microns using the Herschel Space Observatory. The [O I] emission at 145 microns, the [C II] emission at 158 microns, the high-J CO emissions, and the warm water emissions were not detected. Continuum emission was detected at 1.2 mm. The continuum from the near- to the far-infrared and the [O I] emission are well explained by the emission from a compact ($R_{out}$ = 10–15 AU) hydrostatic disc model with a gas mass of $5 \times 10^{-6} M_{\odot}$, 100 times that of the solid mass. However, this model fails to match the continuum millimeter flux, which hints at a cold outer disc with a mass in solids of ~$10^{-6} M_{\odot}$ or free-free emission from a photoevaporative disc wind. This outer disc can either be devoid of gas and/or is too cold to emit in the [O I] line. A very flat extended disc model ($R_{out}$ = 400 AU) with a fixed vertical structure and dust settling matches all photometric points and most of the [O I] flux.

**Conclusions.** The observations can be explained by an extended flat disc where dust grains have settled. However, a flat gas disc cannot be reproduced by hydrostatic disc models. The low mass of the 51 Oph inner disc in gas and dust may be explained either by the fast dissipation of an initial massive disc or by a very small initial disc mass.

**Key words.** astrochemistry – protoplanetary disks – line: identification

## 1. Introduction

The large variety of extrasolar planetary systems requires a theoretical and observational understanding of protoplanetary discs where the exoplanets are formed. Observing protoplanetary discs at different stages of evolution makes it possible to test the various models of disc dissipation, planet-formation, and planet-migration. Discs are characterized by their gas and solid contents and their geometry. The two main categories of discs are the young massive and flaring ones ($M_{disc} = 10^{-2}-10^{-3} M_{\odot}$ in gas and dust, typical size of 100–500 AU), characterized by a strong IR excess in their spectral energy distribution (SED); and the older debris discs, where there is no gas left and where

the small dust grains result from the collisions of planetesimals. The dust discs in the latter category are geometrically flat and the SED shows weak excess in the far-IR. The passage from one to the other categories seems rapid because there are very few of the so-called transition discs. The definition of a transition disc is not commonly accepted yet, but discs with a lack of near-IR excess because of an inner hole caused by a planet, but that have strong far-IR excess can be considered as transitional discs. The rapidly rotating southern B9.5IV-B9.5V star 51 Oph (a.k.a. HR 6591, HD 158643, SAO 185470, IRAS 17283-2355) located at ~131 pc is a rare example of a pre-main-sequence star that is surrounded by a disc where most or all of the primordial gas has probably dissipated. The SED of the 51 Oph disc suggests a new category of disc. It does show a near-IR excess, but the main property of the SED is a drop in the continuum flux beyond 20 microns. No observations existed beyond 60 microns.

<footnote n="1">¹ UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique (IPAG), UMR 5274, 38041 Grenoble, France
e-mail: Wing-Fai.Thi@obs.ujf-grenoble.fr</footnote>

<footnote n="2">² UMI - LFCA, CNRS/INSU France, and Dept. de Astronomia y Obs. Astronomico Nacional, Universidad de Chile, Casilla 36-D, Correo Central, Santiago (UMI 3386), Chile</footnote>

<footnote n="3">³ Dep. de Física Teórica, Fac. de Ciencias, UAM Campus Cantoblanco, 28049 Madrid, Spain</footnote>

<footnote n="4">⁴ Centro de Astrobiologa Depto. Astrofsica (CSICINTA), ESAC Campus, PO Box 78, 28691 Villanueva de la Canada, Spain</footnote>

<footnote n="5">⁵ Kapteyn Astronomical Institute, PO Box 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="6">⁶ SUPA, School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK</footnote>

<footnote n="7">⁷ Department of Physics and Astronomy, Clemson University, USA</footnote>

<footnote n="8">⁸ ALMA, Avda Apoquindo 3846, Piso 19, Edicio Alsacia, Las Condes, Santiago, Chile</footnote>

<footnote>* Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>
The absence of strong far-IR continuum dust emission points to a small amount of dust grains and a non-flaring disc. The 10-microns feature is consistent with emission by amorphous silicates, and there are almost no features that can be attributed to polycyclic aromatic hydrocarbon (PAH) emissions (Keller et al. 2008). Stark et al. (2009) proposed a two-optically-thin-disc model to explain the combination of MIDI-VLT, Keck nulling, and Spitzer-IRS data. Their outer dust disc extends to 1200 AU and is flat.

The evolutionary status of 51 Oph is debated. The low dust mass favours the idea that 51 Oph is a young debris disc, as reported in Stark et al. (2009). A debris disc would have a lower gas-to-dust mass ratio than a young disc. One of the key questions is therefore how much gas (atomic + molecular) is still present in the 51 Oph disc? The total atomic + molecular gas mass of the 51 Oph disc has never been determined. In addition to the mass, the nature of the gas in the 51 Oph disc is also disputed. Is the gas atomic or molecular? A predominately atomic disc would favour the scenario where most of the gas comes from the evaporation of infalling comets without primordial gas left, whereas significant amounts of molecular gas ($H_2$) as traced by rotational emissions by CO and its isotopologues would indicate that the gas is a leftover from planet-formation and/or disc dissipation. Absorption studies in the ultraviolet with the FUSE telescope (Grady & Silvis 1993; Lecavelier Des Etangs et al. 1997; Roberge et al. 2002) suggest infalling atomic gas and no molecular gas, similar to the case of the $\beta$ Pic debris disc. On the other hand, spectra in the infrared show ro-vibrational emission lines from molecules such as CO, $CO_2$, and $H_2O$ (e.g. van den Ancker et al. 2001). The presence of molecules suggests that the gas may be dense and warm enough so that efficient gas-phase neutral-neutral reactions can produce molecules such as water. Neutral carbon (C I) has been detected in absorption with the FUSE telescope, and the C I/dust column density ratio is low, similar to the value found in the line-of-sight of $\beta$ Pic. Hot carbon monoxide (CO bandhead emission) has been observed by several groups and was assigned to the emission from a hot gas in a disc seen close to edge-on (Thi et al. 2005; Berthoud et al. 2007; Tatulli et al. 2008), but no rotational emission from cold ($T = 10-30$ K) CO in the sub-millimeter range is available. In this paper, we present new observations in the millimeter and far-IR obtained with the IRAM 30 m telescope and the Herschel-PACS Space Observatory in Sects. 2 and 3. The 51 Oph Herschel data is part of the HerbigAe and debris disc survey from the GASPS large open time programme (P.I. Dent; Dent et al. 2013). The new observations (continuum and line) combined with archival data are modelled to derive an estimate of the disc dust and gas masses in Sect. 4. In the discussion, we consider the nature and the different possible origins of the 51 Oph disc (Sect. 6).

## 2. Observations and data reduction

The source 51 Oph was observed for continuum emission at 1.2 mm using the MAMBO2 bolometer array (Kreysa et al. 1998) on the IRAM 30 m telescope at Pico Veleta, Spain. Observations were conducted during the bolometer pool in 2008 (Nov. 13–18, 23). Zenith opacity for our observations was typically ~0.2–0.3. Observations were first carried out to a target $1\sigma$ sensitivity of 1 mJy (20 min on-source) and then repeated on two subsequent nights to reach 5 sigma detection. A total of 60 min on source was achieved, in an ON-OFF pattern of 1 minute on target followed by 1 min on sky, with a throw of 32". Flux calibration was carried out using Mars, Saturn, or Jupiter (depending on availability), and local pointing and

<table><caption>**Table 1.** Lines and photometric data obtained by Herschel-PACS.</caption><thead><tr><th>Line</th><th>λ<br>(μm)</th><th>Cont. flux<br>(Jy)</th><th>Line flux<br>(10<sup>-18</sup> W m<sup>-2</sup>)</th></tr></thead><tbody><tr><td>[O I]<sup>3</sup>P<sub>1</sub> → <sup>3</sup>P<sub>2</sub></td><td>63.183</td><td>1.09 ± 0.19</td><td>53.3 ± 2.5</td></tr><tr><td>[O I]<sup>3</sup>P<sub>0</sub> → <sup>3</sup>P<sub>1</sub></td><td>145.525</td><td>0.11 ± 0.10</td><td>&lt;5.1</td></tr><tr><td>[C II]<sup>2</sup>P<sub>3/2</sub> → <sup>2</sup>P<sub>1/2</sub></td><td>157.74</td><td>0.07 ± 0.14</td><td>&lt;6.8</td></tr><tr><td>CO J = 36 → 35</td><td>72.85</td><td>0.71 ± 0.20</td><td>&lt;11.6</td></tr><tr><td>CO J = 33 → 32</td><td>76.36</td><td>0.52 ± 0.28</td><td>&lt;16.1</td></tr><tr><td>CO J = 29 → 28</td><td>90.16</td><td>0.50 ± 0.19</td><td>&lt;8.0</td></tr><tr><td>CO J = 18 → 17</td><td>144.78</td><td>0.12 ± 0.07</td><td>&lt;4.2</td></tr><tr><td>OH 1/2 – 3/2</td><td>79.11</td><td>0.54 ± 0.19</td><td>&lt;11.4</td></tr><tr><td>OH 1/2 – 3/2</td><td>79.18</td><td>0.50 ± 0.19</td><td>&lt;11.3</td></tr><tr><td>o-H<sub>2</sub>O 8<sub>18</sub>–7<sub>07</sub></td><td>63.32</td><td>1.08 ± 0.15</td><td>&lt;8.3</td></tr><tr><td>o-H<sub>2</sub>O 7<sub>07</sub>–6<sub>16</sub></td><td>71.94</td><td>0.71 ± 0.23</td><td>&lt;16.1</td></tr><tr><td>p-H<sub>2</sub>O 3<sub>22</sub>–2<sub>11</sub></td><td>89.99</td><td>0.45 ± 0.24</td><td>&lt;10.3</td></tr><tr><td>p-H<sub>2</sub>O 4<sub>13</sub>–3<sub>22</sub></td><td>144.51</td><td>0.12 ± 0.08</td><td>&lt;4.3</td></tr><tr><td>p-H<sub>2</sub>O 3<sub>31</sub>–4<sub>04</sub></td><td>158.31</td><td>0.09 ± 0.15</td><td>&lt;7.3</td></tr><tr><td>o-H<sub>2</sub>O 2<sub>12</sub>–1<sub>01</sub></td><td>179.52</td><td>0.24 ± 0.19</td><td>&lt;7.0</td></tr><tr><td>o-H<sub>2</sub>O 2<sub>21</sub>–2<sub>12</sub></td><td>180.42</td><td>0.13 ± 0.16</td><td>&lt;5.7</td></tr><tr><td>CH<sup>+</sup> J = 5–4</td><td>72.14</td><td>0.64 ± 0.3</td><td>&lt;19.3</td></tr><tr><td>CH<sup>+</sup> J = 4–3</td><td>90.01</td><td>0.45 ± 0.24</td><td>&lt;10.3</td></tr><tr><td colspan="4" style="text-align: center;">Herschel-PACS continuum data</td></tr><tr><td>...</td><td>70.0</td><td>0.873 ± 0.023</td><td>...</td></tr><tr><td>...</td><td>100.0</td><td>0.347 ± 0.010</td><td>...</td></tr><tr><td>...</td><td>160.0</td><td>0.090 ± 0.006</td><td>...</td></tr></tbody></table>

**Notes.** The errors are $1\sigma$ and the upper limits are $3\sigma$. The calibration error adds an extra ~10% uncertainty.

secondary flux calibration was carried out using IRAS 16293-2422B. The data were reduced using the facility reduction software, MOPSIC<fnref n="1" />.

We obtained photometric and line observations with the PACS instrument (Poglitsch et al. 2010) on board the Herschel Space Observatory (Pilbratt et al. 2010). We obtained PACS scan map observations at 70, 100, and 160 microns. At 100 microns, 51 Oph was observed at two angles (70 and 110 degrees) to improve the noise suppression, while at 70 microns, we only observed one scan direction. Both the 70 and 100 microns scans include the red band at 160 microns. The PACS photometric data were reduced with the mini scanmap pipeline in HIPE version 8.1.0, (calTree version 32), and we constructed mosaics of the two scan maps at 100 and 160 microns. The continuum flux was extracted using an aperture of 12" (sky between 30 and 35") at 70 microns, of 15" (sky between 33 and 38") at 100 microns, of 20" (sky between 40 and 50") at 160 microns. We obtained fluxes of 873 ± 23 mJy at 70, 347 ± 10 mJy at 100 and 90 ± 6 mJy at 160 microns (see Table 1). The calibration errors haven been included. The spectroscopic data were reduced with the official version 8.0.1 of the Herschel Interactive Processing Environment (HIPE; Ott 2010), using standard tasks provided in HIPE. These include bad pixel flagging; chop on/off subtraction; spectral response function division; rebinning with oversample =2 and upsample =1, corresponding to the native resolution of the instrument; spectral flatfielding and finally averaging of the two nod positions. To conserve the best signal and to not introduce additional noise, we only extracted the central spaxel and corrected for the flux loss with an aperture correction provided by the PACS instrument team. We extracted the fluxes

<footnote n="1">¹ http://www.iram.es/IRAMES/mainWiki/CookbookMopsic</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Herschel spectra towards 51 Oph. The [O I] line is spectrally unresolved. The line peak wavelengths are shown.</figcaption></figure>

of the detected lines using a Gaussian fit to the emission lines with a first-order polynomial to the continuum, using the rms on the continuum (excluding the line) to derive a $1\sigma$ error on the line by integrating a Gaussian with height equal to the continuum RMS and width of the instrumental FWHM. In case of a non-detection, we give a $3\sigma$ upper limit, also calculated on the continuum RMS. The measured line fluxes are listed in Table 1.

Pipeline-processed FEROS<fnref n="2" /> ESO 2.2 m high-resolution ($R \sim 45\,000$) optical spectra (4000–9000 Å) were obtained from the Advanced Data Products ESO archive<fnref n="3" />. Details of the FEROS data reduction pipeline can be found in the ESO HARPS and FEROS pipeline web-site<fnref n="4" />. The spectra are displayed at the observed wavelength and no correction was performed for the radial velocity or motion of the earth around the sun. Finally, we complemented the new data with archival photometry.

### 3. Observational results

We detected the spatially unresolved (on a scale of 10′′) emission at 1.2 mm with a flux of $5 \pm 1$ mJy (1$\sigma$ error) towards 51 Oph using the MAMBO bolometer at the IRAM 30 m telescope, consistent with the upper limit of 42 mJy obtained at 1.1 mm by

Sylvester et al. (1996). The continuum was detected by PACS in all settings. The total SED from the optical to the millimeter is shown in Fig. 1. Fajardo-Acosta et al. (1993) detected the 10 µm silicate feature at low spectral resolution. The feature was later confirmed by van den Ancker et al. (2001). The millimeter flux translates into a dust grain (with radii <1 mm) mass of (1–2) × 10⁻⁶ M☉ assuming $T_{\text{dust}} = 50$ K and the standard dust opacity $\kappa_v = 2.0(v/v_0)^\beta g^{-1} \text{cm}^2$ where $v_0 = 230.769$ GHz and $\beta = 0.6$ (Beckwith & Sargent 1991). This derived value is an upper limit to the dust mass since we cannot estimate the contribution from free-free emission at 1.2 mm due to the lack of data at longer wavelengths. Alternatively, pebbles with size greater than 1 cm will emit as blackbodies in the millimeter range. Assuming a gas-to-dust mass ratio of 100, the disc gas mass is (1–2) × 10⁻<fnref n="4" /> M☉.

We complemented the SED with archival photometric data and the Spitzer-IRS spectrum (see Fig. 3). We fitted the SED with a Kurucz stellar model in addition to a modified blackbody of the form

$$vF_v = B_v(T)(\lambda/\lambda_0)^{-\beta-1}, \quad (1)$$

where $B_v(T)$ is the Planck function at the temperature $T$. In the formula, we assumed that the grains will emit with $\kappa = \kappa_0 g^{-1} \text{cm}^{-1}$ for $\lambda < \lambda_0$ and with $\kappa_v = \kappa_0(\lambda/\lambda_0)^{-\beta} g^{-1} \text{cm}^2$ for $\lambda \ge \lambda_0$. The value of $\lambda_0$ is set by the criterion $\lambda_0 > 2\pi\langle a \rangle$, where $\langle a \rangle$ is the mean grain radius. We adopted $\lambda_0 = 100$ µm ($\langle a \rangle = 16$ µm). The fits to the SED for three values of $\beta$ ($\beta = 1, 0, -0.3$) are shown in Fig. 3. A value for $\beta$ of 1 is typical for grains

<footnote n="2">² The Fiber-fed Extended Range Optical Spectrograph, Kaufer et al. (1999), The Messenger 95, 8.</footnote>

<footnote n="3">³ http://archive.eso.org/</footnote>

<footnote n="4">⁴ http://www.eso.org/sci/facilities/lasilla/instruments/harps/tools/software.html/</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Normalized optical spectrum of 51 Oph around the Hα emission obtained by FEROS.</figcaption></figure>

found in protoplanetary discs (Draine 2006). A pure blackbody emission has $\beta = 0$. Alternatively, a grain size powerlaw distribution with an index $p = 3$ will result in a quasi-constant $\kappa$ with wavelength (Draine 2006). The model with the standard value for protoplanetary discs $\beta = 1$ fails to match the millimeter flux. A pure blackbody emission also underpredicts the millimeter flux, but the prediction is still consistent with the observations within $3\sigma$. The best fit is obtained with $\beta = -0.3$. A value for $\beta$ of $-0.3$ is unphysical and points to the emission by an extra component, either a cold outer disc/envelope or free-free emission.

In addition, far-IR spectral line observations were obtained using the Herschel Space Observatory. The only detected line is the [O I] emission at 63 microns. The [O I] emission is confined to the central spaxels. This line is always the strongest in the two [O I] and [C II] surveys of the GASPS project (Meeus et al. 2012). No high-J CO or water lines were detected. The line flux and three sigma upper limits are summarized in Table 1. Finally, we show the archival Hα line profile in Fig. 2. The profile shows a double-peaked structure on the top of the stellar Hα absorption. The origin of this profile is discussed in Sect. 4.

## 4. SED and gas modelling

We modelled the continuum and gas emission with the codes MCFOST and PRODIMO. MCFOST is a continuum and NLTE line Monte Carlo radiative transfer code that can handle 2D and 3D geometries (Pinte et al. 2006, 2009). PRODIMO is a radiative chemo-physical code that solves the continuum radiative transfer, the gas heating and cooling balance, the gas chemistry, and the disc vertical hydrostatic structure self-consistently with the gas pressure (Woitke et al. 2009). New features of the codes are described in Kamp et al. (2010), Thi et al. (2011), and Woitke et al. (2011). The 2D continuum radiative transfer module in PRODIMO and the code MCFOST has been benchmarked against other codes (Pinte et al. 2009).

The modelling strategy is twofold. First we assumed a fixed vertical structure disc with dust settling. The disc extends from 1.5 AU to 400 AU as suggested by the model of Stark et al. (2009). We employed an automatic fit procedure based on the Simplex optimization algorithm. The free parameters of the model (model a) are the total disc gas+dust mass, the surface density profile index, the gas-to-dust ratio, the disc gas scale-height, and the flaring index. In model b, the free parameters are

<table><caption>**Table 2.** Common disc parameters.</caption><tr><td>Stellar mass</td><td>M<sub>s</sub></td><td>3.8 M<sub>☉</sub></td></tr><tr><td>Stellar luminosity</td><td>L<sub>s</sub></td><td>260 L<sub>☉</sub></td></tr><tr><td>Effective temperature</td><td>T<sub>eff</sub></td><td>10 000 K</td></tr><tr><td>Distance</td><td>d</td><td>130 pc</td></tr><tr><td>Disc inclination</td><td>i</td><td>80° (0° face-on)</td></tr><tr><td>Disc inner radius</td><td>R<sub>in</sub></td><td>1.5 or 2.0 AU</td></tr><tr><td>Disc outer radius</td><td>R<sub>out</sub></td><td>15 or 400 AU</td></tr><tr><td>Vert. column density index</td><td>ε</td><td>1</td></tr><tr><td>Inner rim soft edge</td><td></td><td>on</td></tr><tr><td>Gas-to-dust mass ratio</td><td>δ</td><td>100</td></tr><tr><td>Dust grain material mass density</td><td>ρ<sub>dust</sub></td><td>3.5 g cm<sup>-3</sup></td></tr><tr><td>Minimum dust particle size</td><td>a<sub>min</sub></td><td>0.1 µm</td></tr><tr><td>Maximum dust particle size</td><td>a<sub>max</sub></td><td>5000 µm</td></tr><tr><td>Composition</td><td></td><td>ISM silicate</td></tr><tr><td>Dust size distribution power law</td><td>p</td><td>3.5</td></tr><tr><td>H<sub>2</sub> cosmic ray ionization rate</td><td>ζ<sub>CR</sub></td><td>1.7 × 10<sup>-17</sup> s<sup>-1</sup></td></tr><tr><td>ISM UV field w.r.t. Draine field</td><td>χ</td><td>1.0</td></tr><tr><td>Abundance of PAHs relative to ISM</td><td>f<sub>PAH</sub></td><td>10<sup>-3</sup></td></tr><tr><td>α viscosity parameter</td><td>α</td><td>0.0</td></tr><tr><td>Turbulence width</td><td>δv</td><td>0.15 km s<sup>-1</sup></td></tr><tr><td>Total disc mass</td><td>M<sub>disc</sub></td><td>5 × 10<sup>-6</sup>–1 × 10<sup>-4</sup> M<sub>☉</sub></td></tr><tr><td>Dust settling power law</td><td>s</td><td>0, 1.0</td></tr><tr><td>Minimum grain radius for settling</td><td>a<sub>s</sub></td><td>0.5 µm</td></tr></table>

**Notes.** Parameters are defined in Woitke et al. (2009).

the same except for the gas-to-dust ratio, which is fixed at 100. The dust scale height including settling is

$$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 fixed parameters: $s$ is the power-law index and $a_s$ is the minimum radius for a grain to be affected by settling. The settling law is described in details in Woitke et al. (2009) and has previously been used in modelling the HD 163296 disc (Tilling et al. 2012). We adopted a minimum dust radius for settling $a_s$ of 0.5 µm and $s = 1$. The disc density and gas-to-dust ratio structure are shown in Fig. 4.

Second, we relaxed the assumption of fixed disc vertical structure by letting the PRODIMO code calculate the gas scale-height that would be consistent with the local gas pressure as determined by the gas heating and cooling balance.

The disc dust and gas parameters are gathered in Table 2. The stellar parameters are taken from Montesinos et al. (2009). The input stellar spectrum is taken from the Kurucz database of models. The disc is seen close to edge-on between 80°–90° (Thi et al. 2005; Berthoud et al. 2007; Tatulli et al. 2008). The optical photometry shows that the star is seen without strong circumstellar extinction, therefore the dust disc has to be flat. On the other hand, the gas disc should be more extended vertically such that absorption lines can be detected (Grady & Silvis 1993; Lecavelier Des Etangs et al. 1997; Roberge et al. 2002). Therefore we tested disc models where the gas and dust are well mixed and where the dust grains have settled towards the midplane.

The accretion rate of the 51 Oph disc is disputed. Assuming a magnetospheric origin, Garcia Lopez et al. (2006) and Brittain et al. (2007) estimated a gas accretion rate of (1–2) × 10<sup>-7</sup> M<sub>☉</sub> yr<sup>-1</sup> from the hydrogen Bry line. However, Mendigutía et al. (2011) found that the Hα 10% peak width for 51 Oph is much lower than the expected minimum value required by the magnetospheric accretion model (Fig. 2). Therefore, they argued that the inner disc of 51 Oph is decoupled from the star and that the mass accretion rate of 10<sup>-7</sup> M<sub>☉</sub> yr<sup>-1</sup> is an upper
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Modified blackbody fits to the 51 Oph SED. The error bars on the photometric points are $3\sigma$. The orange line is the *Spitzer*-IRS spectrum taken from the NASA archive. The best fit with $\beta = -0.3$ is unphysical.</figcaption></figure>

limit. In a subsequent paper, Mendigutía et al. (2012) derived an empirical relation between the disc mass $M_{\text{disc}}$ in $M_{\odot}$ and the disc mass accretion $\dot{M}_{\text{acc}}$ in $M_{\odot} \text{yr}^{-1}$ valid for both T Tauri and Herbig Ae stars

$$ \log \dot{M}_{\text{acc}} = 1.1(\pm 0.3) \times \log M_{\text{disc}} - 5.0(\pm 0.5). \quad (3) $$

Using this relation and upper limit to the gas disc mass of $2 \times 10^{-4} M_{\odot}$ derived from the 1.2 mm flux assuming a gas-to-dust mass ratio of 100, we obtained an upper limit to the mass accretion rate of $\sim 8.5 \times 10^{-10} M_{\odot} \text{yr}^{-1}$. In this paper, we modelled the disc around 51 Oph as passive ($\alpha$ viscosity parameter is 0 in Table 2). The $H\alpha$ emission most likely does not trace magnetospheric accretion. The origin of the double-peaked $H\alpha$ emission profile can be ascribed to a stellar/disc wind (Kraus et al. 2008). A low turbulent width of $0.15 \text{ km s}^{-1}$ was assumed.

The inner disc radius is consistent with the location of the sublimation radius for a grain size distribution that follows a power-law with index 3.5 ($a_{\text{min}} = 0.1 \mu\text{m}$ and $a_{\text{max}} = 2000 \mu\text{m}$). The dust grains are composed of pure amorphous silicate (astronomical silicates, Laor & Draine 1993). The choice of amorphous silicate is motivated by the absence of strong crystalline silicate features in the infrared (Keller et al. 2008). The dust opacities computed using the Mie theory for compact spherical particles are shown in Fig. 7. Except from 10 and 100 $\mu\text{m}$, the opacity is dominated by the scattering term.

Even a perfect fit to the SED alone does not allow one to constrain the disc size, which is a key parameter in many physical phenomena: dust coagulation rate, dust-gas interaction, etc (Armitage 2010). One of the main uncertainties is the disc outer radius. We adopted two values for the outer disc: a small disc with $R_{\text{out}} = 15 \text{ AU}$ and a large disc with $R_{\text{out}} = 400 \text{ AU}$. The high value for the outer radius at 400 AU was motivated by the model of Stark et al. (2009) while a compact disc is suggested by the drop in the continuum flux at wavelengths beyond 20 micron. We also modelled various discs with a sub-standard gas-to-dust ratio (models 9 to 16).

The OI emission at 63.18 $\mu\text{m}$ probes only relatively warm and dense gases because of its upper energy level at 227 K and

<table><caption>**Table 3.** Best-fit model parameters for the fixed vertical disc structure model.</caption><tr><td>Total disc gas mass</td><td>$M_{\text{gas+dust}} (M_{\odot})$</td><td>$1.1 \times 10^{-4}$</td></tr><tr><td>Gas-to-dust mass ratio</td><td></td><td>29, 100</td></tr><tr><td>Abundance of PAHs relative to ISM</td><td>$f_{\text{PAH}}$</td><td>$10^{-3}$</td></tr><tr><td>Surface density profile index</td><td>$\epsilon$</td><td>0.4</td></tr><tr><td>Reference gas scale height</td><td>$H_0$ (AU)</td><td>0.14</td></tr><tr><td>Reference radius</td><td>$R_{\text{ref}}$ (AU)</td><td>5.0</td></tr><tr><td>Flaring index</td><td>$\gamma$</td><td>0.55</td></tr></table>

a critical density of $\sim 4.7 \times 10^5 \text{ cm}^{-3}$. In the absence of detected lines that probe cold gases, we cannot constrain the gaseous content beyond 10–20 AU well.

One of the major unknowns in the gas line modelling is the amount of PAHs in the disc. The *Spitzer*-IRS spectrum in Fig. 3 does not show the typical IR features of PAHs (Keller et al. 2008; Stark et al. 2009). We ran the models with the PAH abundance depleted with respect to the interstellar value of $3 \times 10^{-7}$ by a factor $10^{-3}$.

The second source of uncertainties are the elemental abundances relative to H-nuclei. In young protoplanetary discs where the gas comes for the collapse of a prestellar core, the elemental abundances are the same as for molecular clouds. For debris discs, the gas comes from the evaporation of icy planetesimals and the gas is enriched in elements heavier than hydrogen and helium. The evolutionary state of 51 Oph is disputed. The presence of significant amount of hot CO as evidenced by the bandhead emission, of water vapour, and of carbon dioxide support the hypothesis that the oxygen and carbon elemental abundance with respect to hydrogen are close to the molecular cloud values.

The chemical network was limited to 71 gas and ice species including formation of $H_2$ on grain surfaces and fast neutral reaction involving vibrationally excited $H_2$. Most of the reaction rates are taken from the UMIST database (Woodall et al. 2007). We adopted a standard value of $1.7 \times 10^{-17} \text{ s}^{-1}$ for the cosmic-ray ionization rate.

The code takes cooling by water lines from the near-IR to the far-IR into account. The collisional data were taken from the Leiden-Lambda database (Schöier et al. 2005). The original references for the experimental or theoretical rates are given in Appendix A. The observed [OI] flux is best reproduced by a compact (15 AU) disc with mass of $5 \times 10^{-6} M_{\odot}$. The disc gas-to-dust mass ratio is consistent with the interstellar value of 100. All other modelled fluxes are below the 3 sigma upper limits.

We ran two groups of vertical hydrostatic models: one with a compact disc and one with a large disc. The parameters of the models are listed in Table 4.

## 5. Model results

We first present the results from the Simplex automatic fitting using the fixed vertical structure disc model. The fitting procedure leads to an extremely flat disc for the gas ($H/r = 2\%$). The dust disc is even flatter because of settling (upper panels of Fig. 4). The best model has a gas mass of $1 \times 10^{-4} M_{\odot}$ and a dust mass of $1 \times 10^{-6} M_{\odot}$, similar to the values derived from the simple analysis of the millimeter emission. The best model fits to the SED and the continuum emission cumulative plots of 15% and 85% are shown in the lower panel of Fig. 4. The cumulative plots show that the continuum emission up to $\sim 100 \mu\text{m}$ comes
<table><caption>**Table 4.** Model parameters and line fluxes with 3 σ upper limits.</caption><thead><tr><th>Model</th><th><i>R</i><sub>out</sub><br>(AU)</th><th><i>M</i><sub>disc</sub><br>(<i>M</i><sub>☉</sub>)</th><th>Gas-to-dust<br>ratio</th><th>Settling</th><th>[O I] 63 µm</th><th>[O I] 145 µm<br>(10<sup>-18</sup> W m<sup>-2</sup>)</th><th>[C II] 157 µm</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Observations</td></tr><tr><td colspan="8" style="text-align:center;">53.3 ± 2.5 &lt;5.1 &lt;6.8</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Fixed vertical structure model</td></tr><tr><td>a</td><td>400</td><td>1.0 × 10<sup>-4</sup></td><td>29</td><td>yes</td><td>35.8</td><td>1.3</td><td>0.17</td></tr><tr><td>b</td><td>400</td><td>1.0 × 10<sup>-4</sup></td><td>100</td><td>yes</td><td>65.4</td><td>3.6</td><td>0.02</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Vertical hydrostatic models</td></tr><tr><td>1</td><td>15</td><td>5 × 10<sup>-6</sup></td><td>100</td><td>yes</td><td>47.2</td><td>2.8</td><td>0.02</td></tr><tr><td>2</td><td>15</td><td>5 × 10<sup>-6</sup></td><td>100</td><td>no</td><td>49.4</td><td>2.9</td><td>0.03</td></tr><tr><td>3</td><td>400</td><td>5 × 10<sup>-6</sup></td><td>100</td><td>yes</td><td>22.2</td><td>0.6</td><td>1.2</td></tr><tr><td>4</td><td>400</td><td>5 × 10<sup>-6</sup></td><td>100</td><td>no</td><td>22.2</td><td>0.6</td><td>1.2</td></tr><tr><td>5</td><td>400</td><td>5 × 10<sup>-5</sup></td><td>100</td><td>yes</td><td>65.3</td><td>2.8</td><td>1.1</td></tr><tr><td>6</td><td>400</td><td>5 × 10<sup>-5</sup></td><td>100</td><td>no</td><td>64.9</td><td>2.7</td><td>1.0</td></tr><tr><td>7</td><td>400</td><td>1 × 10<sup>-4</sup></td><td>100</td><td>yes</td><td>94.6</td><td>4.2</td><td>1.2</td></tr><tr><td>8</td><td>400</td><td>1 × 10<sup>-4</sup></td><td>100</td><td>no</td><td>93.7</td><td>4.2</td><td>1.2</td></tr><tr><td>9</td><td>15</td><td>5 × 10<sup>-7</sup></td><td>10</td><td>yes</td><td>38.2</td><td>1.1</td><td>0.03</td></tr><tr><td>10</td><td>15</td><td>5 × 10<sup>-7</sup></td><td>10</td><td>no</td><td>40.8</td><td>1.2</td><td>0.03</td></tr><tr><td>11</td><td>15</td><td>5 × 10<sup>-8</sup></td><td>1</td><td>yes</td><td>5.3</td><td>0.1</td><td>0.005</td></tr><tr><td>12</td><td>15</td><td>5 × 10<sup>-8</sup></td><td>1</td><td>no</td><td>5.3</td><td>0.1</td><td>0.005</td></tr><tr><td>13</td><td>400</td><td>5 × 10<sup>-5</sup></td><td>10</td><td>yes</td><td>144.7</td><td>4.2</td><td>9.8</td></tr><tr><td>14</td><td>400</td><td>5 × 10<sup>-5</sup></td><td>10</td><td>no</td><td>134.6</td><td>3.9</td><td>8.4</td></tr><tr><td>15</td><td>400</td><td>5 × 10<sup>-6</sup></td><td>1</td><td>yes</td><td>39.5</td><td>1.0</td><td>2.5</td></tr><tr><td>16</td><td>400</td><td>5 × 10<sup>-6</sup></td><td>1</td><td>no</td><td>39.5</td><td>1.0</td><td>2.5</td></tr></tbody></table>

from within the first 20 AU of the disc. Most of the [O I] emission is confined within 20 AU as well.

For the vertical hydrostatic discs, none of the models, either compact ($R_{out}$ = 15 AU) or extended ($R_{out}$ = 400 AU), managed to simultaneously fit the SED from near-IR to the millimeter flux density. The best simultaneous fit to the SED from the near- to the far-IR, excluding the data point at 1.2 mm and [O I] flux, is obtained for a compact disc (model 1 and 2) with a disc gas mass of $5 \times 10^{-6} M_{☉}$ and a dust mass of $5 \times 10^{-8} M_{☉}$. The fits to the SED with and without dust settling do not fit the observed continuum emission at 1.2 mm. The models underpredict the 1.2 mm flux because of the small dust mass in the models. Models 7 and 8, both with dust mass of $10^{-6} M_{☉}$ as suggested by the simple analysis of the 1.2 mm flux, overpredict all observed fluxes at wavelengths shorter than 1.2 mm. The gas scale-heights in the hydrostatic models are much higher than the value suggested by the best model with a fixed vertical structure disc. As a result, even with dust settling, the hydrostatic models of 400 AU discs produce far too strong far-IR continuum fluxes.

In the absence of spatially resolved image at 1.2 mm, we can constrain neither the location of the cold dust grains nor their geometry precisely. However, an extended flat disc has been advocated by Stark et al. (2009), who fitted the SED and mid-IR interferometry data simultaneously. An alternative explanation for the 1.2 mm emission is stellar free-free emission, which can also explain the $H\alpha$ and Bry emission.

Hydrostatic disc models with dust settling provide a better fit of the mid- and far-IR photometric points while the model without dust settling performs better in the near-IR but less well in the far-IR. This conclusion on the settling applies only for the assumption that the grains are composed of pure compact spherical amorphous silicate grains. Other grain compositions and shapes are possible. Grains with amorphous carbon inclinations emit more efficiently in the near-IR. The [O I] flux from the unsettle-model is closer to the observed value but both flux

estimates are well within the observational uncertainties. We tested the effect of well-mixed gas and dust (no settling) model and all other parameters the same. The fit to the SED of the well-mixed model is slightly worse (overpredicting the flux in the far-IR) than in the model with settling and the [O I] emissions are comparable, albeit different (Table 4).

The disc vertical hydrostatic structure shown in Fig. 5 is complex but clearly drops beyond 8 AU due a drop in the gas temperature in both the settled and well-mixed models. The dust density structure follows that of the gas structure (Fig. 6). The dust grains are warmer farther out in the settled-dust model, which is expected since the grains are on average smaller in the disc upper layers. We also show the gas temperature in Fig. 9, which is much higher than the dust temperature in the disc atmosphere.

The [O I] flux from models 1 and 2 slightly unpredicts the observed flux. The [O I] flux emission area is shown in Fig. 10. The [O I] emission at 63 micron arises from the disc midplane within the first 10–15 AU. The [O I] line does probably not probe the same circumstellar region as the 1.2 mm continuum emission. It is clear that the [O I] line cannot be used to probe cold gas emission.

In all models, atomic oxygen is by far the main oxygen carrier in the disc, followed by CO. The disc midplane is molecular and CO reaches its maximum abundance of $1.3 \times 10^{-4}$. In our models the [O I] emission traces the disc molecular gas. The disc around 51 Oph is mostly molecular ($H_2$, CO, etc.) and poor in $C^+$ apart from the upper disc atmospheres. The resulting C II fluxes are low, consistent with the observed upper limit. There is an intermediate layer where neutral atomic carbon is abundant. The oxygen 63 micron line is emitted within 15 AU of the star regardless of the physical size of the disc.

The main gas coolant in the [O I] emitting region are the CO ro-vibrational lines. The [O I] line probes a tiny fraction of the total disc gas mass only. In the absence of PAHs the main heating
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Disc gas density, gas-to-dust ratio, SED, and localization of the continuum emission for the fixed vertical structure disc model. The disc is gas-rich above the mid-plane and gas-poor in the mid-plane. For clarity a couple of observed photometric points are not shown in the SED. The lower-right panel shows the cumulative continuum emission for 15% and 85% for each wavelength. The continuum emission for wavelengths greater than 500 µm comes from the outer disc ($R > 100$ AU).</figcaption></figure>

is the CO IR-pumping and background heating by atoms and ions (FeII, SiII, CII, CI) in the first 5 AU and dust photoelectric on dust grains beyond 5 AU. Water is a very minor component, and its ro-vibrational lines do not contribute to the gas cooling.

A larger disc ($R_{\text{out}} = 400$ AU) with the same gas mass as the best compact model predicts the continuum and [O I] line flux. The modelling confirms that we cannot constrain the actual disc outer radius. However, the location of the O I emission is confined to the first 10–15 AU for both the compact and the more extended disc. The more massive large discs (models 7 and 8) do not match the SED and overpredict the [O I] flux. A few models with sub-standard gas-to-dust mass ratio have [O I] flux close to the observed values. However, like all other hydrostatic models, the fits to the SED are poor.

**6. Discussion**

What is the origin of the gas and dust grains in the inner disc? Do the grains come from the coagulation of interstellar grains,

or are they the results of collisions between planetesimals scattered from the outer belt reservoir by unseen companions? Is the gas the remnant of an initial massive disc? 51 Oph is young: ~0.7 Myr (Montesinos et al. 2009), while the youngest stars with established debris discs are about 10 Myr old (β Pic has an age of 12 Myr). Could the disc around 51 Oph be considered as a debris disc at its age?

The best fixed vertical structure model is a flat gas disc whereas the gas scale height derived from the hydrostatic models is much more extended vertically. The comparison between the fixed vertical structure disc model and the hydrostatic models suggests that there may be a mechanism to flatten the gas disc around 51 Oph. The vertical hydrostatic structure is derived from balancing the thermal pressure and the gravity forces.

The disc average gas-to-dust mass ratio is ~30–100 for the best model, consistent with the standard interstellar value and with the ratio in other Herbig Ae discs: HD 169142 (Meeus et al. 2010) and HD 163296 (Tilling et al. 2012). A gas-to-dust mass
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Disc gas density structure from hydrostatic models 1 and 2 ($R_{\text{out}} = 2$ AU, $M_{\text{disc}} = 5 \times 10^{-6} M_{\odot}$) between 2 and 15 AU. The left panel corresponds to the model with dust settling (model 1), the right panel to the model without dust settling (model 2). The dust temperature contours are overplotted in black. The red contours correspond to the Stokes numbers ($S_t = 0.1, 1$, and 10 contour values).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Disc dust density structure from hydrostatic models 1 and 2 ($R_{\text{out}} = 2$ AU, $M_{\text{disc}} = 5 \times 10^{-6} M_{\odot}$) between 2 and 15 AU. The left panel corresponds to the model with dust settling (model 1), the right panel to the model without dust settling (model 2). The dust temperature contours are overplotted in black. The Stokes parameters are overplotted in red ($S_t = 0.1, 1$, and 10 contour values).</figcaption></figure>

ratio of 100 is much higher than the ratio in the debris discs (e.g., Lebreton et al. 2012), where the dust grains come from shattering of planetesimals and the gas is absent. Nevertheless, the 51 Oph disc gas mass ($5 \times 10^{-6} M_{\odot}$) is much lower than the gas mass of the aforementioned Herbig Ae discs, which can reach $0.01 M_{\odot}$. If the disc around 51 Oph is primordial, as suggested by its young age, then either the initial disc had a low mass or the initial massive disc has been actively dissipating either by photoevaporation, accretion onto the star, and/or planet formation. As discussed above, the current mass accretion rate onto 51 Oph is difficult to assess, not to mention the accretion rate history. In discs where the gas is replenished by planetesimal evaporation, the gas-to-dust ratio would be closer to unity. Hydrodynamic simulations show that dust grains can gather into

eccentric rings as observed in near-IR images of debris discs without the need of the presence of a planet (Lyra & Kuchner 2012).

The high luminosity of the star ($L_* = 260–300 L_{\odot}$) compared with other Herbig Ae stars ($L_* = 20–40 L_{\odot}$) may explain a fast disc photoevaporation or increases in the efficiency of radiative pressure on grains resulting in a radially extended disc. Disc photo-evaporation occurs when hot gases have enough energy to escape the star+disc gravitational potential entraining the dust grains that are dynamically coupled to the gas. The upper disc atmosphere outside the critical radius of ~6.8 AU for 51 Oph is unbound and flows as a wind (Hollenbach et al. 1994). The mass-loss rate through disc evaporative wind is $\dot{M}_w \approx 4 \times 10^{-10} M_{\odot} \text{yr}^{-1}$ assuming a stellar
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Grain mass absorption ($\kappa_{abs}$ in red dotted-line), scattering ($\kappa_{sca}$ in black dashed-line), and extinction ($\kappa$ in black full line) opacity for the adopted dust grain properties listed in Table 2. The slope in the upper-right corner shows a hypothetical opacity with $\kappa \sim \lambda^{-1}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Fit to the SED for model 1 with and without dust settling. The four far-IR photometric points are Herschel-PACS observations. The 1.2 mm photometry was obtained with MAMBO at the IRAM 30 m telescope.</figcaption></figure>

ionizing flux of $\Phi = 10^{42}$ ionizing photons per second, valid for a low-mass solar-luminosity star (Alexander & Armitage 2009; Pascucci et al. 2011). Scaling to the luminosity of 51 Oph, the wind mass-loss rate reaches $1.2 \times 10^{-7}$ yr$^{-1}$. A detailed disc-wind model is needed to test whether the mass-loss rate of $\sim 10^{-7}$ yr$^{-1}$ can explain the observed H$\alpha$ and Bry line fluxes and profiles. In 0.7 Myr, the disc around 51 Oph would have lost $0.1 M_{\odot}$ of gas. If photoevaporation is occurring at a constant rate, the initial disc gas mass around 51 Oph would have been $\sim 0.1 M_{\odot}$ compared with the current stellar mass of $2.8 M_{\odot}$. Since the inner disc is warm, the initial disc would have been gravitationally stable. The current disc within 10–15 AU will quickly disappear at this high photoevaporation rate. A mass-loss rate of $\sim 10^{-7}$ yr$^{-1}$ is very high compared with the mass in the disc ($\sim 5 \times 10^{-6} M_{\odot}$) and may be overestimated. A detailed disc-wind model is warranted to better constrain the photoevaporative mass-loss rate.

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Gas temperature structure for model 1. The gas temperature is higher than 250 K in the inner compact disc.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** [O I] emission from the disc for model 1. The upper panel shows the continuum at 63 microns and the line vertical optical depth. The middle panel shows the cumulative [OI] 63 micron emission from the inner to the outer disc. The lower panel exhibits the atomic oxygen density distribution in the disc. The contour encloses the 70% × 70% emission area.</figcaption></figure>

The dust mass in the small compact disc around 51 Oph derived from the simultaneous fit to the SED (near- to far-IR photometry data) using the hydrostatic disc model is small at $5 \times 10^{-8} M_{\odot}$ (0.017 $M_{\oplus}$), two orders of magnitude lower than if all the 1.2 mm flux were attributed to dust emission. However, this dust mass is still ten times the value derived by van den Ancker et al. (2001) Their value was derived from a fit to photometric points at wavelengths shorter than 100 microns only. Our derived dust mass is five times higher than the mass found by Tatulli et al. (2008).
The discrepancies can also be attributed to differences in dust opacities. The continuum emission at 63 micron is optically thin, thus it probes the amount of cool dust in the disc.

The best models with a fixed vertical structure or with a disc in vertical hydrostatic equilibrium suggest that the gas-to-dust mass ratio is ~30-100, but are the gas and dust well-mixed or have the dust grains settled to the midplane? We used the model results to compute the Stokes parameter (Weidenschilling 1977),

$$S_t = \frac{\Omega_K}{2\pi} \frac{a_d m_d}{c_s \rho_g}, \quad (4)$$

where $\Omega_K = (GM_*/R^3)^{1/2}$ is the Keplerian angular velocity at radius $R$, $a_d$ the grain radius in cm, $m_d$ is the volume mass density of a dust grain, which is ~1 g cm⁻³ for a silicate grain, $c_s$ the sound speed in cm s⁻¹, and $\rho_g = n_H \times a$. a.u. the gas mass density in g cm⁻³ (a.u. is the atomic mass unit in grams). The Stokes parameter varies throughout the disc. It measures the level of dust grain coupling with the gas via drag force. For $S_t \gg 1$, the gas drag is inefficient and the dust grains quickly settle towards the disc midplane. For $S_t \ll 1$, the gas and dust grains are well mixed and are collocated. The contours for $\log_{10}$ the Stokes parameters (values of -10, 1, and 10) in the disc models with and without settling for the mean grain size are plotted over the density distribution in Fig. 5 and over the dust density distribution in Fig. 6. The Stokes parameter is small throughout the disc in our models because the gas is dense ($n_H > 10^6$ cm⁻³) and warm ($T_{gas} > 250$ K).

Fits to the SED and gas lines do not strongly constrain the presence of dust settling. However, the SED together with VLT-MIDI and Keck nulling data support models with dust settling (Stark et al. 2009). A more detailed modelling is needed to compute the exact dust settling, coagulation, and migration, which is beyond the scope of this paper (e.g. Barrière-Fouchet et al. 2005).

Interestingly, dust settling does not much affect the gas vertical hydrostatic structure for two reasons First, photoelectric-heating by dust grains is not the only heating process in the disc layers located above the midplane. Heating by PAHs, which are small enough to be coupled to the gas, is important, but their abundance is very low. Second, sub-micron grains, which are more strongly coupled to the gas and remain in the disc atmosphere, provide most of the photoelectrons. Likewise, the gas line collisional excitations do not change much in our model of the 51 Oph disc whether the dust grains have settled or not.

If we can attribute the millimeter emission to a cold outer disc instead of free-free emission and the near- to mid-IR continuum emission and the [O I] to a warm compact inner disc, the circumstellar environment of 51 Oph resembles a young and massive counterpart to main-sequence stars with a warm exozodiacal dust ring and a cold outer planetesimal belt (Defrère et al. 2012). The main difference between 51 Oph and β Pic is that the dust grains in an exozodiacal ring around β Pic have to be replenished, while part the dust grains around 51 Oph may still be primordial.

## 7. Conclusions

We detected [O I] far-IR and 1.2 mm continuum emission towards 51 Oph. The [O I] probes the warm gas in a compact inner circumstellar disc while the 1.2 mm emission suggests a cold outer disc. An extended flat disc with dust settling can explain

the continuum and line data consistent with previous model. Hydrostatic disc models predict discs that are too extended in the vertical direction for the gas and the dust. The failure of the hydrostatic models to explain the observations suggests that mechanism(s) other than thermal pressure and gravity forces are at work to keep the gas and dust disc flat. The gas-to-dust mass ratio is consistent with the standard value of 100. Future observations in the (sub)-millimeter with ALMA will help in constraining the cold gas and dust around 51 Oph.

*Acknowledgements.* We thank the Agence Nationale pour la Recherche (contracts ANR-07-BLAN-0221, ANR-2010-JCJC-0504-01, and ANR-2010-JCJC-0501-01), CNES, PNPS of CNRS/INSU, France, and support from the Millenium Science Initiative (Chilean Ministry of Economy), through grant "Nucleus P10-022-F" for support. We acknowledge funding from the EU FP7-2011 under Grant Agreement nr. 284405 (PERG06-GA-2009-256513). Computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI) on the super-computer Fostino funded by ANR. C. Eiroa, G. Meeus, and B. Montesinos are partly supported by Spanish grant AYA 2011-26202. We finally thank the referee for his/her comments that helped improving the manuscript.

## Appendix A: Collisional data

The original articles for the line frequencies, Einstein coefficients, and collisional rates are CO (Flower 2001; Jankowski & Szalewicz 2005; Yang et al. 2006; Wernli et al. 2006), H₂O (Barber et al. 2006; Dubernet & Grosjean 2002; Faure et al. 2004, 2007; Daniel et al. 2011), O I (Abrahamsson et al. 2007; Bell et al. 1998; Chambaud et al. 1980; Jaquet et al. 1992; Launay & Roueff 1977), C II (Flower & Launay 1977; Launay & Roueff 1977; Wilson & Bell 2002), CH⁺ (Müller 2010; Lim et al. 1999; Hammami et al. 2009; Turpin et al. 2010), OH (Offer et al. 1994).

## References

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

Alexander, R. D., & Armitage, P. J. 2009, ApJ, 704, 989

Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge University Press)

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

Barrière-Fouchet, L., Gonzalez, J.-F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250

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

Berthoud, M. G., Keller, L. D., Herter, T. L., Richter, M. J., & Whelan, D. G. 2007, ApJ, 660, 461

Brittain, S. D., Simon, T., Najita, J. R., & Rettig, T. W. 2007, ApJ, 659, 685

Chambaud, G., Levy, B., Millie, P., et al. 1980, J. Phys. B Atom. Mol. Phys., 13, 4205

Daniel, F., Dubernet, M.-L., & Grosjean, A. 2011, A&A, 536, A76

Defrère, D., Lebreton, J., Le Bouquin, J.-B., et al. 2012, A&A, 546, L9

Dent, W. R. F., Thi, W. F., Kamp, I., et al. 2013, PASP, 125, 477

Draine, B. T. 2006, ApJ, 636, 1114

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

Fajardo-Acosta, S. B., Telesco, C. M., & Knacke, R. F. 1993, ApJ, 417, L33

Faure, A., Gorfinkiel, J. D., & Tennyson, J. 2004, MNRAS, 347, 323

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

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

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

Garcia Lopez, R., Natta, A., Testi, L., & Habart, E. 2006, A&A, 459, 837

Grady, C. A., & Silvis, J. M. S. 1993, ApJ, 402, L61

Hammami, K., Owono Owono, L. C., & Stäuber, P. 2009, A&A, 507, 1083

Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654

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

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

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

Keller, L. D., Sloan, G. C., Forrest, W. J., et al. 2008, ApJ, 684, 411
Kraus, S., Hofmann, K.-H., Benisty, M., et al. 2008, A&A, 489, 1157

Kreysa, E., Gemuend, H.-P., Gromke, J., et al. 1998, in SPIE Conf. Ser. 3357, ed. T. G. Phillips, 319

Laor, A., & Draine, B. T. 1993, ApJ, 402, 441

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

Lebreton, J., Augereau, J.-C., Thi, W.-F., et al. 2012, A&A, 539, A17

Lecavelier Des Etangs, A., Vidal-Madjar, A., Backman, D. E., et al. 1997, A&A, 321, L39

Lim, A. J., Rabadán, I., & Tennyson, J. 1999, MNRAS, 306, 473

Lyra, W., & Kuchner, M. J. 2012 [arXiv:1204.6322]

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

Meeus, G., Montesinos, B., Mendigutía, I., et al. 2012, A&A, 544, A78

Mendigutía, I., Calvet, N., Montesinos, B., et al. 2011, A&A, 535, A99

Mendigutía, I., Mora, A., Montesinos, B., et al. 2012, A&A, 543, A59

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

Müller, H. S. P. 2010, A&A, 514, L6

Offer, A. R., van Hemert, M. C., & van Dishoeck, E. F. 1994, J. Chem. Phys., 100, 362

Ott, S. 2010, in Astronomical Data Analysis Software and Systems XIX, eds. Y. Mizumoto, K.-I. Morita, & M. Ohishi, ASP Conf. Ser., 434, 139

Pascucci, I., Sterzik, M., Alexander, R. D., et al. 2011, ApJ, 736, 13

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

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

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Roberge, A., Feldman, P. D., Lecavelier des Etangs, A., et al. 2002, ApJ, 568, 343

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

Stark, C. C., Kuchner, M. J., Traub, W. A., et al. 2009, ApJ, 703, 1188

Sylvester, R. J., Skinner, C. J., Barlow, M. J., & Mannings, V. 1996, MNRAS, 279, 915

Tatulli, E., Malbet, F., Ménard, F., et al. 2008, A&A, 489, 1151

Thi, W.-F., van Dalen, B., Bik, A., & Waters, L. B. F. M. 2005, A&A, 430, L61

Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711

Tilling, I., Woitke, P., Meeus, G., et al. 2012, A&A, 538, A20

Turpin, F., Stoecklin, T., & Voronin, A. 2010, A&A, 511, A28

van den Ancker, M. E., Meeus, G., Cami, J., Waters, L. B. F. M., & Waelkens, C. 2001, A&A, 369, L17

Weidenschilling, S. J. 1977, MNRAS, 180, 57

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

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

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

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

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2006, J. Chem. Phys., 124, 104304