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

HAL
open science

# A gas density drop in the inner 6 AU of the transition disk around the Herbig Ae star HD 139614

A. Carmona, W. Thi, I. Kamp, C. Baruteau, A. Matter, M. van den Ancker, C. Pinte, Á. Kóspál, M. Audard, A. Liebhart, et al.

► To cite this version:

A. Carmona, W. Thi, I. Kamp, C. Baruteau, A. Matter, et al. A gas density drop in the inner 6 AU of the transition disk around the Herbig Ae star HD 139614. Astronomy & Astrophysics - A&A, 2017, 598, pp.A118. 10.1051/0004-6361/201628472. hal-02369665

HAL Id: hal-02369665

https://hal.science/hal-02369665v1

Submitted on 27 Jun 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A gas density drop in the inner 6 AU of the transition disk around the Herbig Ae star HD 139614

## Further evidence for a giant planet inside the disk?★

A. Carmona<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, W. F. Thi<fnref n="4" />, I. Kamp<fnref n="5" />, C. Baruteau<fnref n="1" />, A. Matter<fnref n="6" />, M. van den Ancker<fnref n="7" />, C. Pinte<fnref n="8" />,<fnref n="9" />, A. Kóspál<fnref n="2" />, M. Audard<fnref n="10" />, A. Liebhart<fnref n="11" />, A. Sicilia-Aguilar<fnref n="12" />, P. Pinilla<fnref n="13" />, Zs. Regály<fnref n="2" />, M. Güdel<fnref n="11" />, Th. Henning<fnref n="14" />, L. A. Cieza<fnref n="15" />, C. Baldovin-Saavedra<fnref n="11" />, G. Meeus<fnref n="3" />, and C. Eiroa<fnref n="3" />,<fnref n="16" />

Received 9 March 2016 / Accepted 15 September 2016

## ABSTRACT

**Context.** Quantifying the gas surface density inside the dust cavities and gaps of transition disks is important to establish their origin. **Aims.** We seek to constrain the surface density of warm gas in the inner disk of HD 139614, an accreting 9 Myr Herbig Ae star with a (pre-)transition disk exhibiting a dust gap from 2.3 ± 0.1 to 5.3 ± 0.3 AU.

**Methods.** We observed HD 139614 with ESO/VLT CRIRES and obtained high-resolution (R ~ 90 000) spectra of CO ro-vibrational emission at 4.7 µm. We derived constraints on the disk’s structure by modeling the CO isotopolog line-profiles, the spectroastrometric signal, and the rotational diagrams using grids of flat Keplerian disk models.

**Results.** We detected $v = 1 \rightarrow 0$ $^{12}$CO, $2 \rightarrow 1$ $^{12}$CO, $1 \rightarrow 0$ $^{13}$CO, $1 \rightarrow 0$ C$^{18}$O, and $1 \rightarrow 0$ C$^{17}$O ro-vibrational lines. Lines are consistent with disk emission and thermal excitation. $^{12}$CO $v = 1 \rightarrow 0$ lines have an average width of 14 km s$^{-1}$, $T_{gas}$ of 450 K and an emitting region from 1 to 15 AU. $^{13}$CO and C$^{18}$O lines are on average 70 and 100 K colder, 1 and 4 km s$^{-1}$ narrower than $^{12}$CO $v = 1 \rightarrow 0$, and are dominated by emission at R ≥ 6 AU. The $^{12}$CO $v = 1 \rightarrow 0$ composite line-profile indicates that if there is a gap devoid of gas it must have a width narrower than 2 AU. We find that a drop in the gas surface density ($\delta_{gas}$) at R < 5–6 AU is required to be able to simultaneously reproduce the line-profiles and rotational diagrams of the three CO isotopologs. Models without a gas density drop generate $^{13}$CO and C$^{18}$O emission lines that are too broad and warm. The value of $\delta_{gas}$ can range from $10^{-2}$ to $10^{-4}$ depending on the gas-to-dust ratio of the outer disk. We find that the gas surface density profile at $1 < R < 6$ AU is flat or increases with radius. We derive a gas column density at $1 < R < 6$ AU of $N_H = 3 \times 10^{19}–10^{21}$ cm$^{-2}$ ($7 \times 10^{-5}–2.4 \times 10^{-3}$ g cm$^{-2}$) assuming $N_{CO} = 10^{-4}N_H$. We find a $5\sigma$ upper limit on the CO column density $N_{CO}$ at R ≤ 1 AU of $5 \times 10^{15}$ cm$^{-2}$ ($N_H \le 5 \times 10^{19}$ cm$^{-2}$).

**Conclusions.** The dust gap in the disk of HD 139614 has molecular gas. The distribution and amount of gas at R ≤ 6 AU in HD 139614 is very different from that of a primordial disk. The gas surface density in the disk at R ≤ 1 AU and at $1 < R < 6$ AU is significantly lower than the surface density that would be expected from the accretion rate of HD 139614 ($10^{-8} M_{\odot}$ yr$^{-1}$) assuming a standard viscous $\alpha$-disk model. The gas density drop, the non-negative density gradient in the gas inside 6 AU, and the absence of a wide (>2 AU) gas gap, suggest the presence of an embedded <2 $M_J$ planet at around 4 AU.

**Key words.** protoplanetary disks – stars: pre-main sequence – planets and satellites: formation – techniques: spectroscopic – stars: variables: T Tauri, Herbig Ae/Be

<footnote n="1">¹ Université de Toulouse, UPS-OMP, IRAP, 14 avenue E. Belin, Toulouse, 31400, France
e-mail: andres.carmona@irap.omp.eu</footnote>

<footnote n="2">² Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, PO Box 67, 1525 Budapest, Hungary</footnote>

<footnote n="3">³ Departamento de Física Teórica, Universidad Autónoma de Madrid, Campus Cantoblanco, 28049 Madrid, Spain</footnote>

<footnote n="4">⁴ Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching bei München, Germany</footnote>

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

<footnote n="6">⁶ Laboratoire Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Boulevard de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France.</footnote>

<footnote n="7">⁷ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany</footnote>

<footnote n="8">⁸ UMI-FCA, CNRS/INSU, France (UMI 3386), and Dept. de Astronomía, Universidad de Chile, 1058 Santiago, Chile</footnote>

<footnote n="9">⁹ Univ. Grenoble Alpes, IPAG; CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="10">¹⁰ Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland</footnote>

<footnote n="11">¹¹ University of Vienna, Department of Astronomy, Türkenschanzstrasse 17, 1180 Vienna, Austria</footnote>

<footnote n="12">¹² SUPA, School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK</footnote>

<footnote n="13">¹³ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="14">¹⁴ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="15">¹⁵ Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejercito 441, 1058 Santiago, Chile</footnote>

<footnote n="16">¹⁶ Unidad Asociada Astro-UAM/CSIC, 28850 Madrid, Spain</footnote>

<footnote>* Based on CRIRES observations collected at the VLTI and VLT (European Southern Observatory, Paranal, Chile) with program 091.C-0671(B).</footnote>

<footnote>** Part of this research has been done by A. Carmona under the frame of ESO’s scientist visitor program during November 2013 and at Université Grenoble Alpes, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France. Current address: Institut de Recherche en Astrophysique et Planétologie (IRAP), 14 avenue E. Belin, Toulouse, 31400, France.</footnote>
# 1. Introduction

Transition disks are protoplanetary disks that exhibit a deficit of continuum emission at near- and/or mid-IR wavelengths in their spectral energy distribution (for a recent review, see Espaillat et al. 2014). This deficit of emission is commonly interpreted as evidence of a dust gap, a dust cavity, or a dust hole inside the disk<fnref n="1" />. Sub-mm interferometry observations have confirmed the existence of dust cavities by spatially resolving the thermal emission from cold large (~mm) grains at tens of AU in transition disks (e.g., Piétu et al. 2006; Brown et al. 2009; Andrews et al. 2011; Cieza et al. 2012; Casassus et al. 2013; Pérez et al. 2014). Observations of scattered light in the near-IR using adaptive optics have further confirmed the dust cavities in micron-sized dust grains. These high-spatial resolution observations show that the cavity size in small grains can be smaller than that in large grains (e.g., Muto et al. 2012; Garufi et al. 2013; Follette et al. 2013; Pinilla et al. 2015a). Furthermore, near-IR scattered light imaging and sub-mm interferometry observations have revealed that a large fraction of transition disks has asymmetries in the dust distribution (e.g. spirals, blobs, and horsehoe shapes), although, the presence and shape of asymmetries appear to be different depending on the wavelength of the observations and thus the dust sizes traced (e.g., Muto et al. 2012; van der Marel et al. 2013; Isella et al. 2013; Pérez et al. 2014; Benisty et al. 2015; Follette et al. 2015).

The origin of the dust cavities and gaps in transition disks is a matter of intense debate in the literature: scenarios such as grain growth (e.g., Dullemond & Dominik 2005; but see Birnstiel et al. 2012), size-dependent dust radial drift (e.g., Pinte & Laibe 2014), dust dynamics at the boundary of the dead-zone (Regály et al. 2012), photoevaporation (e.g., Clarke et al. 2001; Alexander & Armitage 2007; Owen et al. 2012), giant planet(s) (e.g., Marsh & Mahoney 1992; Lubow et al. 1999; Rice et al. 2003; Quillen et al. 2004; Varnière et al. 2006; Zhu et al. 2011), dynamical interactions in multiple systems (e.g., Artymowicz & Lubow 1996; Ireland & Kraus 2008; Fang et al. 2014), and magneto-hydrodynamical phenomena (Chiang & Murray-Clay 2007) have all been proposed.

Accretion signatures in many transition disks (e.g., Fang et al. 2009; Sicilia-Aguilar et al. 2013; Manara et al. 2014) and emission of warm (e.g., Bary et al. 2003, Pontoppidan et al. 2008, 2011, Salyk et al. 2009, 2011) and cold (Casassus et al. 2013; Bruderer et al. 2014; Perez et al. 2015; Canovas et al. 2015; van der Marel et al. 2015b, 2016) molecular gas indicate that the dust cavities in accreting transition disks contain gas. Radiative transfer modeling of CO ro-vibrational emission (Carmona et al. 2014) and CO pure rotational emission (Bruderer 2013; Perez et al. 2015; van der Marel et al. 2015b, 2016) further suggests a gas surface density drop ($\delta_{gas}$) inside the dust cavity, with $\delta_{gas}$ values varying from 0.1 up to $10^{-5}$ (see Table C.1). Some of the transition disks are not accreting and thus do not seem to have gas (Sicilia-Aguilar et al. 2010). There is also a substantial difference in the global structure and/or disk mass between accreting and non-accreting transition disks, with the non-accreting disks being significantly more evolved (lower

masses, flatter disks) as seen with Herschel (Sicilia-Aguilar et al. 2015).

The different spatial locations of dust grains of different sizes, the gas inside the sub-mm dust cavities, together with the different surface density profiles of gas and dust strongly favor the planet(s) scenario. However, we probably witness several coexisting mechanisms, because planet formation might affect the dynamics of the dust in the disk (e.g., Rice et al. 2003; Zhu et al. 2011; Pinilla et al. 2012, 2015b) or favor the onset of photoevaporation, when the accretion rate has decreased (e.g., Rosotti et al. 2013; Dittkrist et al. 2014). A large portion of studies of transition disks have focused on investigating disks that are bright in the sub-mm and that have large dust cavities of tens of AU (e.g., Andrews et al. 2011; van der Marel et al. 2015a). Because a single Jovian planet interacting with the disk is expected to open a gap only a few AU wide (e.g., Kley 1999; Crida & Morbidelli 2007), multiple (unseen) giant planets have been postulated as a possible explanation for the observed large dust cavities (Zhu et al. 2011; Dodson-Robinson & Salyk 2011).

In a recent near- and mid-IR interferometry campaign, Matter et al. (2014, 2016) have revealed that the 9 Myr old (Alecian et al. 2013) accreting ($10^{-8} M_{\odot}$/yr, Garcia Lopez et al. 2006) Herbig A7Ve star HD 139614 has a transition disk with a narrow dust gap extending from 2.3 ± 0.1 to 5.3 ± 0.3 AU<fnref n="2" />, and a dust density drop $\delta_{dust}$ at $R < 6$ AU of $10^{-4}$ (see Table 1 for a summary of the stellar properties). HD 139614 is one of the first objects with a spatially resolved dust gap with a width of only a few AU, thus it might be the case of a transition disk where the dust gap has been opened by a single giant planet.

HD 139614 is located within the Sco OB2-3 association (Acke et al. 2005) at a distance of 131±5 pc (Gaia Collaboration 2016). HD 139614 has peculiar chemical abundances in its photosphere (Folsom et al. 2012), with depletions of heavier refractory elements, while C, N, and O are approximately solar. HD 139614 belongs to the group I Herbig Ae stars according to the Spectral Energy Distribution (SED) classification scheme of Meeus et al. (2001), which suggests that its outer disk is flared. Matter et al. (2016) derived a dust disk mass of $10^{-4} M_{\odot}$ based on a fit to the SED. The Spitzer mid-IR spectra of HD 139614 exhibit a weak amorphous silicate feature at 10 µm (Juhász et al. 2010) and Polycyclic Aromatic Hydrocarbons (PAH) emission (Acke et al. 2010). The disk's mid-IR continuum has been spatially resolved at 18 µm (FWHM of 17 ± 4 AU) but it is not resolved at 12 µm (Mariñas et al. 2011). Kóspál et al. (2012) reported that the ISOPHOT-S, Spitzer and TIMMI-2/ESO 3.6m mid-infrared spectra taken at different epochs agree within the measurement uncertainties, thus suggesting that there is no strong mid-IR variability in the source. Emission from cold CO gas in the outer disk of HD 139614 has been reported in JCMT single-dish observations by Dent et al. (2005) and Panić & Hogerheijde (2009). Emission of [O I] at 63 µm from the disk has been detected by Herschel (Meeus et al. 2012; Fedele et al. 2013). The [O I] 63 µm line flux of HD 139614 is among the weakest of the whole Herbig Ae sample observed by Herschel. No emission of [O I] at 145 µm, [C II] at 157 µm, CO, H₂O, OH or CH⁺ in the 50–200 µm region was detected by Herschel (Meeus et al. 2012, 2013; Fedele et al. 2013).

In this paper we present the results of high-resolution spectroscopy observations of CO ro-vibrational emission at 4.7 µm

<footnote n="1">¹ We call a dust hole, when no dust emission is detected inside a determined radius in the disk at *all* wavelengths. We call a dust cavity, a region where there is a drop in the dust density. Inside the dust cavity radius dust is still present (i.e. continuum emission is detected inside the cavity radius, for instance at IR wavelengths). We call a dust gap when continuum emission is detected at radii smaller and larger than the location of the gap. A dust cavity can have a dust gap inside it.</footnote>

<footnote n="2">² The dust gap limits derived in Matter et al. (2016) are 2.5 ± 0.1 to 5.7 ± 0.3 AU. They were calculated using a distance of 140 pc. The values in the text are the values corrected by the new Gaia distance. Both values are consistent within the uncertainties.</footnote>
<table><caption>**Table 1.** Stellar properties.</caption><thead><tr><th>Star</th><th>Sp. type</th><th><i>T</i><sub>eff</sub><br>[K]</th><th><i>d</i><br>[pc]</th><th>Mass<br>[M<sub>☉</sub>]</th><th>Radius<br>[R<sub>☉</sub>]</th><th>RV<br>[km s<sup>-1</sup>]</th><th><i>W</i><sub>2</sub><br>[mag]</th><th>Age<br>[Myr]</th><th><i>i</i><sub>disk</sub><br>[°]</th><th><i>L</i><sub>x</sub><br>[erg s<sup>-1</sup>]</th><th><i>M</i><sup>*</sup><br>[M<sub>☉</sub> yr<sup>-1</sup>]</th></tr></thead><tbody><tr><td rowspan="2">HD 139614</td><td rowspan="2">A7Ve<sup>a</sup></td><td>7600 ± 300<sup>b</sup></td><td>131 ± 5<sup>c</sup></td><td>1.76<sup>+0.15b</sup><sub>-0.08</sub></td><td>2.06 ± 0.42<sup>b</sup></td><td>0.3 ± 2.3<sup>b</sup></td><td>5.1<sup>d</sup></td><td>8.8<sup>+4.5b</sup><sub>-1.9</sub></td><td>20<sup>e</sup></td><td>1.2 × 10<sup>29f</sup></td><td>10<sup>-8g</sup></td></tr><tr><td>7850<sup>a</sup></td><td></td><td>1.7 ± 0.3<sup>h</sup></td><td>1.6<sup>h</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**References.** (a) van Boekel et al. (2005); (b) Folsom et al. (2012), Alecian et al. (2013); (c) Gaia Collaboration (2016); (d) 4.6 µm, WISE satellite release 2012 (Cutri 2012); (e) Matter et al. (2016); (f) Güdel et al. (in prep.) see Sect. 5.4; (g) Garcia Lopez et al. (2006); (h) van Boekel et al. (2005), stellar properties used in Matter et al. (2016).

<table><caption>**Table 2.** Log of the science and calibrator observations.</caption><thead><tr><th>Star</th><th>UT Date Obs.<br>[y-m-d]</th><th><i>t</i><sub>exp</sub><br>[s]</th><th>Airmass</th><th>Seeing<br>["]</th><th>RV<sub>bary</sub><sup>a</sup><br>[km s<sup>-1</sup>]</th><th>PSF<sub>FWHM</sub><sup>b</sup><br>[mas]</th><th><i>S</i>/<i>N</i><sup>b,c</sup></th><th>Sensitivity 3σ<sup>b,d</sup><br>[10<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup>]<br>3.3 km s<sup>-1</sup> 20 km s<sup>-1</sup></th></tr></thead><tbody><tr><td>HD 139614</td><td>2013-06-15</td><td>2400</td><td>1.07–1.13</td><td>0.93–1.23</td><td>9.74 ± 0.02</td><td>178 ± 10</td><td>160–100</td><td>0.2–0.3</td></tr><tr><td>CAL HIP 76829</td><td>2013-06-15</td><td>320</td><td>1.16–1.18</td><td>0.87–1.06</td><td>9.36 ± 0.01</td><td>172 ± 10</td><td>310–200</td><td>1.2–2.0</td></tr></tbody></table>

**Notes.** (a) Radial velocity due to the rotation of the Earth, the motion of the Earth about the Earth-Moon barycenter, and the motion of the Earth around the Sun; (b) measured in one nod position; (c) for the science spectra the S/N is measured in the telluric-corrected spectrum, note that the S/N decreases from chip 1 to chip 4; (d) integrated flux sensitivity limits are given for a spectrally unresolved line of width 3.3 km s⁻¹ and a line of width 20 km s⁻¹.

towards HD 139614 obtained with the ESO/VLT CRIRES instrument (Käufli et al. 2014). Our aim is to use CO isotopolog spectra to constrain the warm gas content in the inner disk of HD 139614 and address the following questions: What is the gas distribution in the inner disk of HD 139614? Does the HD 139614 disk have a gas-hole, a gas-density drop or a gap in the gas? How does the gas distribution compare with the dust distribution? What is the most likely explanation for the observed gas and dust distributions in HD 139614?

The paper is organized as follows. We start by describing the observations and data reduction in Sect. 2. In Sect. 3, we present the observational results. In Sect. 4, we derive the CO-emitting region, the average temperature and column density of the emitting gas, and the gas surface density and temperature distribution. In Sect. 5, we discuss our results in the context of the proposed scenarios for the origin of transition disks and compare HD 139614 with other transition disks. Section 6 summarizes our work and provides our conclusions.

## 2. Observations and data reduction

HD 139614 was observed with the high-resolution near-IR spectrograph CRIRES at the ESO Very Large Telescope at Cerro Paranal Chile in June 2013. CRIRES has a pixel scale of 0.086 arcsec/pixel in the spatial direction (11 AU at 131 pc) and 2.246 × 10⁻⁶ µm in the wavelength direction (0.14 km s⁻¹ at 4.7 µm). Observations were performed with a 0.2" slit oriented north-south using adaptive optics and the target as a natural guide star. Observations were performed in the CRIRES ELEV mode, which maintains the slit at the same north-south position angle during the whole observing sequence. A standard ABBA nodding sequence was executed using a nodding throw of 12" along the slit and two ABBA nodding cycles. Observations used a wavelength setting centered on 4.780 µm, covering a wavelength range from 4.713 µm to 4.818 µm. The telluric standard star HIP 76829 was observed immediately following the science observations. We provide a summary of the observations in Table 2.

We reduced the data with the CRIRES pipeline version 2.3.1³ and a custom set of IDL routines for improved 1D spectrum-merging from the two nodding positions, accurate telluric correction and wavelength calibration. Nodding sequences were corrected for non-linear effects, flat-fielded, and combined using the CRIRES pipeline. A combined 2D spectrum for the nod A and nod B positions was generated individually. Each combined 2D spectrum was corrected for combination residuals (due to small fluctuations in the sky brightness between nods) by subtracting a background spectrum at each position. This residual background spectrum was obtained by computing at each wavelength the median of two background windows each 20 pixels wide at either sides of the PSF. Before subtraction, the residual background spectrum was smoothed in the wavelength direction with a three-pixel box.

A 1D spectrum was extracted from each combined 2D spectrum of nod A and nod B using the optimal extraction method implemented within the CRIRES pipeline. Bad pixels and cosmic rays in the 1D spectrum of each nod were removed manually using the information of the 1D spectrum of the other nod. The 1D spectra of both nods were merged taking their average. Before merging, the 1D spectrum of nod B was shifted a fraction of a pixel such that the cross-correlation between the 1D spectrum of nod A and nod B was maximized. This was done to correct for small sub-pixel differences in wavelength that are due to the tilt of the spectra in the spatial direction. The merged 1D spectrum was wavelength calibrated using the telluric absorption lines by cross-correlation with a HITRAN atmospheric spectrum of Paranal. The accuracy in the wavelength calibration is 0.15–0.2 km s⁻¹.

A 1D telluric standard star spectrum was obtained from the telluric standard observation following the same procedure as was used for the 1D science spectrum. The science 1D spectrum was then corrected for telluric absorption by dividing it by the 1D spectrum of the standard star. Two adjustments in the 1D standard star spectrum were performed before the telluric

<footnote n="3">³ https://www.eso.org/sci/software/pipelines</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 1. Composite normalized spectrum of the $v = 1 \rightarrow 0$ $^{12}$CO, $^{13}$CO, and $^{18}$O lines. Error bars are $1\sigma$ in each spectrum.</figcaption></figure>

correction. First, the 1D standard star spectrum was shifted in the wavelength direction a fraction of a pixel, such that the cross-correlation with the science spectrum was maximized. Second, the differences on the depth of the telluric lines of the 1D standard star with respect to the 1D science spectrum spectrum were corrected. For this we found the parameter $f$ in

$$I_{\text{STD corrected}} = I_0 e^{-f \tau} \quad (1)$$

that gave the smallest $\chi^2$ statistic between the normalized science spectrum and normalized standard star spectrum. The factor $f$ controls the depth of the atmospheric absorption. The goal is to find the value of $f$ that makes the depth of telluric lines of the standard star the same as in the science spectrum. The $\chi^2$ between the normalized science spectrum and normalized standard star spectrum (thus the value of the factor $f$) was calculated for each chip independently. A region with several unsaturated sky absorption lines and without CO ro-vibrational emission was selected in each chip for this. The optical depth $\tau$ was estimated using

$$\tau = -\ln(I_{\text{STD observed}}/I_0) \quad (2)$$

$$I_0 = \bar{I}_{\text{STD observed}} + 3\sigma. \quad (3)$$

Here $\bar{I}_{\text{STD observed}}$ is median of the standard star flux at the wavelengths within 95% and 100% of the atmospheric transmission and $\sigma$ is the noise in the standard star spectrum in the same wavelength range.

The telluric corrected 1D spectrum was flux calibrated by first normalizing it with a polynomial fit to the continuum and then multiplying the normalized spectrum by the expected flux of the WISE W2 (4.6 µm) magnitude of HD 139614 (5.1 mag, WISE release 2012, Cutri 2012). To convert the magnitude into flux we used the 4.7 µm photometry and the zero points of Johnson (1966)<fnref n="4" />. Errors in the final flux-calibrated spectra are

dominated by slit losses and systematic errors in the telluric correction and are around 20%. Finally, the flux-calibrated 1D spectrum was corrected for the radial velocity (RV) of the star (0.3 ± 2.3 km s⁻¹, Alecian et al. 2013) and the radial velocity due to the rotation of Earth, the motion of Earth around the Sun, and the motion of Earth about the Earth-Moon barycenter, using the velocities given by the IRAF task *rvcorrect* (RV$_{bary}$ = −1 × $V_{\text{helio}}$). Integrated line fluxes, line-profile centers and FWHM were measured in the telluric-corrected 1D spectrum using a Gaussian fit to the line-profiles. The errors on these quantities are the errors on the Gaussian fit.

To produce a merged 2D spectrum, we employed the following procedure. The 2D nod A and nod B spectra were corrected for the tilt of the PSF along the wavelength axis using a second-degree polynomial. The 2D nod B spectrum was shifted by a fraction of a pixel in the wavelength direction with a value equal to the shift found for the 1D spectrum extracted for nod B. A 2D section of ±20 pixels from the PSF center was extracted from the nod A and nod B 2D spectra, and both sections were averaged to obtain a merged 2D spectrum. The merged 2D spectrum was corrected for telluric absorption by diving it by the 1D spectrum of the standard star.

The photocenter (i.e., spectro-astrometric signature) was calculated from the merged 2D spectrum by employing the formalism described by Pontoppidan et al. (2011). The PSF-FWHM as a function of the wavelength was calculated by fitting a Gaussian in the spatial direction of the merged 2D spectrum. We calculated the composite 1D line-profiles, photocenter and PSF-FWHM for each isotopolog by averaging the data of individual detected lines. This was done to increase the signal of the CO line with respect to the continuum. The averaging procedure was performed using the velocity as wavelength scale. The theoretical wavelength center of each transition was used as $v = 0$ km s⁻¹ velocity reference. We selected only emission lines that were not blended with other transitions. For each velocity channel we selected the data in regions with atmospheric transmission higher than 20%, and calculated the average flux when at least three data points were available. Channels with fewer than three data points were defined as NaN to exclude data from regions of poor atmospheric transmission. The error in each channel was defined as the standard deviation of the values in each channel. For further analysis, the composite data were recentered such that the center of the 1D spectrum was at $v = 0$ km s⁻¹, and the 1D spectrum was continuum subtracted and normalized by the peak flux (median of the flux within ±2 km s⁻¹).

## 3. Observational results

We have detected $v = 1 \rightarrow 0$ $^{12}$CO, $^{13}$CO, $^{18}$O, $^{17}$O, and $v = 2 \rightarrow 1$ $^{12}$CO emission lines. The $v = 3 \rightarrow 2$ $^{12}$CO emission lines are not detected in the spectrum. We display a summary of the CO lines detected together with the atmospheric transmission in Fig. 2. In Table A.1, we summarize the observed lines, their centers, integrated fluxes, FWHM and the average line ratios with respect to $1 \rightarrow 0$ $^{12}$CO emission. To keep the notation short in the remaining of the paper, we mean by $^{12}$CO, $^{13}$CO, $^{18}$O emission $v = 1 \rightarrow 0$ $^{12}$CO, $^{13}$CO, $^{18}$O emission unless otherwise specified.

We reached a $3\sigma$ sensitivity of $2 \times 10^{-16}$ erg s⁻¹ cm⁻² for a line width of 3.3 km s⁻¹, and $1.2 \times 10^{-15}$ erg s⁻¹ cm⁻² for a line width of 20 km s⁻¹ (equivalent widths of 0.01 and 0.075 Å respectively). We achieved a spectral resolution of ~3.3 km s⁻¹ ($R \sim 10^5$) as measured in an unresolved unsaturated sky-absorption line. The centers of the CO emission lines in the

<footnote n="4">⁴ The WISE and Johnson (1966) zero-points differ by 10%, which is a value lower than the uncertainties due to slit losses and systematic errors.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Examples of the $v = 1 \rightarrow 0$ $^{12}\text{CO}$, $^{13}\text{CO}$, $\text{C}^{18}\text{O}$, $\text{C}^{17}\text{O}$ and $v = 2 \rightarrow 1$ $^{12}\text{CO}$ lines observed. The *lower panels* display the normalized spectrum of the target (in red) and the spectrum of the telluric standard (in black). The spectra are presented corrected by the radial velocity of HD 139614 and the barycentric velocity. The references for $v = 0 \text{ km s}^{-1}$ are the theoretical wavelengths of each of the transitions. Note that the flux scale is larger for the $v = 1 \rightarrow 0$ $^{12}\text{CO}$ lines. Error bars are $3\sigma$. Several $v = 3 \rightarrow 2$ $^{12}\text{CO}$ lines were covered in the spectra but none were detected. See Table A.1 for a summary of the centers, fluxes, flux upper limits and FWHM of the lines.</figcaption></figure>

barycentric and radial-velocity-corrected spectra are located on average at $v = 2 \pm 1 \text{ km s}^{-1}$ (see Table A.1). As this value is close to zero and is lower than the uncertainty of $\pm 2.3 \text{ km s}^{-1}$ in the radial velocity (Alecian et al. 2013), we conclude that the CO emission is at the stellar velocity and therefore most likely originates in the disk. The $^{12}\text{CO}$ composite line-profile has a flat top and does not display evidence of asymmetries$^5$. The composite $^{13}\text{CO}$ and $\text{C}^{18}\text{O}$ lines are single peaked. Some asymmetric sub-structures are present in both lines but they are consistent with noise.

$\text{C}^{18}\text{O}$ emission is, at the $2\sigma$ level, $4 \text{ km s}^{-1}$ narrower than $^{12}\text{CO}$ emission. $^{13}\text{CO}$ and $2 \rightarrow 1$ $^{12}\text{CO}$ lines are $1 \text{ km s}^{-1}$ narrower than the $^{12}\text{CO}$ line, at the $1\sigma$ level. To further test whether the $^{12}\text{CO}$, $^{13}\text{CO}$ and $\text{C}^{18}\text{O}$ line-profiles are different, we ran a two-sample Kolmogorov-Smirnov (K-S) test$^6$ on the composite 1D spectra in the $\pm 15 \text{ km s}^{-1}$ interval, after normalization by the line peaks. The K-S significance between the $^{12}\text{CO}$

and $\text{C}^{18}\text{O}$ line-profiles is 8%, between the $^{12}\text{CO}$ and $^{13}\text{CO}$ line-profile is 30%, and between the $\text{C}^{18}\text{O}$ and $^{13}\text{CO}$ line-profile is 97%. The K-S test indicates that the $^{12}\text{CO}$ and $\text{C}^{18}\text{O}$ profiles are different (the $\text{C}^{18}\text{O}$ is narrower), and that statistically the $^{13}\text{CO}$ profile resembles the $\text{C}^{18}\text{O}$ profile more than the $^{12}\text{CO}$ profile.

The $1\sigma$ average error obtained in the stacked photocenter is 0.06 pixels, which is equivalent to 5 mas or 0.7 AU at $d = 131$ pc (see Fig. 3). We note that different channels have different error bars and the $1\sigma$ error quoted is an average value. No displacement of the photocenter centroid is detected at the position of the $^{12}\text{CO}$ lines (the interpretation of this constraint requires modeling and is discussed in the next section).

The single-nod PSF-FWHM continuum of HD 139614 ($178 \pm 10$ mas) and the telluric standard ($172 \pm 10$ mas) are consistent within the errors, which means that there is no evidence of extended continuum emission at $4.7 \mu$m. We measured a stacked continuum PSF-FWHM of $2.40 \pm 0.05$ pixels ($1\sigma$), equivalent to $206 \pm 4$ mas or $27 \pm 0.6$ AU at $d = 131$ pc (29 AU at 140 pc) (see Fig. 3). The difference of 30 mas (~1/3 pixel) between the stacked continuum PSF-FWHM and the single-nod PSF-FWHM corresponds to systematic errors introduced during the merging of the 2D nod A and nod B spectra, and to small differences between the PSF-FWHM at the location of the continuum of the different CO transitions. The composite PSF-FWHM at the

<footnote>$^5$ We note that the $1\sigma$ error in the flux in the $^{12}\text{CO}$ composite line-profile is slightly larger at negative velocities. This is because the left side of the line is located in a region with a lower atmospheric transmission. The small differences in the flux between negative and positive velocities in the line are mostly due to differences in the atmospheric transmission.</footnote>

<footnote>$^6$ KSTWO function in IDL.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Observed composite of the normalized $v = 1 \rightarrow 0$ $^{12}$CO line-profile, photocenter, and PSF-FWHM. In red we show the same quantities produced by a flat disk model with a power-law intensity with $R_{in} = 1.2$ AU and $R_{out} = 15$ AU (black cross in Fig. 4). Error bars on the composite line-profile are $1\sigma$. Horizontal dotted lines are the average $1\sigma$ errors in the +20 to +40 km s$^{-1}$ region. Note that this plot assumes a distance of 140 pc for HD 139614. Using the recently announced Gaia distance of $131 \pm 5$ pc, the mean of PSF-FWHM is 27 AU.</figcaption></figure>

location of the line appears constant as a function of the wavelength. This directly indicates that there is no $^{12}$CO emission extending to spatial scales larger than ~30 AU. More stringent limits are deduced in the next section.

## 4. Analysis

We derived constraints on the disk structure from our CRIRES data using models with an increasing complexity. First, we deduce the extent of the CO-emitting region from the composite $^{12}$CO spectrum and spectro-astrometric signature, using a flat Keplerian disk with a parametric power-law intensity. Then, we constrain the average column density of the gas and the temperature of each isotopolog from the rotational diagrams, using an 1D local thermodynamic equilibrium (LTE) slab model with

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** $\chi^2_{red}$ contour plot for the grid of flat disk Keplerian models with a power-law intensity. The black cross displays the model with the lowest $\chi^2_{red}$ (0.35). The yellow and red curves show, for each $R_{in}$, the value of $R_{out}$ that would generate a $1\sigma$ spectro-astrometric signal or a $1\sigma$ PSF-FWHM broadening, respectively.</figcaption></figure>

single temperature and single column density. Finally, we derive the column density and temperature distribution of the gas as a function of the radius from the simultaneous fit of the line-profiles and rotational diagrams of the three CO isotopologs. For this, we use a large grid of 1D flat Keplerian LTE disk models with a power-law temperature and column density distribution. A summary of our analysis strategy is given in Table 3.

### 4.1. Extent of the CO ro-vibrational emitting region

The simplest way to model a line-profile and spectro-astrometric signature and deduce the emitting area is to assume a flat Keplerian disk with a power-law intensity as a function of the radius:

$$I(R) = I_0(R/R_{in})^\alpha, \quad (4)$$

extending from the an inner radius $R_{in}$ to an outer radius $R_{out}$, where $I_0$ is the intensity at $R_{in}$ which is assumed initially to be 1. The exponent $\alpha$ is obtained for each pair of $R_{in}$ and $R_{out}$ such that $I(R_{out}) = 0.01 \times I_0$. In this model, all the physics of the excitation of the line is in the exponent $\alpha$. The 1% limit on the intensity was chosen because the line-profile does not change significantly when integrating to a lower percentage.

We modeled the composite $^{12}$CO line-profile, photocenter, and PSF-FWHM with this simple flat Keplerian disk with parametrized intensity. We provide the details of the model in Sect. B of the Appendix. The model includes the effect of the disk inclination, the effects of the slit width, the spectral broadening due to the CRIRES resolution, and the spatial resolution during the observations. In the models, we used a central stellar mass of $1.7 M_\odot$, an inclination $i = 20^\circ$, a PA = 292$^\circ$ (Matter et al. 2016), and a north-south slit orientation. Models assume a distance of 140 pc as calculations were performed before the recent Gaia distance measurement of $131 \pm 5$ pc.

We calculated a grid of disk models varying $R_{in}$ between 0.1 and 70 AU and $R_{out}$ between 0.2 and 100 AU. In Fig. 4 we present the contour plots of the $\chi^2_{red}$ reduced statistic (assuming
<table><caption>**Table 3.** Types of models used to interpret the observations.</caption><thead><tr><th>Model</th><th>Data modeled</th><th>Constraint</th></tr></thead><tbody><tr><td>1. flat disk with a power-law intensity</td><td>composite <sup>12</sup>CO line-profile, photocenter and <i>PSF-FWHM</i></td><td><ul><li>emitting region</li><li>limits on the gap width</li></ul></td></tr><tr><td>2. 1D LTE slab with single <i>N</i><sub>H</sub> and single <i>T</i><sub>gas</sub></td><td><sup>12</sup>CO, <sup>13</sup>CO, and <i>C</i><sup>18</sup>O rotational diagrams<br><sup>12</sup>CO and <i>ν</i> = 2 → 1 <sup>12</sup>CO rotational diagrams</td><td><ul><li>average <i>N</i><sub>H</sub> and <i>T</i><sub>gas</sub> for each isotopolog</li><li>CO excitation mechanism</li></ul></td></tr><tr><td>3. flat disk with a power-law column density, temperature distribution, and LTE excitation</td><td><sup>12</sup>CO, <sup>13</sup>CO, <i>C</i><sup>18</sup>O rotational diagrams and <sup>12</sup>CO P(9), <sup>13</sup>CO R(4), and <i>C</i><sup>18</sup>O R(6) line-profiles</td><td><ul><li>column density distribution</li><li>temperature distribution</li><li>depth of the gas density drop</li><li>limits on the gas gap width and depth</li></ul></td></tr></tbody></table>

three free parameters: $R_{in}$, $R_{out}$, and $I_0$) of the model of the composite $^{12}$CO line-profile. Disk models with $0.9 < R_{in} < 1.8$ AU, and $13 < R_{out} < 20$ AU gave the best fit to the $^{12}$CO composite line-profile. The model that displays the smallest $\chi^2_{red}$ has $R_{in} = 1.2$ AU and $R_{out} = 15$ AU and an $\alpha$ exponent of the intensity $-1.8$.

Although disk models with $R_{out}$ as large as 20 AU provide a good fit to the $^{12}$CO composite line-profile, star + disk models with a CO-emitting region with $R_{out} > 18$ AU generate a photocenter displacement at the position of the CO line that is larger than the $1\sigma$ limit of the observations (see the yellow curve in Fig. 4, which displays for each $R_{in}$ the value of $R_{out}$ that would generate a displacement of the photocenter by $1\sigma$). In a similar way, star + disk models with a CO-emitting region with $R_{out} > 15$ AU generate a PSF-FWHM broadening at the position of the line that is $1\sigma$ larger than the observations (see orange line in Fig. 4, which displays for each $R_{in}$ the $R_{out}$ that generates a PSF-FWHM larger than $1\sigma$ at the line position). The non-detections of the photocenter displacement and the PSF-FWHM broadening constrain $R_{out}$ to less than 15 AU. The model that best fits the $^{12}$CO line-profile is compatible with the non-detection of the astrometric signature and PSF broadening at the position of the line (see Fig. 3).

Our simple flat-disk model accurately describes the overall line-profile, the line width, and the line wings (emission at $5 < \nu < 15$ km s$^{-1}$). However, the model appears to slightly underpredict the emission at velocities near zero. This suggests that a weak emission component at large radii might be present. Nevertheless, if present, this component does not generate a detectable spectrometric signature. The zero-velocity component could be an additional emission component from the outer disk at $R \ge 6$ AU that is not captured in our simple flat-disk model, or a disk wind emission component as seen in CO ro-vibrational in other protoplanetary disks (e.g., Pontoppidan et al. 2011; Hein Bertelsen et al. 2016). Although presence of a wind cannot be ruled-out completely, the symmetry of the line (i.e. the lack of an emission shoulder in the blue), the lack of a spectroastrometric signature, and the very fact that the line-profile is well described by a disk model suggest that the observed CO emission is consistent with disk emission.

### 4.1.1. The inner radius of the CO emission

The models that best reproduce the $^{12}$CO composite line-profile have an inner radius around 1 AU. However, some models with smaller $R_{in}$ are also compatible with the data. In Figs. B.1a,b we display the line-profiles expected for disks with $R_{in}$ ranging from 0.1 to 1.2 AU. In panel (a) we show the results of the models with $R_{out}$ fixed to 15 AU ($\alpha$ adjusted such that $I(R_{out}) = 0.01 \times I(R_{in})$).

In panel (b) the line-profiles with $\alpha$ fixed to $-1.8$ ($R_{out}$ set such that $I(R_{out}) = 0.01 \times I(R_{in})$). Depending on the value of $\alpha$, models with $R_{in}$ as low as 0.3 AU can be compatible with the observed $^{12}$CO composite line-profile.

In all our power-law intensity models, we have assumed a sharp inner edge, thus an abrupt increase in the intensity from zero to $I_0$ at $R_{in}$. If instead we assume a soft inner edge, thus a smooth increase of the intensity from $R_{in}$ up to the radius of the maximum intensity $R_{I_{max}} = 1.2$ AU, then $R_{in}$ can be as small as 0.01 AU and the line profile would still be compatible with the data (see Fig. B.1c). The $R_{in}$ constraint from a power-law intensity model with a sharp inner edge corresponds to the radius of the maximum intensity. CO gas can still be present farther in if the inner edge is soft. In Sect. 4.5 we provide upper limits to the gas column density at $R < 1$ AU based on the $^{12}$CO line-profile shape.

### 4.1.2. A continuous or a gapped gas distribution?

The $^{12}$CO line-profile data is well described by a continuous and smooth intensity profile from 1.2 AU up to 15 AU. As mentioned in the introduction, Matter et al. (2014, 2016) resolved a gap in the dust from 2.5 AU to 6 AU based on near- and mid-IR VLT interferometric observations. This raises the question whether the $^{12}$CO line-profile could be described by an intensity distribution with a gap.

The $^{12}$CO line-profile clearly indicates that there is emission at $R < 6$ AU, otherwise the line-profile would have been much narrower (see the blue line in panel a of Fig. 5). As a consequence, an inner gas hole of 6 AU radius is ruled out. Furthermore, the line-profile rules out a CO-emitting region confined to a narrow ring between 1.2 and 2.5 AU, otherwise the line would have been much broader (see the yellow curve in panel a of Fig. 5).

We have tested the scenario in which the intensity distribution of the best solution of the power-law intensity model has a gap (i.e., no emission) between 2.5 AU and 6 AU. In this case (Fig. 5a), the velocity channels between 3 and 8 km s$^{-1}$ are not well reproduced. Such a large gap of 3.5 AU is not compatible with the observed $^{12}$CO line. If the gap in the gas is smaller than 2 AU, then the line-profile could be consistent with the observations (Fig. 5b). Given the CRIRES resolution, a small gap of 1–2 AU in the intensity would not be detectable. In Sect. 4.4 we derive constraints on the CO column density inside a potential gap.

### 4.2. Average temperature and column density

The detection of ro-vibrational emission of the CO isotopologs $C^{18}$O and $C^{17}$O indicates that the emitting medium must be
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Line-profiles predicted for models with a gap in the intensity distribution: **a)** line-profiles expected for a intensity distribution with a gap devoid of gas from 2.5 AU to 6 AU, in orange the contribution from 1.2 to 2.5 AU, in blue the contribution from 6 to 15 AU, and in red the total line-profile; **b)** similar plot but for a gap devoid of gas of width 1 AU extending from 5 to 6 AU. Error bars on the spectrum are $1\sigma$.</figcaption></figure>

dense and warm. To constrain the average temperature and column density of the gas probed in the CRIRES spectra, we modeled the ¹²CO, ¹³CO and C¹⁸O line fluxes using a simple semi-infinite slab model in LTE with a single gas temperature and column density. We did not model the C¹⁷O observations because only two line fluxes are available and no reliable estimate of the temperature can be derived. We wrote a CO LTE slab model using the frequencies, energy levels, and Einstein coefficients from Chandra et al. (1996).

We generated a grid of LTE slab models by varying the hydrogen-nuclei column density $N_H$ from $10^{18}$ to $10^{25}$ cm⁻² with steps of 0.25 dex, and the gas temperature $T_g$ from 100 to 1000 K with steps of 25 K. We assumed a turbulent line broadening of 0.1 km s⁻¹. We used a ¹²CO abundance of $10^{-4}$ ($N_{12\text{CO}}$ = 1.0 × $10^{-4}N_H$), a ¹²CO/¹³CO ratio of 100 and a ¹²CO/C¹⁸O ratio of 690 following Smith et al. (2009)<fnref n="7" />.

The output of a slab model is in units of line flux per stardian and the optical depth of each transition. To compare slab calculations with the observed line fluxes, an average solid angle of the emitting region needs to be prescribed. A model with a single temperature and a single column density model is equivalent to assuming a disk model with constant intensity with radius (i.e., $\alpha = 0$). We tested Keplerian disk models with a constant intensity and found that the line-profile, photocenter and PSF-FWHM could be reproduced by a flat disk of $R_{\text{in}} = 0$, $R_{\text{out}} = 6.5$ AU. We therefore used an average emitting region radius of 6.5 AU and a distance of 140 pc<fnref n="8" /> to determine the solid angle for the 1D slab model<fnref n="9" />.

For each $N_H$ and $T_{\text{gas}}$ slab model, a rotational diagram for each CO isotopolog was calculated. We used the X and Y coordinates $Y = \ln(F_{\text{ul}}/\text{vg}_u A_{\text{ul}})$ and $X = E_u$. Here $F_{\text{ul}}$ is the line flux of the transition between the upper level u and the lower level l, $g_u$ is the degeneracy of the upper level (2J+1), $A_{\text{ul}}$ is the Einstein coefficient of transition and $E_u$ the upper energy level of the transition. For the rotational diagram of each isotopolog we calculated the statistical quantity $\chi_{\text{red}}^2 = \frac{1}{N-4} \sum_i (Y_{\text{slab } i} - Y_{\text{obs } i})^2 / \sigma_{Y_{\text{obs } i}}^2$. Here $\sigma_{Y_{\text{obs } i}}$ is the difference between Y calculated using the observed line flux and Y calculated using the observed line flux plus $3\sigma$. The $N-4$ corresponds to three degrees of freedom ($N_H$, $T_{\text{gas}}$, $\Omega$), and $N$ the number of data points (8 for ¹²CO, 7 for ¹³CO and 9 for C¹⁸O). In Fig. 6 we display the $\chi_{\text{red}}^2$ contour plots for each CO isotopolog and one $\chi_{\text{red}}^2$ combining the data of the three isotopologs.

We find that the best fit to the rotational diagram is different for each CO isotopolog. ¹²CO is best described with $T_{\text{gas}} \sim 450$ K and $N_H > 10^{20}$ cm⁻². ¹³CO is best reproduced by lower temperatures ($T_{\text{gas}} \sim 380$ K) and $N_H$ column densities of at least $5 \times 10^{22}$ cm⁻². The C¹⁸O is best fit by even lower temperatures ~350 K and $N_H$ column densities higher than $5 \times 10^{22}$ cm⁻². The errors on $T_{\text{gas}}$ are 10–20 K. The colder temperatures for ¹³CO and C¹⁸O emission indicate that they are produced at larger radii than the ¹²CO emission or at lower vertical scale heights.

The best-fit combined model for the rotational diagrams of the three isotopologs emission is a model with $N_H = 5 \times 10^{21}$ cm⁻² and $T_{\text{gas}} = 450$ K. In this model the ¹²CO emission is optically thick and the ¹³CO and C¹⁸O optically thin (Fig. 7). The combined solution only satisfactorily describes the ¹²CO rotational diagram, however. The slopes of the ¹³CO and C¹⁸O rotational diagrams are not well reproduced. The curvature on the rotational diagrams suggests that there ¹³CO and C¹⁸O emissions are optically thick or that there is a gradient in temperature. This is also suggested by the $\chi_{\text{red}}^2$ contours, because solutions are degenerate with respect to $N_H$, giving only lower limits for the column density.

<footnote n="7">⁷ The ¹²CO/¹³CO and ¹²CO/C¹⁸O ratios of Smith et al. (2009) were deduced from high-resolution spectroscopy of CO ro-vibrational lines detected in absorption toward VV CrA, a binary T Tauri star in the Corona Australis molecular cloud. The ¹²C/¹³C from Smith et al. is nearly twice the expected interstellar medium (ISM) ratio, and the ¹²CO/C¹⁸O ratio is ×1.4 the ISM ratio. We used the ¹²CO/¹³CO and ¹²CO/C¹⁸O ratios of Smith et al. instead of the ISM ratios because we consider that they are</footnote>

<footnote>more representative of the isotopologs ratios that would be expected in the inner disk of HD 139614.</footnote>

<footnote n="8">⁸ Models were calculated before the recent Gaia distance measurement of 131 ± 5 pc, the difference in distance of 5–9 pc does not change the conclusions reached.</footnote>

<footnote n="9">⁹ We note that the $R_{\text{out}}$ of 15 AU found in Sect. 4.1 cannot be used here because this emitting region was found with a decreasing power-law intensity, which is equivalent to assuming a radial decreasing temperature and column density distribution. The 1D slab model has a single temperature and column density.</footnote>
<figure><img src="image_9.png" /></figure>

**Fig. 6.** $\chi^2_{\text{red}}$ contour plot of the modeled rotational diagrams for $^{12}\text{CO}$, $^{13}\text{CO}$, and $\text{C}^{18}\text{O}$ using a single temperature and surface density slab model. The combined $\chi^2_{\text{red}}$ contour plot is obtained by simultaneously using all the data of the three isotopologs. The cross indicates the location of the $\chi^2_{\text{red}}$ minimum in each panel. The numbers inside the contour indicate n-times the value of the minimum $\chi^2_{\text{red}}$.

### 4.2.1. Thermally excited or UV-pumped emission?

The detection of $v = 2 \rightarrow 1^{12}\text{CO}$ emission raises the question whether the observed CO ro-vibrational emission is thermally excited or if it is due to UV-pumping as observed in some other Herbig Ae/Be stars (e.g., Brittain et al. 2007; van der Plas et al. 2015). The $v = 2 \rightarrow 1^{12}\text{CO}$/$v = 1 \rightarrow 0^{12}\text{CO}$ average line ratio of ~0.2 is lower than in Herbig Ae/Be stars with flared disks, but it is at the higher end of the T Tauri sample with single-component CO ro-vibrational emission (see Table 3 in Banzatti & Pontoppidan 2015). The non-detection $v = 3 \rightarrow 2^{12}\text{CO}$ emission and the upper limit of 0.04 on the $v = 3 \rightarrow 2^{12}\text{CO}$/$v = 1 \rightarrow 0^{12}\text{CO}$ line ratio indicate that the UV-pumping, if present, is not as strong as in other Herbig Ae/Be stars with CO fluorescent emission. For example, in the case of HD 100546, van der Plas et al. (2015) measured an average $v = 3 \rightarrow 2^{12}\text{CO}$/$v = 1 \rightarrow 0^{12}\text{CO}$ line ratio of 0.3. The A7V spectral type of HD 139614 is later than most of the Herbig Ae/Be stars studied in Brittain et al. (2007) and van der Plas et al. (2015). This might in part explain the weaker effect of UV-pumping in HD 139614 with respect to other Herbig Ae/Be stars studied previously.

We explored models around the best solution of the combined fit to the $v = 1 \rightarrow 0^{12}\text{CO}$, $^{13}\text{CO}$, and $\text{C}^{18}\text{O}$ rotational diagrams and calculated the LTE $v = 2 \rightarrow 1^{12}\text{CO}$ and $v = 3 \rightarrow 2^{12}\text{CO}$ emission. We found that the observed $v = 2 \rightarrow 1^{12}\text{CO}$ emission and the non-detection of the $v = 3 \rightarrow 2^{12}\text{CO}$ emission can be well described by an LTE model, with $T_{\text{gas}} = 460$ K and $N_{\text{H}} = 10^{22} \text{ cm}^{-2}$, in which the rotational temperature is equal to the vibrational temperature$^{10}$. We

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

**Fig. 7.** Rotational diagrams and optical depths of the slab model ($N_{\text{H}} = 5 \times 10^{21} \text{ cm}^{-2}$, $T_{\text{gas}} = 450$ K) with the lowest $\chi^2_{\text{red}}$ combining all the $^{12}\text{CO}$, $^{13}\text{CO}$, and $\text{C}^{18}\text{O}$ data (cross in Fig. 6). In the left panels observations are shown in black and models in red. In the right panels the optical depths of the observed transitions are shown in black, other transitions are plotted in red. A single-temperature and single-density slab does not correctly describe the rotational diagrams of the three CO isotopologs simultaneously. Error bars on the rotational diagram are $1\sigma$.

show this model in Fig. 8, where we display the rotational diagrams of $v = 1 \rightarrow 0^{12}\text{CO}$ and $v = 2 \rightarrow 1^{12}\text{CO}$ emission of the model and the observations. This model generates $v = 3 \rightarrow 2^{12}\text{CO}$ emission lines with integrated fluxes on the order of $10^{-18}-10^{-17} \text{ erg s}^{-1} \text{ cm}^{-2}$, which is consistent with the non-detection of these lines in our CRIRES spectra. We conclude that the observed CO emission is most likely thermally excited.

### 4.3. Deriving the surface density and temperature distribution

The different temperatures and column densities of each CO isotopolog and the difference on line widths between the CO isotopologs, suggest that the emission of each isotopolog is produced at different radial distances and/or vertical heights. The narrower and colder $^{13}\text{CO}$ and $\text{C}^{18}\text{O}$ lines pose an interesting puzzle for the interpretation of the data. If the CO ro-vibrational emission is modeled with a 1D power-law column density and temperature distribution ($N_{\text{H}} \propto R^{\alpha_{\text{H}}}$ and $T_{\text{gas}} \propto R^{\alpha_{\text{gas}}}$) and both distributions have no discontinuities (i.e. no gaps nor density drops or jumps in the temperature), and if the $^{12}\text{CO}/^{13}\text{CO}$ and $^{12}\text{CO}/\text{C}^{18}\text{O}$ abundance ratios are constant, then, a disk model able to reproduce the $^{12}\text{CO}$ line-profile and $^{12}\text{CO}$ rotational diagram would generate $^{13}\text{CO}$ and $\text{C}^{18}\text{O}$ lines that are too broad, too strong, and too

<footnote>$^{10}$ The model with $T_{\text{gas}} = 450$ K and $N_{\text{H}} = 10^{22} \text{ cm}^{-2}$ of Fig. 7 generated slightly weaker $v = 2 \rightarrow 1^{12}\text{CO}$ emission.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 8. Rotational diagrams of the $v = 1 \to 0$ $^{12}$CO (squares) and $v = 2 \to 1$ $^{12}$CO (triangles) observed emission. Overplotted is a slab model with $T_{\text{rot}} = T_{\text{vib}} = 460$ K and $N_{\text{H}} = 10^{22} \text{ cm}^{-2}$, in red for $v = 1 \to 0$ $^{12}$CO emission and in green for $v = 2 \to 1$ $^{12}$CO emission.</figcaption></figure>

warm to be consistent with the observations (see Fig. 9). In the next sub-section we provide the details of the model. Since $^{13}$CO and $^{18}$O emission is optically thin up to relatively high column densities ($N_{\text{H}} \sim 10^{22} \text{ cm}^{-2}$ and $N_{\text{H}} \sim 10^{23} \text{ cm}^{-2}$ respectively), the 1D single power-law models generate $^{13}$CO and $^{18}$O lines that are too broad, too strong, and too warm because the column of gas at small radii is too large.

### 4.3.1. Power-law temperature and column density Keplerian disk model with a depleted inner region

Near- and mid-IR interferometry and the SED (Matter et al. 2016) indicate a depletion of dust mass on the order of $10^{-4}$ in the inner 6 AU with respect to the extrapolated surface density at $R > 6$ AU. If a similar behavior were also be followed by the gas, the narrower and colder $^{18}$O and $^{13}$CO profiles could be explained as the result of a gas density drop in the inner 6 AU of the disk. A gas density drop can cause the $^{18}$O and $^{13}$CO lines to be optically thin in the inner 6 AU and cause their emission be dominated by the contribution at $R > 6$ AU where the density is higher and the gas colder.

To test the hypothesis of the gas density drop, we modeled the observed CO ro-vibrational emission with a flat disk in Keplerian rotation, in which the gas column density and temperature are described by a power-law distribution:

$$N_{\text{H}}(R \ge 6 \text{ AU}) = N_{\text{H}}(R = 6 \text{ AU}) \left(\frac{R}{6 \text{ AU}}\right)^{\alpha_{\text{N}_{\text{H}} \text{ outer}}} \quad (5)$$

$$T_{\text{gas}}(R \ge 1 \text{ AU}) = T_{0}(R = 1 \text{ AU}) \left(\frac{R}{1 \text{ AU}}\right)^{\alpha_{T_{\text{gas}}}}. \quad (6)$$

We allowed the models to have a reduced surface density at $R < 6$ AU by a factor $\delta_{\text{gas}}$

$$N_{\text{H}}(R < 6 \text{ AU}) = \delta_{\text{gas}} \cdot N_{\text{H}}(R = 6 \text{ AU}) \left(\frac{R}{6 \text{ AU}}\right)^{\alpha_{\text{N}_{\text{H}} \text{ inner}}}. \quad (7)$$

A reference radius of 6 AU was set for the gas column density to compare the gas distribution with that of the dust. A reference radius of 1 AU was selected for the temperature, given that the modeling of the $^{12}$CO line with a power-law intensity suggested that 1 AU is the radius of the maximum intensity.

We chose $\alpha_{N_{\text{H} \text{ inner}}}$ to range from -2.5 up to +3 to cover a wide range of possible surface density distributions in the inner disk$^{11}$. As Matter et al. (2016) found a dust depletion factor of $10^{-4}$, we tested models with $\delta_{\text{gas}}$ ranging from 1 to $10^{-4}$.

The choice of modeling the temperature as a power-law instead of using a full radiative transfer calculation (e.g., Woitke et al. 2009; Thi et al. 2013; Carmona et al. 2014; Bruderer 2013; Bruderer et al. 2014) enables us to describe the CO temperature in the emitting region independently of that of the dust with a minimum number of free parameters. This permitted us to explore a large portion of the parameter space.

The disk was modeled with a flat geometry using a radial and azimuthal grid. The model is analogous to the model described in Sect. 4.1 and Sect. B, with the difference that the intensity at each radius was calculated using the local $T_{\text{gas}}$ and $N_{\text{H}}$ using the CO slab model previously described,

$$I(R) = I(T_{\text{gas}}, N_{\text{H}})_{\text{slab}}. \quad (8)$$

The local broadening of the line is the convolution of the turbulent broadening ($0.1 \text{ km s}^{-1}$), the local thermal broadening and the spectral resolution$^{12}$. We recall that the CO slab model assumes LTE excitation, which is a good approximation because CO emission is most likely thermally excited. We set the outer radius equal to 30 AU. Observations set an upper limit of 15 AU to the emitting region, but, we used a larger outer radius to permit some combinations of $N_{\text{H}}(R)$ and $T_{\text{gas}}(R)$ inside the grid to have a sufficiently large radial extent to let the intensity decrease to low levels. A model in which the radial calculation grid ends artificially early would generate a line-profile and a line flux that does not correctly represent the selected $N_{\text{H}}(R)$ and $T_{\text{gas}}(R)$. In the flat-disk parametric intensity models that describe the CO emission in HD 139614, we saw no significant change in the line-profiles with or without slit (lines are dominated by the contribution in the inner 15 AU). Therefore, we did not include the slit effects in our model to enable the calculation of a large number of models.

A model has six free parameters: $N_{\text{H}}(R = 6 \text{ AU})$, $\alpha_{N_{\text{H} \text{ inner}}}$, $\alpha_{N_{\text{H} \text{ outer}}}$, $\delta_{\text{gas}}$, $T_{(R = 1 \text{ AU})}$ and $\alpha_{T_{\text{gas}}}$. Each model produces integrated line-fluxes and rotational diagrams for the three CO isotopologs and synthetic line-profiles at the CRIRES resolution for the $^{12}$CO P(9) line at 4745.13 nm, $^{13}$CO R(4) at 4730.47 nm, and the $^{18}$O R(6) line at 4724.03 nm. These three CO transitions were selected because their line-profiles have a high S/N and are less affected by telluric absorption. The merged composite line-profile was not used because the model line predictions needed to be compared with a line-profile in flux units. In the merged composite spectra the flux information is lost.

To find the models that best describe the observations, we ran a uniform grid of 81 000 models covering the parameter space described in Table 4. To calculate the most probable values of

<footnote>$^{11}$Previous models of CO ro-vibrational emission in transition disks have assumed a flat surface density (Pontoppidan et al. 2008, 2011) or a surface density decreasing with a -3/2 exponent (Salyk et al. 2009). Carmona et al. (2014) found an increasing surface density profile for the gas with exponent +0.2 in the inner disk of the transition disk HD 135344B. Andrews et al. (2011), Bruderer et al. (2014), and van der Marel et al. (2015b) modeled SMA or ALMA observations of transition disks assuming the same surface density exponent of -1 for the inner and the outer disk. Matter et al. (2016) found that the surface density of the dust in the inner 2.5 AU of HD 139614 increases radially with an exponent +0.6.</footnote>

<footnote>$^{12}$Including the instrument resolution in the convolution kernel of the thermal and turbulent broadening enable saving hundreds of convolutions per model, $10^7$ convolutions in the grid, and it is equivalent to convolving the final data-cube by the instrument resolution.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Examples of the predicted ¹³CO R(4) and C¹⁸O R(6) emission for a continuous power-law distribution of temperature and column density that describe the ¹²CO P(9) line-profile and the ¹²CO rotational diagram. Colors in the spectra and rotational diagrams correspond to the different column density distributions in the first panel. A disk model with a single power-law surface density cannot simultaneously reproduce the ¹²CO, ¹³CO, and C¹⁸O line-profiles and rotational diagrams.</figcaption></figure>

<table><caption>**Table 4.** Parameter space of the power-law $N_H$ and $T_{gas}$ models.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Values</th></tr></thead><tbody><tr><td>$N_H$ (R = 6 AU)</td><td>[cm<sup>-2</sup>]</td><td>10<sup>20</sup>, 10<sup>21</sup>, 10<sup>22</sup>, 10<sup>23</sup>, 10<sup>24</sup>, 10<sup>25</sup></td></tr><tr><td>$\alpha_{N_{H,\text{inner}}}$(R &lt; 6 AU)</td><td></td><td>-2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr><tr><td>$\alpha_{N_{H,\text{outer}}}$(R &gt; 6 AU)</td><td></td><td>-2.5, -2.0, -1.5, -1.0, -0.5</td></tr><tr><td>$\delta_{\text{gas}}$(R = 6 AU)</td><td>[K]</td><td>1, 10<sup>-1</sup>, 10<sup>-2</sup>, 10<sup>-3</sup>, 10<sup>-4</sup></td></tr><tr><td>$T_0$ (R = 1 AU)</td><td></td><td>550, 575, 600, 625, 650, 675, 700, 725, 750</td></tr><tr><td>$\alpha_{T_{gas}}$</td><td></td><td>-0.40, -0.35, -0.30, -0.25, -0.20</td></tr><tr><td>Total</td><td></td><td>81 000 models</td></tr></tbody></table>

the parameters in the grid, we used a Bayesian approach (see for example Pinte et al. 2007, and references therein). We provide details of the calculation of the Bayesian probabilities in Appendix C.

### 4.3.2. Grid results

Figure 11 displays the Bayesian probability distribution diagrams of the grid. The first panel in the upper left displays the probability of models with and without a gas density drop. The diagram clearly shows that a gas column density drop in the inner 6 AU of the disk is required to simultaneously reproduce the CO ro-vibrational line-profiles and the rotational diagrams of the three isotopologs. A $\delta_{gas} = 10^{-2}$ in the column density appears as the most likely value. To directly illustrate this, we show the progression from a model without a column density drop to the best-fit grid model which has a column density drop of $10^{-2}$ in Fig. 12. The empirical evidence of the gas density drop emerges from both the line-profile shapes and the rotational diagrams. Models without a gas density drop generate ¹³CO and C¹⁸O lines that are too broad, strong and warm (too much warm

¹³CO and C¹⁸O emitted at small radii) to be consistent with the observations.

The column density of gas traced by CO at R = 6 AU is well constrained to $N_H \sim 10^{23}$ cm⁻². This gas column density is similar to the dust column density at R = 6 AU found by Matter et al. (2016)<fnref n="13" />. The gas column density at R = 6 AU is higher than the column density found in the single $T_{gas} - N_H$ slab model of the three CO isotopologs of Sect. 4.2. This is because higher column densities describe the C¹⁸O and ¹³CO rotational diagrams better. In fact, the emission of C¹⁸O and ¹³CO is optically thick in the 6–10 AU region where 80 to 90 % of the line flux is emitted (see Fig. 13).

CO ro-vibrational emission traces the gas in regions where the dust is optically thin or the disk upper layers where $T_{gas} > T_{dust}$. In Fig. 14 we display the dust optical depth at 4.7 µm from the best-fit model of the SED and IR-interferometry data of HD 139614 from Matter et al. (2016). At R < 6 the dust is optically thin at 4.7 µm down to the disk mid-plane, at R ≥ 6 AU the dust is optically thick at 4.7 µm (except in the disk surface layer). As a consequence, in the inner disk at R < 6 AU, the gas column density traced by CO ro-vibrational emission should be a good estimate of the total column density of gas. Our models suggest that the column density of gas at 1 < R < 6 AU ranges between $N_H = 3 \times 10^{19}$ cm⁻² and $N_H = 10^{21}$ cm⁻² ($7 \times 10^{-5} - 2.4 \times 10^{-3}$ g cm⁻²). In the outer disk at R > 6 AU, the gas column density traced by CO ro-vibrational emission is a lower limit, as we trace only the gas in the disk surface where the dust is thin and $T_{gas} > T_{dust}$. If we assume a gas-to-dust mass ratio of 100 for the outer disk and use the dust column density at R ≥ 6 of Matter et al. (2016), then the total column of gas at R ≥ 6 AU can be up to a factor 100 higher than the column density traced by CO ro-vibrational emission. Therefore, the gas density drop $\delta_{gas}$ could be as large as $10^{-4}$ depending on the total gas mass of the outer disk.

<footnote n="13">¹³ We fixed the density drop location at 6 AU in the models, see the discussion section for the models with a varying radius of the density drop.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Schematic illustration of the free parameters in the flat Keplerian disk model with a power-law column density and temperature. The grid includes models with decreasing ($\alpha_{\text{inner}} < 0$), flat ($\alpha_{\text{inner}} = 0$) and increasing ($\alpha_{\text{inner}} > 0$) surface density distributions in the inner 6 AU. The anchor point for the surface density is at 6 AU and for the temperature is at 1 AU. In all the models the cavity radius is at 6 AU.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Bayesian probability distributions of the grid of the $N_H$ and $T_{\text{gas}}$ power-law models calculated (81 000 models). The $\chi^2$ statistic used to calculate the probability considers the rotational diagrams and line-profiles of the three isotopologs simultaneously.</figcaption></figure>

The Bayesian probability distributions indicate that the preferred values for the surface density exponent in the inner 6 AU are flat or positive, with $\alpha_{N_{\text{H,inner}}}$ ranging from 0.0 up to high values such as +3.0. A few models with negative inner disk power-law exponents can describe the data, but, the largest fraction of models describing the observations have flat or positive surface density exponents.

The probability plots show that increasingly negative exponents of the density in the outer disk $\alpha_{N_{\text{H,outer}}}$ have a higher probability. This behavior is due to the fact that the models that best reproduce the ¹³CO and C¹⁸O emission are those in which a large portion of the line-flux (~60%) is produced between 6 and 10 AU (see Fig. 13). This suggests that the emission of ¹³CO and C¹⁸O is most likely dominated by the contribution of the inner rim of the outer disk. This naturally explains the narrower line-profiles of these two isotopologs.

Figure 13 displays the surface density profile, temperature profile, optical depth, flux density, cumulative line flux,

rotational diagrams, and the ¹²CO P(9), ¹³CO R(4) and C¹⁸O R(6) line-profiles of the model in the grid that shows the best combined fit to the data. This model has an $N_H$ at R = 6 AU of 10²³ cm⁻², $\alpha_{N_{\text{H,inner}}}$ (R < 6 AU) = +2.0, $\alpha_{N_{\text{H,outer}}}$ (R ≥ 6 AU) = -2.5, $\delta_{\text{gas}}$ (R = 6 AU) = 10⁻², $T_0$ (R = 1 AU) = 675 K, and $\alpha_{T_{\text{gas}}} = -0.35$. We note that the solution is not unique, and models with other parameters can still provide a satisfactory fit to the data. The Bayesian probability distribution diagrams show which values of the model parameters are the most likely based on the observations.

From the pure modeling point of view, the higher probability of models with a flat or with an increasing surface density in the inner disk can be easily explained. The emission of ¹³CO and C¹⁸O is optically thin at R < 6 AU (see upper-right panel in Fig. 13). Thus to have a weaker flux at v > 10 km s⁻¹, a small column of gas is needed at small radii. A surface density with a decreasing profile in the inner 6 AU, even with a low column density (see Fig. 9), produces ¹³CO and C¹⁸O line-profiles with
<figure><img src="image_15.png" /><figcaption>2. Same model but with $\delta_{\text{gas}} = 10^{-1}$: $^{13}$CO and $^{18}$O emission are weaker, but still too strong</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>3. Same model but with $\delta_{\text{gas}} = 10^{-2}$: $^{13}$CO and $^{18}$O lines have a flux closer to observations, but $^{13}$CO emission is still too broad because the amount of gas at small radii.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>4. Change of the $\alpha$ exponent of the inner disk from -1.0 to 0 (flat inner disk): $^{13}$CO and $^{18}$O line-profiles and rotational diagrams are more similar to the observations.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>5. Change of the $\alpha$ exponent of the inner disk to +1.0: improvement of the fit of the $^{13}$CO and $^{18}$O line-profiles and rotational diagrams.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>6. Change of the $\alpha$ exponent of the inner disk to +2.0: improvement of the fit of the wings of the $^{12}$CO line and the $^{12}$CO rotational diagram.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>7. Change of the $\alpha$ exponent of the outer disk from -1.0 to -2.5: better description of the peak flux and rotational diagram of $^{18}$O emission.</figcaption></figure>

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

**Fig. 12.** Example of a progression of disk models taken from the grid to illustrate the effect of the change of parameters in the line-profiles and rotational diagrams. The observed data are plotted in black, and the model predictions in red. Error bars on the line-profiles are $3\sigma$. The dashed line is the dust surface density from Matter et al. (2016). The model starts with a continuous gas disk that follows the same power-law as the dust in the outer disk, and it is refined until the best description of the CO ro-vibrational observations is reached. The two branches seen in the rotational diagram correspond to the R and P branches of CO ro-vibrational emission.

line wings that are too strong. The effect is also seen in the $^{13}$CO and $^{18}$O rotational diagrams. Flat and increasing surface density profiles provide colder rotational diagrams that better describe the observations. A visualization of this is provided in Fig. 12, which illustrates the effect of changing the exponent of the surface density profile in the inner disk from -1.0 to +2.0 while keeping a constant $\delta_{\text{gas}} = 10^{-2}$ at R = 6 AU. As soon as $\alpha_{N_{\text{H inner}}} = 0$ is reached, the strong high-velocity line wings of the $^{13}$CO and lines $^{18}$O disappear. The $^{12}$CO line-profile hardly changes in all the models with $\alpha_{N_{\text{H inner}}} -1.0$ to +1.0 because it is optically thick. However, when $\alpha_{N_{\text{H inner}}}$ is higher than +1.5, $^{12}$CO

becomes optically thin at R < 2 AU, the wings of the $^{12}$CO line become weaker, and the $^{12}$CO line-profile is well fitted.

4.4. Quantitative constraints on the width and column density depth of a gas gap

As discussed in Sect. 4.1.2, the $^{12}$CO ro-vibrational composite line-profile does not display evidence of a gap (i.e. a zone devoid of gas) of size larger than 2 AU. To derive quantitative limits on the gap width and the column density of the gas that could be
BEST-FIT GRID MODEL: $\alpha_{N_{H inner}} = +2.0 \ \delta_{gas} = 10^{-2}$

<figure><img src="image_22.png" /><figcaption>Fig. 13. Gas column density, temperature, CO optical depth, flux density, cumulative line flux, rotational diagrams, and line-profiles of the ¹²CO P(9), ¹³CO R(4), and C¹⁸O R(6) emissions of the best-fit grid model. The model is shown in red, and the observations in black. Observed line-profiles are displayed in flux units after continuum subtraction with 3σ error bars. The two branches seen in the rotational diagram correspond to the R and P branches of CO ro-vibrational emission. The rightmost panels compare the normalized theoretical line-profiles with the observed composite line-profile of each CO isotopolog with a 1σ error bar.</figcaption></figure>

inside a potential gap, we calculated the expected ¹²CO P(9) line for a series of models around the best-fit grid model for gaps of increasing width (i.e., from 5 to 6 AU, from 4 to 6 AU, from 3 to 6 AU, and from 2 to 6 AU) and varying $N_H$ inside the gap (at R = 6 AU) from $10^{21}$ to $10^{17}$ cm⁻². The normalized theoretical ¹²CO P(9) spectrum was compared with the high S/N composite ¹²CO spectrum.

Figure 15 displays the results. The models confirm the suggestion of the simple power-law intensity model (Sect. 4.1.2). A gap of 2 AU or smaller remains undetected. Gaps of width larger than 2 AU and with an $N_H$ inside the gap lower than $10^{18}$ cm⁻² (yellow and red, $\delta_{gas} < 10^{-5}$) would have been seen in the ¹²CO composite spectrum as a line-profile with shoulders at $\pm 15$ km s⁻¹.

### 4.5. Upper limits to the gas column density at R < 1 AU

¹²CO ro-vibrational emission requires relatively low column densities ($N_{CO} \sim 10^{15}$ cm⁻²) to be optically thick. Therefore, if the gas is sufficiently warm, CO ro-vibrational emission is relatively easy to detect. As the gas in the inner 1 AU of a

disk around a Herbig Ae star has a temperature warmer than 300 K, the lack of strong CO ro-vibrational emission in the inner 1 AU of HD 139614 suggests a low column density of gas at R ≤ 1 AU. To derive upper limits to the column of gas at R ≤ 1 AU in HD 139614, we calculated the expected emission from gas between 0.1 and 1.0 AU assuming a flat surface density profile and the temperature profile of the best-grid model. We found that gas column densities higher than $N_H = 5 \times 10^{19}$ cm⁻² ($1.2 \times 10^{-4}$ g cm⁻²) would have produced line-profile wings ($v > 15$ km s⁻¹) stronger than 5σ the noise of the ¹²CO P(9) line (Fig. 16). As we assume standard abundances, our CRIRES observation set a 5σ upper limit to the CO column at R ≤ 1 AU of $N_{CO} = 5 \times 10^{15}$ cm⁻².

UV photodissociation could be responsible for the destruction of CO in the inner 1 AU of the disk around HD 139614. In the absence of dust, CO self-shields against photodissociation in the vertical and radial direction if $N_{CO} > 10^{15}$ cm⁻² (van Dishoeck & Black 1988). Therefore, from the self-shielding perspective, the absence of CO ro-vibrational emission from R < 1 AU suggests that $N_H \le 10^{19}$ cm⁻² at R ≤ 1 AU, assuming standard abundances. This value is consistent
<table><caption>**Table 5.** Parameters of the best-fit grid model.</caption><thead><tr><th>Model</th><th>N<sub>H</sub> (R=6 AU) [cm<sup>-2</sup>]</th><th>α<sub>N<sub>H</sub>inner</sub></th><th>α<sub>N<sub>H</sub>outer</sub></th><th>δ<sub>gas</sub></th><th>T<sub>(R=1 AU)</sub> [K]</th><th>α<sub>T<sub>gas</sub></sub></th></tr></thead><tbody><tr><td>Best-fit grid model</td><td>10<sup>23</sup></td><td>+2.0</td><td>-2.5</td><td>10<sup>-2</sup></td><td>675</td><td>-0.35</td></tr><tr><td>Best-fit grid model if α<sub>N<sub>H</sub>inner</sub> ≤ +1.0</td><td>10<sup>23</sup></td><td>0.0</td><td>-2.0</td><td>10<sup>-2</sup></td><td>650</td><td>-0.35</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Fig. 14.** Vertically integrated optical depth of the dust at 4.7 µm from the HD 139614 model of Matter et al. (2016), which best describes the SED and VLTI IR-interferometry data. The dust is optically thin at 4.7 µm in the inner 6 AU of the disk. At R > 6 AU the dust is optically thick except in the uppermost layers of the disk.</figcaption></figure>

with the upper limit derived from the CO ro-vibrational line-profile modeling.

## 4.6. Uncertainties, limitations and robustness tests

Optical depth effects explain the higher probabilities of models with $\alpha_{N_{H}inner}$ of +2 or +3. However, such high $\alpha_{N_{H}inner}$ are hard to justify physically<fnref n="14" />. To test the robustness of the modeling conclusions, we recalculated the Bayesian probability plots for a sub-grid selecting only models with $\alpha_{N_{H}inner} \le +1$ (54 000 models, Fig. C.1). We found that a gas density drop and gas surface density at R < 6 AU that is flat are still needed to explain the observations. In Fig. C.2 and Table 5, we display the characteristics of the best-fit model of the $\alpha_{N_{H}inner} \le +1$ sub-grid. The model has $\delta_{gas} = 10^{-2}$ and has a flat ($\alpha_{N_{H}inner} = 0.0$) surface density at R < 6 AU<fnref n="15" />.

C¹⁸O emission is detected at lower S/N than ¹²CO and ¹³CO emission. To further test the reliability of the gas density drop, we recalculated the Bayesian probability diagrams only taking into account the ¹²CO and ¹³CO emission (Fig. C.3). We retrieved the same results: the models with the highest probabilities are those with a gas density drop of at least a factor 100 in the inner 6 AU and that have a surface density at R < 6 AU that is flat or increasing with radius.

Our models assume constant ¹²CO/H₂, ¹²CO/¹³CO and ¹²CO/C¹⁸O ratios with radius. The CO ro-vibrational lines trace warm CO, therefore freeze-out onto dust surfaces and fractionation reactions are not a concern. Photodissociation by UV photons might be relevant, because the self-shielding of ¹³CO and C¹⁸O requires higher column densities than for ¹²CO. Models of disks including selective photodissociation (e.g., Miotello et al. 2014) showed that the gas masses in the outer disk can be underestimated by up to an order of magnitude<fnref n="16" />. However, an underestimation of the column density that is due to selective photodissociation by a factor of ten in the inner disk would not change the conclusion that a surface density drop is required to describe the CO ro-vibrational data. Moreover, van der Marel et al. (2016) recently modeled ALMA CO sub-mm rotational emission from the gas in transition disks with large dust cavities. They found that selective-photodissociation does not significantly affect the CO isotopologs rotational emission from gas inside the dust cavity.

We have used a single vertical temperature for each radius, but disks have a vertical gradient of temperature. At R < 6 AU, the dust (see Fig. 14) and the ¹³CO and C¹⁸O lines (see Fig. 13) are optically thin<fnref n="17" />. Therefore, at R < 6 AU, the ¹³CO and C¹⁸O transitions trace the whole vertical column of gas. Although in a disk the temperature increases from the mid-plane to the surface, in the inner 6 AU, the ¹³CO and C¹⁸O lines are not dominated by the hottest gas (T > 500–1000 K) located in the upper most layers near the surface because the amount of gas in those upper regions is very small. The ¹³CO and C¹⁸O ro-vibrational lines are emitted lower down, in the region where the CO gas is the densest and where the gas temperature in the vertical direction varies by a few 10 K at most down to the mid-plane<fnref n="18" />. The ¹²CO emission is optically thick at R < 6 AU, which means that on average, it is emitted higher up in the disk. The dominant emitting regions for ¹²CO, ¹³CO and C¹⁸O ro-vibrational emission have differences in vertical height but, in fact, they overlap in the vertical direction. The temperature profile in our simple 1D models should be understood as a representative *average* vertical temperature<fnref n="19" />. The column of C¹⁸O (thus $N_{H}$) could be underestimated because the 1D models use a higher average temperature.

<footnote n="14">¹⁴ X-ray photoevaporation models (e.g., Owen et al. 2012; Mordasini et al. 2012) evolve to a flat surface density. Planet-disk interaction models (e.g., Crida & Morbidelli 2007) can give a positive density gradient with radius in the inner disk, but only in the innermost radii.</footnote>

<footnote n="15">¹⁵ We note, however, that models with flat surface density profiles produce a poorer fit to the observations than models with positive $\alpha_{N_{H}inner}$ because of the stronger ¹²CO high-velocity wings (Fig. 12).</footnote>

<footnote n="16">¹⁶ For an application of selective photodissociation to the Solar Nebula see Lyons & Young (2005).</footnote>

<footnote n="17">¹⁷ An additional argument for an optically thin inner disk comes from the detailed radiative transfer modeling of CO ro-vibrational emission in the Herbig F4V pre-transition disk HD 135344B in Carmona et al. (2014). We found that, to have CO ro-vibrational emission from inside the dust cavity, the dust in the inner-most disk should be optically thin.</footnote>

<footnote n="18">¹⁸ See Fig. A.4 in Carmona et al. (2014) for the temperature distribution of the gas inside an optically thin dust cavity of a transition disk.</footnote>

<footnote n="19">¹⁹ The best-fit grid model has a radial temperature that is a compromise between the ¹²CO that requires a higher temperature (emitted higher up) and ¹³CO and C¹⁸O emission that require lower temperatures (emitted farther down). The best-fit model combining ¹²CO, ¹³CO, and C¹⁸O has a gas $T_0$ (R = 1 AU) = 675 K. This temperature is in between the temperature of the best-fit grid model for ¹²CO alone ($T_0$ (R=1 AU) = 725 K) and the temperature of the best-fit grid model for C¹⁸O emission alone ($T_0$ (R = 1 AU) = 625 K).</footnote>
<figure><img src="image_24.png" /><figcaption>**Fig. 15.** $^{12}$CO P(9) line-profiles predicted for different gap sizes and diverse column densities inside the gap, for the best-fit grid model (see Fig. 13), compared with composite $^{12}$CO line-profile. Observed and modeled line-profiles are normalized such that the continuum is at zero and the line peak at 1. The color code in the surface density panel corresponds to the color in the line-profile plot. The solution without gap is displayed in black. The temperature profile is kept constant in all the models. Error bars are $1\sigma$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 16.** Predicted $^{12}$CO P(9) line-profile for a disk extending from 0.1 to 1.0 AU with a flat surface density $N_H = 5 \times 10^{19} \text{ cm}^{-2}$ and the extrapolated temperature profile from the best-fit grid model. Error bars are $3\sigma$. The horizontal dashed line is the $5\sigma$ limit.</figcaption></figure>

However, the column density in the inner 6 AU should be lower than $N_H = 10^{21} \text{ cm}^{-2}$. Higher $N_H$, even with $T_0 (R = 1 \text{ AU}) = 575 \text{ K}$ (100 K lower than the best-fit grid model), would generate $^{13}$CO and $C^{18}$O lines with high-velocity wings which would be too strong to be compatible with the CRIRES spectrum.

We assumed a smooth temperature profile with radius, without bumps or discontinuities. However, as the dust density is much lower at $R < 6 \text{ AU}$, the temperature at the inner rim of the outer disk might be higher, as seen in thermochemical disks models (Thi et al. 2013; Carmona et al. 2014; Woitke et al. 2016; Hein Bertelsen 2015). The question in HD 139614 is which fraction of the column of the emitting CO around 6 AU is at higher temperatures. The smoothness and width of the $^{12}$CO line-profile, and the fact that the average temperatures of $^{13}$CO and $C^{18}$O lines (380 and 350 K respectively) are similar to that of the power-law temperature at $6 < R < 10 \text{ AU}$ of the

best-fit grid model (350–300 K, where most of the $^{13}$CO and $C^{18}$O flux is produced) suggest that the column of CO at temperatures much higher than 400 K in the inner rim of the outer disk AU should be small. We conclude that a smooth-temperature power-law describes the temperature of the largest column of gas emitting the CO ro-vibrational lines.

We assumed that the gas density drop occurs at the same radius as the dust density drop. But the gas density drop does not have to occur at 6 AU. We have tested a sub-grid of models (28 800 models) in which we varied the radius of the gas density drop between 4.0 and 6.0 AU (see Fig. C.4)$^{20}$. The grid shows that the most likely value for the gas density drop is 6.0 AU. Models with a gas density drop down to 5.0 AU can also describe the data but are less likely$^{21}$.

The model grid was calculated assuming $R_{\text{in}} = 1 \text{ AU}$ because the power-law intensity model indicated that the radius of the maximum intensity is close to 1 AU (Sect. 4.1). We have tested models with gas down to 0.1 AU extending the power-law temperature and density profile. The fit was satisfactory for surface density exponents at $R < 6 \text{ AU}$ ($\alpha_{N_{\text{H inner}}}$) between +1 and +3. Models with $\alpha_{N_{\text{H inner}}}$ smaller than +1 gave too strong line wings for the $^{12}$CO emission. The best-fit grid model, which has $\alpha = +2.0$ surface density exponent, gave a good fit when we extended it down to 0.1 AU (Fig. C.5). This model is consistent with the upper limits on the surface density at $R < 1 \text{ AU}$ derived in Sect. 4.5.

In the grid of models we did not fit the $v = 2 \rightarrow 1^{12}\text{CO}$ lines. We checked the predicted $v = 2 \rightarrow 1^{12}\text{CO P(3)}$ and P(4) lines for the best-fit grid model. We found that the model is able to reproduce the observed FWHM of the $v = 2 \rightarrow 1$ lines, but the

<footnote n="20">20 That dust is seen at $R \sim 6 \text{ AU}$ and that the dust requires a scale height to fit the near-IR data (SED and visibilities) implies that gas should be present at the location of the dust. For this reason larger radii than 6 AU for the gas density drop were not considered in the modeling.</footnote>

<footnote n="21">21 The change in the gas-density drop radius can be compensated for by changes in $\delta_{\text{gas}}$ or the surface density or temperature exponent.</footnote>
<figure><img src="image_26.png" /></figure>

**Fig. 17.** Gas and dust surface density. In blue we show the dust surface density derived from modeling the SED + near and mid-IR observations, taken from Matter et al. (2016). In red we plot the gas surface density of the best-fit model of the whole grid. In orange we show the gas surface density of the best-fit grid model restricted to only models with $\alpha_{N_{H\text{inner}}} \le +1.0$. Note that the value of $\delta_{\text{gas}}$ depends on the gas-to-dust ratio assumed for the outer disk. This plot assumes a distance of 140 pc for HD 139614. The new *Gaia* distance of 131 ± 5 pc implies that the location of the gas and dust density drop is 0.4 AU closer in. However, this difference is within the uncertainties of the modeling of the data.

model has weaker $v = 2 \rightarrow 1$ line fluxes than the observations. We explored models around the best-fit grid solution. We found that a model with an $N_H$ at 6 AU three times larger ($N_H$ ($R=6$ AU) = $3 \times 10^{23}$ cm$^{-2}$), the same density profile at $R < 6$ AU (thus $\delta_{\text{gas}} = 3.3 \times 10^{-3}$), and the same temperature profile described the $v = 2 \rightarrow 1$ lines well while having a good fit to the $v = 1 \rightarrow 0$ lines. The model indicates that the $v = 2 \rightarrow 1$ lines are dominated by the contribution at $6 < R < 10$ AU. We present the predicted emission lines of this model in Fig. C.6.

Our model grids were calculated before the *Gaia* distance release, and therefore assumed a distance of 140 pc. The new *Gaia* distance of 131 ± 5 pc translates into a location of gas (and dust) density drop 0.4 AU closer in. This difference is within the uncertainties of the data modeling, and therefore the conclusions we reach are not affected.

# 5. Discussion

The main results of the observations and modeling (see Table 3 for a modeling overview) of the CO ro-vibrational emission lines in the HD 139614 disk are as follows.

1. CO ro-vibrational emission extends from 1 AU to 15 AU in HD 139614, which means that the dust gap observed in IR-interferometry data contains molecular gas.

2. C$^{18}$O lines are a few km s$^{-1}$ narrower than $^{12}$CO lines.

3. $^{13}$CO and C$^{18}$O emission are 50–100 K colder than $^{12}$CO emission.

4. The observed CO emission is very likely thermally excited.

5. A drop of $10^{-2}$ to $10^{-4}$ in the gas column density at $R < 5-6$ AU is required to simultaneously reproduce the line-profiles and rotational diagrams of the three CO isotopologs.

6. The gas surface density $N_H$ at $1 < R < 6$ AU ranges between $3 \times 10^{19}$ and $10^{21}$ cm$^{-2}$ and has a distribution that most likely is flat or that increases with radius.

7. The $5\sigma$ upper limit on the CO column density $N_{\text{CO}}$ at $R \le 1$ is $5 \times 10^{15}$ cm$^{-2}$, which corresponds to a gas column density $N_H < 5 \times 10^{19}$ cm$^{-2}$ if standard abundances are assumed.

8. Our data does not show evidence of a gap (devoid of gas) in the gas distribution. The width of any possible gap is constrained to be smaller than 2 AU.

## 5.1. Gas vs. dust surface density distribution

In Fig. 17 we display the dust surface density (in blue) derived by Matter et al. (2016) from ESO/VLTI near and mid-IR interferometry observations. In the same figure we illustrate the gas surface density derived from the CO ro-vibrational emission. In red we plot the best-fit grid model, in orange the best-fit grid model if $\alpha_{N_{H\text{inner}}} \le 1.0$.

We find that the $\alpha_{\text{inner}} = +0.6$ deduced for the dust is consistent with the range of $\alpha_{N_{H\text{inner}}}$ of the gas observations. The near-IR continuum and the CO ro-vibrational observations suggest gas-to-dust mass ratios for the inner disk at $R < 2.5$ AU ranging from a few up to 100, depending on the gas surface density exponent at $R < 6$ AU. Concerning the depletion levels of the gas and the dust ($\delta_{\text{gas}}$ and $\delta_{\text{dust}}$), if the gas-to-dust mass ratio is 100 in the outer disk, then the level of depletion for the gas and the dust would be similar. However, it is likely that given the age of HD 139614 (9 Myr, Alecian et al. 2013) and the weak [O I] 63 µm line flux ($4.5 \times 10^{-17}$ W m$^{-2}$, among the weakest of the whole Herbig Ae sample of Meeus et al. 2012), that the
gas-to-dust mass ratio in the outer disk is below one hundred<fnref n="22" />. In this case, the density drop in the gas would be less deep than the density drop in the dust.

## 5.2. CO ro-vibrational emission in HD 139614 and other protoplanetary disks

The integrated ¹²CO line fluxes of HD 139614 are on the same order ($10^{-14}$ erg s⁻¹ cm⁻²) as previously observed primordial and transition disks (e.g., Najita et al. 2003; Salyk et al. 2011; Brown et al. 2013; Banzatti & Pontoppidan 2015). The widths of the ¹²CO line-profiles of HD 139614 (and other transition disks) are, however, significantly narrower and lack the high-velocity wings (emission at $v > 15$ km s⁻¹) that are observed in primordial disks<fnref n="23" />. This difference arises because in primordial disks the ¹²CO emission is dominated by gas inside 1 AU, while in HD 139614 (and most transition disks) the CO ro-vibrational emission is dominated by gas at $R > 1$ AU.

Given that ¹²CO ro-vibrational emission becomes optically thick at relatively low column densities ($N_{CO} \sim 10^{15}$ cm⁻²), the lack of high-velocity wings in the ¹²CO line-profiles of HD 139614 (and other transition disks) directly indicates that at $R < 1$ AU the column density of gas is much lower than in primordial disks. The $5\sigma$ upper limit of $N_{H} \sim 5 \times 10^{19}$ cm⁻² on the column of gas at $R < 1$ AU (assuming standard abundances) derived for HD 139614 is significantly lower than the $N_{H} = 10^{24}-10^{26}$ cm⁻² typical of primordial disks.

In HD 139614, the ¹³CO emission is narrower and colder than the ¹²CO emission. Differences between CO isotopologs line widths and temperatures have been reported in a number of primordial and transition disks (e.g., Brown et al. 2013). However, the ¹²CO and ¹³CO line width difference in primordial disks is on the order of tens of km s⁻¹, while in HD 139614 (and other transition disks) it is only a few km s⁻¹. This further suggests that there is a different inner disk gas structure between transition disks and primordial disks (as already pointed out by previous authors, e.g., Brown et al. 2013; Banzatti & Pontoppidan 2015).

The $v = 1 \rightarrow 0$ CO and $v = 2 \rightarrow 1$ CO emission in HD 139614 can be described by thermal emission ($T_{ex} = T_{rot} = T_{vib}$). The average CO temperature of 460 K ($\log T = 2.7$) and the CO inner radius of 1 AU locate HD 139614 in the left side of the UV pumping regime in the recently proposed $T/R$ diagram of Banzatti & Pontoppidan (2015). According to the $T/R$ diagram, HD 139614 belongs to disk category 2: objects with partly devoid disk's gaps. HD 139614 appears in the $T/R$ diagram as an intermediate case between category 1 disks that are primordial, and category 3 that are transition disks that have $v = 2 \rightarrow 1$ CO ro-vibrational emission dominated by UV pumping.

In summary, the properties of CO ro-vibrational emission are clearly different between HD 139614 and primordial disks. The data and models show that although there is molecular gas inside the inner 6 AU of HD 139614, the gas distribution and gas mass is different from that of a primordial disk at $R < 6$ AU.

## 5.3. Comparison to other transition disks with quantified gas surface densities inside the dust cavity

HD 139614 adds to the growing number of transition disks with quantified gas density drops inside the dust cavity. In Table C.1 we provide a summary of the properties of these sources<fnref n="24" />. The largest portion of transition disks shows stronger depletions in the dust than in the gas ($\delta_{dust} < \delta_{gas}$) inside the dust cavity. Some disks have a similar depletion level (e.g., HD 139614), and only one source, J1604-2130, has a depletion level higher in the gas than the dust. Altogether, this suggests that the dust depletes faster than the gas in the inner disk. This behavior might arise because it is hard to stop viscous accretion of gas through the disk unless the mass of the (outer) disk is very low. This seems to be the case for transition disks around low-mass/solar-type stars (Sicilia-Aguilar et al. 2015). If trapping of mm-dust grains occurs at the inner rim of the outer disk, then a mm-dust depletion higher than the gas depletion would be expected inside the dust cavity (e.g., Pinilla et al. 2016).

Table C.1 shows that in a large fraction of objects, the gas density drop radius is smaller than the sub-mm dust density drop radius. This is consistent with the detection of micron-sized particles at radii smaller than the sub-mm dust cavity radius (e.g., Muto et al. 2012; Garufi et al. 2013; Follette et al. 2013; Pinilla et al. 2016) because small dust grains are expected to be coupled to the gas. Furthermore, this also provides observational support for the scenario of dust trapping by pressure bumps due to the presence of (unseen) planets inside the dust cavity (e.g., Rice et al. 2006; Paardekooper & Mellema 2006; Fouchet et al. 2010; Pinilla et al. 2012, 2015b; Gonzalez et al. 2015). In the case of HD 139614, the dust and gas cavity radii are similar, but it is the only object in the sample with a dust cavity radius smaller than 10 AU and the only object where the dust cavity radius is measured in the IR. Future sub-mm observations of HD 139614, will enable us to test whether its sub-mm dust cavity radius is larger than the gas cavity radius.

## 5.4. Origin of the dust and gas density distribution in HD 139614

Photoevaporation has been suggested as a potential mechanism for the origin of the dust cavities in transition disks (see recent reviews by Alexander et al. 2014; Owen 2016). The key factor in this scenario is when the mass accretion rate of the disk becomes lower than the photoevaporation mass-loss rate as a result of to the radiation of the central star. The accretion rate of $10^{-8} M_{\odot}/yr$ in HD 139614 (Garcia Lopez et al. 2006) is at the high end of the mass-loss rates ($10^{-8}-10^{-10} M_{\odot}/yr$) in which photoevaporation starts to become relevant.

The inner radius of the dust gap of ~3 AU is consistent with the critical radius of the gap that is expected to be opened by EUV photoevaporation for a $1.7 M_{\odot}$ star ($R_{c,EUV} \approx 1.8 \times M_{*}/M_{\odot}$ AU, e.g., Alexander et al. 2014). However, the very fact that we detect molecular gas inside 3 AU, the lack of a gap in the gas and the accretion rate of the source do not favor the EUV photoevaporation scenario, unless we are seeing the gap just beginning to form through photoevaporation. This has a very low

<footnote n="22">22 In the detailed model of the transition disk HD 135344B (which has a dust cavity of 30 AU) by Carmona et al. (2014) we found that to reproduce its weak [O I] line ($4.7 \times 10^{-17}$ W m⁻², Meeus et al. 2012) the gas-to-dust mass ratio in the outer disk needed to be much lower than 100 with a best value below 10.</footnote>

<footnote n="23">23 The ¹²CO ro-vibrational emission from primordial disks is two to five times broader than that of transition disks (Brown et al. 2013).</footnote>

<footnote n="24">24 We note that the number of transition disks with detections of gas inside the dust cavities from CO ro-vibrational and CO rotational emission is much larger (e.g., Pontoppidan et al. 2011; Salyk et al. 2011; Brown et al. 2013; Banzatti & Pontoppidan 2015). We discuss here only the sources with published gas and dust surface density profiles.</footnote>
probability given the EUV photoevaporation timescales ($10^5$ yr). Although EUV photoevaporation is probably not the dominant mechanism for the formation of the dust gap and the gas density drop, it cannot be completely excluded because of the age of HD 139614 (~9 Myr).

X-ray photoevaporation has been suggested to be relevant for accreting transition disks because of its high mass-loss rates (e.g., Owen et al. 2011, 2012). Herbig Ae stars are, however, generally weak X-ray emitters (e.g., Stelzer et al. 2006, 2009). HD 139614 has been observed with XMM-Newton in the context of a large program investigating young stars and their protoplanetary disks (PI. M. Güdel). Details of the observations and data reduction will be described in a future publication of that survey. A first analysis of the XMM data indicate an unabsorbed X-ray flux of $5.26 \times 10^{-14}$ erg s$^{-1}$ cm$^{-2}$, which translates into an X-ray luminosity of $1.2 \times 10^{29}$ erg s$^{-1}$. Using this X-ray luminosity and the scaling relation for X-ray photoevaporation in Alexander et al. (2014), we obtain a mass-loss rate of $5.6 \times 10^{-10}$ M$_{\odot}$/yr. Since this value is much lower than the accretion rate, it is likely that X-ray photoevaporation (alone) is not the dominant mechanism responsible for the dust gap and gas density drop in the inner disk of HD 139614. Furthermore, in the X-ray photoevaporation model, a gap in the gas of 5–6 AU width would be expected for a gas density drop $\delta_{\text{gas}}$ ranging from $10^{-2}$ to $10^{-4}$ (e.g., see Fig. 9 in Owen et al. 2011). This X-ray photoevaporation-induced gap is larger than 2 AU upper limit we derive from our CRIRES CO observations$^{25}$. Finally, the lack of [O I] 6300 Å emission (Acke et al. 2005) as a tracer of photoevaporation winds (e.g., Font et al. 2004; Ercolano & Owen 2010; Gorti et al. 2011; Baldovin-Saavedra et al. 2012; Rigliaco et al. 2013) does not support photoevaporation as the main gas-depletion process taking place in HD 139614 (see also Sicilia-Aguilar et al. 2010, 2013, 2015 who argued that accreting transition disks are probably not caused by photoevaporation).

The interaction of a giant planet with its parent disk causes dramatic changes in the distributions of the gas and the dust in the disk, which alters the planet's orbital properties (planetary migration). A recent review is provided in Baruteau et al. (2014). For this discussion, we highlight a few key aspects: (1) a giant planet is expected to open a gap in the gas with a width of a few AU$^{26}$; (2) the gas surface density profile inside the planet orbit could have a variety of behaviors (be lower than in the outer disk or be radially decreasing, flat, or increasing); (3) the gas at radii close to the planet orbit can have velocities that deviate significantly from the Keplerian speed; (4) the presence of giant gap-opening planets will generate pressure bumps in the disk that will trap dust particles (for instance at radii immediately outside of the planet radius).

The width of the gap in the gas (and small dust grains) opened by a planet scales with the Hill radius ($R_H$) of the planet. The Hill radius of a planet of mass $M_p$ at a distance $r_p$ from a star of mass $M_*$ is defined by

$$R_H = r_p \left( \frac{M_p}{3M_*} \right)^{1/3}. \qquad (9)$$

The gap opened in the gas by a planet generally does not exceed five times the Hill radius (e.g., Dodson-Robinson & Salyk 2011). The mass of HD 139614 is $1.7 M_{\odot}$. If we assume that the planet is located at 4.5 AU, according to Eq. (9), planets more massive than $3.7 M_J$ would be expected to open gaps larger than 2 AU in the disk. Matter et al. (2016) presented hydrodynamical simulations adapted to HD 139614 with the objective of exploring a single giant planet as a cause of the dust gap observed in the IR interferometry observations. We refer to that paper for the details of the modeling. The results of the hydrodynamical simulations after 100 000 orbits (1 Myr) are the following: (1) planets of 1.7, 3, and $6.8 M_J$ located at 4.5 AU produce a gas gap of width 2, 3, and 4 AU respectively; (2) the surface density profile in the inner disk ($R < 3$ AU) changes as a result of a $3 M_J$ planet from an initial $R^{-1}$ to $R^{+0.6}$ profile, and features a reduction factor ~10 relative to the outer disk. In the context of the planet-disk simulations, if a planet is responsible for the dust gap and the gas density drop it should have a mass lower than $2 M_J^{27}$.

The exponent of +0.6 observed in the inner disk in the hydro-simulations is compatible with the CO observations. The gas density drop in the hydro-simulations is, however, at least a factor 10 weaker than the $\delta_{\text{gas}}$ suggested from the CO-rovibrational data. But hydrosimulations were run for 1 Myr while HD 139614 has an age ~9 Myr. It is indeed possible that the inner disk has lost a significant fraction of its gas mass due to accretion onto the star and also photoevaporation, a process which is not included in the Matter et al. hydro-simulations$^{28}$.

Moreover, a $1-2 M_J$ giant planet could be responsible for the observed 3.5 AU-wide dust gap, while having a gas gap smaller than 2 AU. A planet can generate pressure bumps in the inner and outer edge of the gap, which could trap particles (e.g., Pinilla et al. 2012). The location of these dust traps are at a radius smaller (for the inner disk) and larger (for the outer disk) than the inner and outer edge of the gas gap opened by the planet (e.g., Pinilla et al. 2016). The gas-dust interaction could thus produce gaps of different width for the dust and for the gas as observed in HD 139614.

Regály et al. (2014) have computed the CO ro-vibrational lines profiles expected for a disk with an embedded giant planet for stars with different masses and planets at different separations. For a $2 M_{\odot}$ star harboring a $10 M_J$ planet at 3 and 5 AU separations, they have predicted asymmetric CO ro-vibrational lines with distortions on the order of 10% with respect to the symmetric line-profiles. The composite $^{12}$CO line-profile is symmetric and distortions are not detected at the $1\sigma$ level. This indicates that the mass of the planet, if present, must be lower than $10 M_J$, which is consistent with the upper limit of $2 M_J$ from the lack of a gap in the gas distribution.

In addition to dynamical interaction with embedded planets and photoevaporation, various mechanisms to trap dust particles have been proposed to explain the gaps, lopsided shapes, and ringed structures observed recently by ALMA in disks (e.g., Pérez et al. 2014; van der Marel et al. 2015b; Brogan et al. 2015; Andrews et al. 2016). Scenarios include a magnetical origin such as radial pressure variations due to MHD turbulence (zonal flows Johansen et al. 2009), radial variations in the disk resistivity (e.g., Flock et al. 2015; Lyra et al. 2015), vortex formation due to instabilities at the edge of the dead-zone (e.g., Regály et al. 2013; Faure et al. 2015). Or a chemical origin such

<footnote n="25">25 We note, however, that the X-ray photoevaporation models described here were mostly developed for 1 M$_{\odot}$ stars and may not be directly applicable to A-type stars.</footnote>

<footnote n="26">26 We note that the gap width and depth depend on the planet mass and disk physical properties such as the temperature and turbulent viscosity near the planet orbit.</footnote>

<footnote n="27">27 We note that the planet mass depends on assumptions such as the disk viscosity. Models of Matter et al. (2016) assumed $\alpha = 0.006$.</footnote>

<footnote n="28">28 For a study of the interplay of photoevaporation and planet formation, see, for example, Rosotti et al. (2013).</footnote>
as the changes in opacity that are due to migrating solids reach-
ing the condensation fronts of volatiles (e.g., Cuzzi & Zahnle
2004; Brauer et al. 2008; Banzatti et al. 2015; Okuzumi et al.
2016). While these scenarios could in part explain the dust fea-
tures discovered in disks, it is not clear whether they are able
to explain the gas density drop that the CO ro-vibrational emis-
sion reveals in HD 139614 (and other transition disks). The low
gas column-density detected inside the dust cavities of transi-
tion disks, combined with the low optical depth of the dust can,
however, impact the ionization structure of the disk and have an
effect on MHD phenomena.

5.5. A dust trap in the inner disk?

The possibility of an increasing gas surface density profile in
the inner 6 AU has interesting consequences from the point of
view of gas and dust evolution. If the gas surface density pro-
file within 6 au of HD 139614 is indeed the result of a giant
planet of mass between 1 and 2 M<sub>J</sub> at 4.5 AU, then pressure
bumps would be present at the two edges of the planetary gap.
Dust can accumulate, grow, and be trapped in these pressure
bumps (e.g., Fouchet et al. 2010; Pinilla et al. 2012). The trap-
ping of the particles depends on how well they are coupled to
the gas and therefore, it depends on their size and local gas sur-
face density. Quantitatively, particles with size $a_{\text{opt}} = 2\Sigma_{\text{gas}}/\pi\rho_{\text{d}}$
drift the fastest toward the regions of high pressure and are the
most efficient particles to trap (e.g., Fouchet et al. 2010). Here
$\rho_{\text{d}}$ is the internal density of the dust grains, typically $1-3 \text{ g cm}^{-3}$,
$\Sigma_{\text{g}}$ is the gas surface density. Using $N_{\text{H}} = 10^{19}-10^{21} \text{ cm}^{-2}$
$(2.4\times 10^{-5}-2.4\times 10^{-3} \text{ g cm}^{-2})^{29}$ for the inner disk of HD 139614,
we obtain that sub-micron and micron-sized grains are trapped
at the inner gap edge. This can lead to an increasing dust sur-
face density with radius in the inner disk, as the analysis of IR
interferometry observations suggested (Matter et al. 2016).

5.6. Gas surface density in the inner disk and accretion rate
onto the star

Many of the bright transition disks have stellar accretion rates ($\dot{M}_{\star}$) between $10^{-9}$ and $10^{-8} M_{\odot} \text{yr}^{-1}$, very similar to accretion rates of classical T-Tauri stars with primordial gas-rich disks (e.g., Manara et al. 2014). However, the surface density of the gas in the cavities of transition disks typically ranges from $10^{-3}$ to $1 \text{ g cm}^{-2}$ (see Table 6), which is several orders of magnitude lower than at similar locations in primordial disks. Also quite surprisingly, transition disks with similar $\dot{M}_{\star}$ can have very different gas surface densities in their inner regions. This is the case of IRS 48 and RXJ1615, both of which have $\dot{M}_{\star} \sim$ a few $\times 10^{-9} M_{\odot} \text{yr}^{-1}$, but estimated gas surface densities inside their cavities that differ by roughly two orders of magnitude, even though the cavities are of similar size in the gas (van der Marel et al. 2015b, 2016). These all seem to be counter intuitive facts when considering that in a steady-state model of a protoplanetary disk, the stellar accretion rate and the disk accretion rate should take similar values throughout the disk.

The disk accretion rate, $\dot{M}_{\text{disk}}$, is related to the surface den-
sity $\Sigma$ and radial velocity $v_R$ of the gas at radius $R$ through
$\dot{M}_{\text{disk}} = -2\pi R v_R(R) \Sigma(R)$. If we now assume that $\dot{M}_{\text{disk}} \sim \dot{M}_{\star}$,
then at a radius of 1 AU and for an accretion rate of $10^{-8} M_{\odot} \text{yr}^{-1}$,
we find that $|v_R|$ should extend from about 70 km s$^{-1}$ down to
$0.07 \text{ km s}^{-1}$ for $\Sigma$ in the range $[10^{-3}-1] \text{ g cm}^{-2}$. At 1 AU, the

Keplerian velocity ($v_K$) is about 30 km s$^{-1}$, which means that for
the transition disks with the lowest densities inside the cavity, the
gas radial velocity can be comparable to the Keplerian velocity.
The above range of radial velocities is at odds with typical val-
ues expected in classical viscous disk models. These models tell
us that $|v_R| \sim \alpha h^2 v_K$, were $h$ is the disk aspect ratio and $\alpha$ is the
dimensionless turbulent viscosity of the disk. If the disk mag-
netic field is able to sustain vigorous turbulence in the cavities
of transition disks, then $|v_R|$ can reach at most $\sim 10^{-5} v_K$. This is
at least two orders of magnitude lower than the above range of
radial velocities obtained when assuming $\dot{M}_{\text{disk}} \sim \dot{M}_{\star}$.

Another way to visualize the problem is to think in terms of
the gas column density expected for the accretion rates reported,
and the CO ro-vibrational line-profiles that such column densi-
ties would produce. For HD 139614, the $\dot{M}_{\text{star}} = 10^{-8} M_{\odot} \text{yr}^{-1}$
would suggest an $N_{\text{H}}$ of $10^{25}-10^{26} \text{ cm}^{-2}$ at $R \le 1$ AU for $\alpha$
ranging from $10^{-2}$ to $10^{-3}$. These gas columns are high enough
for the CO to self-shield against UV photodissociation and pro-
duce strong CO ro-vibrational emission from $R < 1$ AU (not seen
in our CRIRES spectrum). For example, if $N_{\text{H}} \sim 10^{25} \text{ cm}^{-2}$ at
$R < 1$ AU in a disk without dust in the cavity (best case for UV-
photodissociation) around a $T_{\text{eff}} = 10000$ K central star (i.e.,
brighter in UV than HD 139614) $N_{\text{CO}}$ would be $10^{21} \text{ cm}^{-2}$ at
$R < 1$ AU (see Fig. 8 in Bruderer 2013), a value much higher
than the $5\sigma$ limit of $N_{\text{CO}} = 5 \times 10^{15} \text{ cm}^{-2}$ from our CRIRES data.
According to the Bruderer (2013) models (for a $T_{\text{eff}} = 10000$ K),
to have $N_{\text{CO}} < 10^{15} \text{ cm}^{-2}$, $N_{\text{H}}$ needs to below $10^{22} \text{ cm}^{-2}$
at $R < 1$ AU for an inner disk without dust. For HD 139614, which
has a $T_{\text{eff}} \sim 7800$ K (and lower UV, accordingly) and some dust
in the inner disk$^{30}$. $N_{\text{H}}$ should be even lower to enable the effi-
cient photodissociation of CO.

There are probably two ways to solve these problems. The
first is to concede that the inner parts of transition disks do not
have to be in a steady state with $\dot{M}_{\text{disk}} \sim \dot{M}_{\star}$. This implies that
measured accretion rates should not be used as direct proxies to
derive the amount of gas left in the cavities of transition disks.
The second way is to conceive that the radial flow of gas inside
the cavities of transition disks is not necessarily driven by tur-
bulent accretion, but by interactions with one or more massive
companion(s) inside the cavity. HD 142527 may be a pioneering
example of the latter possibility. Casassus et al. (2015) showed
that a highly inclined stellar companion to HD 142527 can gen-
erate fast radial flows inside the cavity of its transition disk.
Pontoppidan et al. (2011) detected in HD 142527 a $^{12}\text{CO}$ ro-
vibrational spectrum and spectro-astrometry signature display-
ing asymmetries, indicating non-Keplerian contributions to the
emission. Additional resolved observations of the gas kinemat-
ics inside the cavities of transition disks will help assess the oc-
currence of this second scenario. In our case of HD 139614, the
$^{12}\text{CO}$ composite line-profile is symmetric, smooth, and consis-
tent with Keplerian motion, which suggests the first scenario,
namely that $\dot{M}_{\text{disk}}$ (1 AU) is most likely different from $\dot{M}_{\star}$.

6. Summary and conclusions

We have obtained VLT/CRIRES high-resolution spectra ($R \sim 90\,000$) of CO ro-vibrational emission at 4.7 µm in HD 139614, an accreting ($10^{-8} M_{\odot} \text{yr}^{-1}$) Herbig Ae star with a (pre-) transition disk that is characterized by a dust gap between 2.3 and 6 AU and a dust density drop $\delta_{\text{dust}}$ of $10^{-4}$ at $R < 6$ (Matter et al. 2016). We have detected $v = 1 \rightarrow 0^{12}\text{CO}, ^{13}\text{CO}, \text{C}^{18}\text{O}, \text{C}^{17}\text{O}$,

<footnote n="29"><sup>29</sup> Here the hydrogen mass fraction of gas with solar composition is assumed (0.7) to calculate the gas surface density in g cm<sup>-2</sup>.</footnote>

<footnote n="30">³⁰ The dust at *R* < 2.5 AU is optically thin at 4.7 µm, but in the UV it has higher opacity, thus helping to shield CO.</footnote>
and $v = 2 \rightarrow 1$ $^{12}$CO ro-vibrational emission. The lines observed are consistent with disk emission and thermal excitation. We find the following:

1. The $v = 1 \rightarrow 0$ $^{12}$CO spectrum indicates that there is gas from 1 AU up to 15 AU, and that there is no gap in the gas distribution. If a gap is present in the gas (i.e., a region devoid of gas) then it should have a width smaller than 2 AU.

2. The spectra of $^{13}$CO and $C^{18}$O $v = 1 \rightarrow 0$ emission are on average colder and emitted farther out in the disk ($R > 6$ AU) than the $^{12}$CO $v = 1 \rightarrow 0$ emission. Keplerian flat-disk models clearly show that a drop in the gas density $\delta_{\text{gas}}$ of a factor of at least 100 at $R < 5-6$ AU is needed to describe simultaneously the line-profiles and rotational diagrams of the three CO isotopologs. Models without a gas density drop produce $C^{18}$O and $^{13}$CO lines that are too wide and warm to be compatible with the data. If the gas-to-dust mass ratio is equal to 100 in the outer disk, the gas depletion factor $\delta_{\text{gas}}$ could be as high as $10^{-4}$. Moreover, we find that the gas surface density profile in the inner 6 AU of the disk is flat or increases with radius.

The presence of molecular gas inside 6 AU and the weak X-ray luminosity do not favor photoevaporation as the main mechanism responsible for the inner disk structure of HD 139614. The gas density drop, a flat or increasing gas surface density profile at $R < 6$ AU, combined with the non-detection of a gap in the gas wider than 2 AU, suggest the presence of a single Jovian-mass planet inside the dust gap. If a giant planet is indeed responsible for the transition disk shape of HD 139614, then its location would be at around 4 AU and its mass would be lower than two Jupiter masses. Furthermore, if a small gap in the gas (due to a planet) were to be present, a gas surface density profile that increases with radius in the inner disk might lead to a dust trap at the gap inner edge for sub-micron and micron sized grains, which could explain that the dust surface density increases with radius at $R < 2.5$ AU, as found in IR interferometry observations.

We constrained the gas column density between 1 and 6 AU to $N_{\text{H}} = 3 \times 10^{19} - 10^{21} \text{ cm}^{-2}$ ($7 \times 10^{-5} - 2.4 \times 10^{-3} \text{ g cm}^{-2}$) assuming $N_{\text{CO}} = 10^{-4} N_{\text{H}}$. We derived a $5\sigma$ upper limit on the CO column density at $R < 1$ AU $N_{\text{CO}} = 5 \times 10^{15} \text{ cm}^{-2}$, which suggests an $N_{\text{H}} < 5 \times 10^{19} \text{ cm}^{-2}$ at $R < 1$ AU. The gas surface density in the disk of HD 139614 at $R \le 1$ AU and at $1 < R < 6$ AU is significantly lower than the surface density that would be expected for the accretion rate of HD 139614, assuming a standard viscous $\alpha$-disk model. Our result, and the low gas surface densities reported in the inner disks of other transition disks, suggests that stellar accretion rates should *not* be used as direct proxies to derive the amount of gas left inside the dust cavities of transition disks. An investigation of the topology of the magnetic fields of young stars with transition disks is needed to help address the question of the differences between the accretion rate and the inner disk gas surface density.

We have discussed the ensemble of transition disks with current constraints for the gas surface density inside the dust cavity. The sample shows that, in the majority of the sources, the drop in the dust density is larger than the drop in the gas density ($\delta_{\text{dust}} < \delta_{\text{gas}}$). This suggests that dust is depleted faster than gas in the inner disk.

The number of transition disks with a complete set of multi-wavelength observations of gas and dust (ALMA, CRIRES, Herschel, Spitzer, VLTI, SEDs, HiCiAO, VLT/SPHERE, and GPI) is growing. A homogenous multi-wavelength and multi-technique modeling of gas and dust observations in transition

disks would be of great help to understand the variety of gas and dust structures that these disks have, and to study the possible links to planet formation. In that respect, it would be of great help to have spatially resolved measurements of the dust and the rotational transitions of CO isotopologs in the sub-mm for HD 139614 (for example with ALMA$^{31}$).

*Acknowledgements.* A. Carmona was partly supported by the Spanish Grant AYA 2011-26202. A. Carmona, A. Kóspal and Zs. Regály were partly supported by the Momentum grant of the MTA CSFK Lendület Disk Research Group of the Hungarian Academy of Sciences. A. Carmona and C. Pinte acknowledge funding from the European Commission's 7th Framework Program (EC-FP7) (contract PERG06-GA-2009-256513) and from Agence Nationale pour la Recherche (ANR) of France under contract ANR-2010-JCJC-0504-01. A. Carmona acknowledges financial support by the European Southern Observatory visitors program. The research leading to these results has received funding from the EC-FP7 under grant agreement no 284405. C. Eiroa is partly supported by the Spanish Grant AYA 2014-55840-P. L.A. Cieza was supported by CONICYT-FONDECYT grant number 1140109 and the Millennium Science Initiative (Chilean Ministry of Economy), through grant "Nucleus RC130007". A.C. would like to thank C. P. Folsom for comments on the manuscript. A.C. and C.B would like to thank G. Lesur and S. Casassus for discussions on protoplanetary disks.

## References

Acke, B., van den Ancker, M. E., & Dullemond, C. P. 2005, A&A, 436, 209

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

Aleanian, E., Wade, G. A., Catala, C., et al. 2013, MNRAS, 429, 1001

Alexander, R. D., & Armitage, P. J. 2007, MNRAS, 375, 500

Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, Protostars and Planets VI, 475

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40

Artymowicz, P., & Lubow, S. H. 1996, ApJ, 467, L77

Baldovin-Saavedra, C., Audard, M., Carmona, A., et al. 2012, A&A, 543, A30

Banzatti, A., & Pontoppidan, K. M. 2015, ApJ, 809, 167

Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15

Baruteau, C., Crida, A., Paardekooper, S.-J., et al. 2014, Protostars and Planets VI, 667

Bary, J. S., Weintraub, D. A., & Kastner, J. H. 2003, ApJ, 586, 1136

Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6

Birnstiel, T., Andrews, S. M., & Ercolano, B. 2012, A&A, 544, A79

Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1

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

Brogan, C. L., Pérez, L. M., et al. (ALMA Partnership) 2015, ApJ, 808, L3

Brown, J. M., Blake, G. A., Qi, C., et al. 2009, ApJ, 704, 496

Brown, J. M., Pontoppidan, K. M., van Dishoeck, E. F., et al. 2013, ApJ, 770, 94

Bruderer, S. 2013, A&A, 559, A46

Bruderer, S., van der Marel, N., van Dishoeck, E. F., & van Kempen, T. A. 2014, A&A, 562, A26

Canovas, H., Schreiber, M. R., Cáceres, C., et al. 2015, ApJ, 805, 21

Carmona, A., Pinte, C., Thi, W. F., et al. 2014, A&A, 567, A51

Casassus, S., van der Plas, G., M, S. P., et al. 2013, Nature, 493, 191

Casassus, S., Marino, S., Pérez, S., et al. 2015, ApJ, 811, 92

Chandra, S., Maheshwari, V. U., & Sharma, A. K. 1996, A&AS, 117, 557

Chiang, E., & Murray-Clay, R. 2007, Nat. Phys., 3, 604

Cieza, L. A., Mathews, G. S., Williams, J. P., et al. 2012, ApJ, 752, 75

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Crida, A., & Morbidelli, A. 2007, MNRAS, 377, 1324

Cutri, R. M., et al. 2012, VizieR Online Data Catalog: II/311

Cuzzi, J. N., & Zahnle, K. J. 2004, ApJ, 614, 490

Dent, W. R. F., Greaves, J. S., & Coulson, I. M. 2005, MNRAS, 359, 663

Dittkrist, K.-M., Mordasini, C., Klahr, H., Alibert, Y., & Henning. T. 2014, A&A, 567, A121

Dodson-Robinson, S. E., & Salyk, C. 2011, ApJ, 738, 131

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

Ercolano, B., & Owen, J. E. 2010, MNRAS, 406, 1553

Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, Protostars and Planets VI, 497

Fang, M., van Boekel, R., Wang, W., et al. 2009, A&A, 504, 461

Fang, M., Sicilia-Aguilar, A., Roccatagliata, V., et al. 2014, A&A, 570, A118

Faure, J., Fromang, S., Latter, H., & Meheut, H. 2015, A&A, 573, A132

<footnote n="31">31 At the time of writing, HD 139614 has not been observed with ALMA.</footnote>
Fedele, D., Bruderer, S., van Dishoeck, E. F., et al. 2013, A&A, 559, A77

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68

Follette, K. B., Tamura, M., Hashimoto, J., et al. 2013, ApJ, 767, 10

Follette, K. B., Grady, C. A., Swearingen, J. R., et al. 2015, ApJ, 798, 132

Folsom, C. P., Bagnulo, S., Wade, G. A., et al. 2012, MNRAS, 422, 2072

Font, A. S., McCarthy, I. G., Johnstone, D., & Ballantyne, D. R. 2004, ApJ, 607, 890

Fouchet, L., Gonzalez, J.-F., & Maddison, S. T. 2010, A&A, 518, A16

Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2

García Lopez, R., Natta, A., Testi, L., & Habart, E. 2006, A&A, 459, 837

Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105

Gonzalez, J.-F., Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015, MNRAS, 454, L36

Gorti, U., Hollenbach, D., Najita, J., & Pascucci, I. 2011, ApJ, 735, 90

Hein Bertelsen, R. 2015, Ph.D. Thesis, Groningen University

Hein Bertelsen, R. P., Kamp, I., van der Plas, G., et al. 2016, MNRAS, 458, 1466

Ireland, M. J., & Kraus, A. L. 2008, ApJ, 678, L59

Isella, A., Pérez, L. M., Carpenter, J. M., et al. 2013, ApJ, 775, 30

Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269

Johnson, H. L. 1966, ARA&A, 4, 193

Juhász, A., Bouwman, J., Henning, T., et al. 2010, ApJ, 721, 431

Kaüfl, H. U., Ballester, P., Biereichel, P., et al. 2014, SPIE, 5492

Kley, W. 1999, MNRAS, 303, 696

Kóspál, Á., Ábrahám, P., Acosta-Pulido, J. A., et al. 2012, ApJS, 201, 11

Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001

Lyons, J. R., & Young, E. D. 2005, Nature, 435, 317

Lyra, W., Turner, N. J., & McNally, C. P. 2015, A&A, 574, A10

Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18

Mariñas, N., Telesco, C. M., Fisher, R. S., & Packham, C. 2011, ApJ, 737, 57

Marsh, K. A., & Mahoney, M. J. 1992, ApJ, 395, L115

Mathews, G. S., Williams, J. P., & Ménard, F. 2012, ApJ, 753, 59

Matter, A., Labadie, L., Kreplin, A., et al. 2014, A&A, 561, A26

Matter, A., Labadie, L., Augereau, J. C., et al. 2016, A&A, 586, A11

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

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

Meeus, G., Salyk, C., Bruderer, S., et al. 2013, A&A, 559, A84

Mendigutía, I., Fairlamb, J., Montesinos, B., et al. 2014, ApJ, 790, 21

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A&A, 547, A112

Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJ, 748, L22

Najita, J., Carr, J. S., & Mathieu, R. D. 2003, ApJ, 589, 931

Okuzumi, S., Momose, M., Sirono, S.-I., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82

Owen, J. E. 2016, PASA, 33, e005

Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13

Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880

Paardekooper, S.-J., & Mellema, G. 2006, A&A, 453, 1129

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

Pérez, L. M., Isella, A., Carpenter, J. M., & Chandler, C. J. 2014, ApJ, 783, L13

Perez, S., Casassus, S., Ménard, F., et al. 2015, ApJ, 798, 85

Piétu, V., Dutrey, A., Guilloteau, S., Chapillon, E., & Pety, J. 2006, A&A, 460, L43

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114

Pinilla, P., de Boer, J., Benisty, M., et al. 2015a, A&A, 584, L4

Pinilla, P., de Juan Ovelar, M., Ataiee, S., et al. 2015b, A&A, 573, A9

Pinilla, P., Klarmann, L., Birnstiel, T., et al. 2016, A&A, 585, A35

Pinte, C., & Laibe, G. 2014, A&A, 565, A129

Pinte, C., Fouchet, L., Ménard, F., Gonzalez, J.-F., & Duchêne, G. 2007, A&A, 469, 963

Pontoppidan, K. M., Blake, G. A., van Dishoeck, E. F., et al. 2008, ApJ, 684, 1323

Pontoppidan, K. M., Blake, G. A., & Smette, A. 2011, ApJ, 733, 84

Quillen, A. C., Blackman, E. G., Frank, A., & Varnière, P. 2004, ApJ, 612, L137

Regály, Z., Juhász, A., Sándor, Z., & Dullemond, C. P. 2012, MNRAS, 419, 1701

Regály, Z., Sándor, Z., Csomós, P., & Ataiee, S. 2013, MNRAS, 433, 2626

Regály, Z., Király, S., & Kiss, L. L. 2014, ApJ, 785, L31

Rice, W. K. M., Wood, K., Armitage, P. J., Whitney, B. A., & Bjorkman, J. E. 2003, MNRAS, 342, 79

Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619

Rigliaco, E., Pascucci, I., Gorti, U., Edwards, S., & Hollenbach, D. 2013, ApJ, 772, 60

Rosotti, G. P., Ercolano, B., Owen, J. E., & Armitage, P. J. 2013, MNRAS, 430, 1392

Salyk, C., Blake, G. A., Boogert, A. C. A., & Brown, J. M. 2009, ApJ, 699, 330

Salyk, C., Blake, G. A., Boogert, A. C. A., & Brown, J. M. 2011, ApJ, 743, 112

Sicilia-Aguilar, A., Henning, T., & Hartmann, L. W. 2010, ApJ, 710, 597

Sicilia-Aguilar, A., Kim, J. S., Sobolev, A., et al. 2013, A&A, 559, A3

Sicilia-Aguilar, A., Roccatagliata, V., Getman, K., et al. 2015, A&A, 573, A19

Sitko, M. L., Day, A. N., Kimes, R. L., et al. 2012, ApJ, 745, 29

Smith, R. L., Pontoppidan, K. M., Young, E. D., Morris, M. R., & van Dishoeck, E. F. 2009, ApJ, 701, 163

Stelzer, B., Micela, G., Hamaguchi, K., & Schmitt, J. H. M. M. 2006, A&A, 457, 223

Stelzer, B., Robrade, J., Schmitt, J. H. M. M., & Bouvier, J. 2009, A&A, 493, 1109

Thi, W. F., Kamp, I., Woitke, P., et al. 2013, A&A, 551, A49

van Boekel, R., Min, M., Waters, L. B. F. M., et al. 2005, A&A, 437, 189

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199

van der Marel, N., Pinilla, P., Tobin, J., et al. 2015a, ApJ, 810, L7

van der Marel, N., van Dishoeck, E. F., Bruderer, S., Pérez, L., & Isella, A. 2015b, A&A, 579, A106

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

van der Plas, G., van den Ancker, M. E., Waters, L. B. F. M., & Dominik, C. 2015, A&A, 574, A75

van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771

Varnière, P., Blackman, E. G., Frank, A., & Quillen, A. C. 2006, ApJ, 640, 1110

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

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47
## Appendix A: Measured line centers, FWHM, integrated fluxes, and upper limits

<table><caption>Table A.1. CO ro-vibrational line fluxes, upper limits, and average line ratios.</caption><thead><tr><th>line</th><th>λ<sub>0</sub><br>[nm]</th><th>λ<sub>obs</sub><sup>a</sup><br>[nm]</th><th>ΔV<br>[km s<sup>-1</sup>]</th><th>Integrated flux<sup>b</sup><br>[10<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup>]</th><th>FWHM<sup>b</sup><br>[km s<sup>-1</sup>]</th><th></th></tr></thead><tbody><tr><td><sup>12</sup>CO 1 → 0 P(6)</td><td>4717.69</td><td>4717.73 ± 0.01</td><td>2.5 ± 0.6</td><td>13.0 ± 0.5</td><td>13.7 ± 0.4</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(7)</td><td>4726.73</td><td>4726.73 ± 0.01</td><td>0.4 ± 0.6</td><td>15.2 ± 0.4</td><td>14.7 ± 0.3</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(9)</td><td>4745.13</td><td>4745.15 ± 0.01</td><td>1.3 ± 0.6</td><td>14.3 ± 0.4</td><td>13.5 ± 0.4</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(10)</td><td>4754.50</td><td>4754.52 ± 0.01</td><td>1.3 ± 0.6</td><td>17.5 ± 0.4</td><td>15.1 ± 0.3</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(11)</td><td>4763.98</td><td>4764.01 ± 0.01</td><td>1.9 ± 0.6</td><td>14.8 ± 0.4</td><td>13.3 ± 0.4</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(12)</td><td>4773.58</td><td>4773.64 ± 0.01</td><td>3.8 ± 0.6</td><td>13.5 ± 0.3</td><td>13.4 ± 0.2</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(13)</td><td>4783.29</td><td>4783.33 ± 0.01</td><td>2.5 ± 0.6</td><td>14.7 ± 0.4</td><td>14.0 ± 0.3</td><td></td></tr><tr><td><sup>12</sup>CO 1 → 0 P(15)</td><td>4803.07</td><td>4803.10 ± 0.01</td><td>1.9 ± 0.6</td><td>13.3 ± 0.5</td><td>13.3 ± 0.4</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td><strong>14.5 ± 0.5</strong></td><td><strong>13.9 ± 0.3</strong></td><td><strong>average</strong></td></tr><tr><td><sup>13</sup>CO 1 → 0 R(6)</td><td>4715.04</td><td>4715.10 ± 0.01</td><td>3.8 ± 0.6</td><td>3.5 ± 0.2</td><td>11.8 ± 0.5</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 R(5)</td><td>4722.70</td><td>4722.72 ± 0.01</td><td>1.3 ± 0.6</td><td>4.3 ± 0.3</td><td>13.1 ± 0.6</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 R(4)</td><td>4730.47</td><td>4730.49 ± 0.01</td><td>1.3 ± 0.6</td><td>4.3 ± 0.2</td><td>13.4 ± 0.5</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 R(2)</td><td>4746.31</td><td>4746.36 ± 0.03</td><td>3.2 ± 1.9</td><td>3.3 ± 4.4</td><td>12.6 ± 2.8</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 R(0)</td><td>4762.56</td><td>4762.60 ± 0.01</td><td>2.5 ± 0.6</td><td>2.4 ± 0.3</td><td>11.2 ± 0.8</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 P(1)</td><td>4779.22</td><td>4779.27 ± 0.01</td><td>3.1 ± 0.6</td><td>2.6 ± 0.2</td><td>12.3 ± 0.9</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 P(2)</td><td>4787.71</td><td>4787.71 ± 0.02</td><td>0.0 ± 1.3</td><td>4.3 ± 0.5</td><td>15.9 ± 1.6</td><td></td></tr><tr><td><sup>13</sup>CO 1 → 0 P(4) +<br><sup>12</sup>CO 2 → 1 P(9)</td><td>4804.99</td><td>—</td><td>—</td><td>8.4 ± 0.8</td><td>—</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td><strong>3.5 ± 0.6</strong></td><td><strong>12.9 ± 0.6</strong></td><td><strong>average</strong></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(7)</td><td>4716.46</td><td>4716.50 ± 0.02</td><td>2.5 ± 1.3</td><td>1.2 ± 0.2</td><td>11.1 ± 1.6</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(6)</td><td>4724.03</td><td>4724.04 ± 0.02</td><td>0.6 ± 1.3</td><td>1.6 ± 0.2</td><td>10.8 ± 1.2</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(5)</td><td>4731.70</td><td>4731.75 ± 0.02</td><td>3.2 ± 1.3</td><td>1.4 ± 0.3</td><td>10.4 ± 1.2</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(3)</td><td>4747.34</td><td>4747.38 ± 0.03</td><td>2.5 ± 1.9</td><td>1.2 ± 0.3</td><td>11.4 ± 2.0</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(2)</td><td>4755.31</td><td>4755.36 ± 0.02</td><td>3.2 ± 1.3</td><td>0.9 ± 0.2</td><td>7.9 ± 1.2</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(1)</td><td>4763.38</td><td>c</td><td>—</td><td>1.5 ± 0.2</td><td>c</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 R(0)</td><td>4771.56</td><td>4771.62 ± 0.03</td><td>3.8 ± 1.9</td><td>0.9 ± 0.3</td><td>10.3 ± 2.2</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 P(1)</td><td>4788.22</td><td>4788.22 ± 0.02</td><td>0.0 ± 1.3</td><td>0.8 ± 0.2</td><td>5.5 ± 1.1</td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 P(3)</td><td>4805.29</td><td>c</td><td>—</td><td>1.1 ± 0.4</td><td>c</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td><strong>1.2 ± 0.1</strong></td><td><strong>9.6 ± 0.8</strong></td><td><strong>average</strong></td></tr><tr><td>C<sup>17</sup>O 1 → 0 R(0)</td><td>4716.96</td><td>—</td><td>—</td><td>&lt;0.5</td><td>—</td><td></td></tr><tr><td>C<sup>17</sup>O 1 → 0 P(1)</td><td>4733.62</td><td>—</td><td>—</td><td>&lt;0.5</td><td>—</td><td></td></tr><tr><td>C<sup>17</sup>O 1 → 0 P(3)</td><td>4750.71</td><td>4750.78 ± 0.04</td><td>4.4 ± 2.5</td><td>0.6 ± 0.2</td><td>9.9 ± 2.3</td><td></td></tr><tr><td>C<sup>17</sup>O 1 → 0 P(4)</td><td>4759.40</td><td>—</td><td>—</td><td>&lt;0.6</td><td>—</td><td></td></tr><tr><td>C<sup>17</sup>O 1 → 0 P(9)</td><td>4804.53</td><td>4804.56 ± 0.02</td><td>1.9 ± 1.2</td><td>0.9 ± 0.3</td><td>7.4 ± 2.3</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td><strong>0.8 ± 0.2</strong></td><td><strong>9 ± 1.6</strong></td><td><strong>average</strong></td></tr><tr><td><sup>12</sup>CO 2 → 1 R(0)</td><td>4715.72</td><td>4715.75 ± 0.01</td><td>1.9 ± 0.6</td><td>&gt;0.7</td><td>d</td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 P(1)</td><td>4732.65</td><td>4732.70 ± 0.01</td><td>3.2 ± 0.6</td><td>1.8 ± 0.2</td><td>12.2 ± 1.0</td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 P(3)</td><td>4750.01</td><td>4750.05 ± 0.01</td><td>2.5 ± 0.6</td><td>2.7 ± 0.2</td><td>11.4 ± 0.7</td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 P(4)</td><td>4758.86</td><td>4758.91 ± 0.01</td><td>3.1 ± 0.6</td><td>4.0 ± 0.2</td><td>12.7 ± 0.5</td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 P(6)</td><td>4776.89</td><td>4776.89 ± 0.03</td><td>0.0 ± 1.9</td><td>&gt;3.0</td><td>e</td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 P(7)</td><td>4786.07</td><td>4786.12 ± 0.01</td><td>3.1 ± 0.6</td><td>5.2 ± 0.3</td><td>14.8 ± 0.8</td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 P(9)</td><td>4804.78</td><td>blend <sup>13</sup>CO 1 → 0 P(4)</td><td>—</td><td>—</td><td>—</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td><strong>3.4 ± 0.7</strong></td><td><strong>12.8 ± 0.7</strong></td><td><strong>average</strong></td></tr><tr><td><sup>12</sup>CO 3 → 2 R(7)</td><td>4718.52</td><td>—</td><td>—</td><td>&lt;0.5</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 R(6)</td><td>4726.31</td><td>—</td><td>—</td><td>&lt;0.6</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 R(5)</td><td>4734.20</td><td>—</td><td>—</td><td>&lt;0.4</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 R(3)</td><td>4750.31</td><td>—</td><td>—</td><td>&lt;0.5</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 R(2)</td><td>4758.52</td><td>—</td><td>—</td><td>&lt;0.8</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 R(0)</td><td>4775.28</td><td>—</td><td>—</td><td>&lt;0.8</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 P(1)</td><td>4792.48</td><td>—</td><td>—</td><td>&lt;0.7</td><td>—</td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 P(2)</td><td>4801.24</td><td>—</td><td>—</td><td>&lt;1.0</td><td>—</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td><strong>&lt;0.6</strong></td><td><strong>—</strong></td><td><strong>average</strong></td></tr><tr><td colspan="7" style="text-align: center;">average ΔV all lines 2.1±1.2</td></tr><tr><td colspan="7" style="text-align: center;">average line ratios</td></tr><tr><td><sup>13</sup>CO 1 → 0 / <sup>12</sup>CO 1 → 0</td><td></td><td></td><td></td><td>0.24 ± 0.05</td><td></td><td></td></tr><tr><td>C<sup>18</sup>O 1 → 0 / <sup>12</sup>CO 1 → 0</td><td></td><td></td><td></td><td>0.08 ± 0.01</td><td></td><td></td></tr><tr><td>C<sup>17</sup>O 1 → 0 / <sup>12</sup>CO 1 → 0</td><td></td><td></td><td></td><td>0.05 ± 0.01</td><td></td><td></td></tr><tr><td><sup>12</sup>CO 2 → 1 / <sup>12</sup>CO 1 → 0</td><td></td><td></td><td></td><td>0.23 ± 0.05</td><td></td><td></td></tr><tr><td><sup>12</sup>CO 3 → 2 / <sup>12</sup>CO 1 → 0</td><td></td><td></td><td></td><td>&lt;0.04</td><td></td><td></td></tr></tbody></table>

**Notes.** Average line ratios were calculated from the average line fluxes. (a) Line centers are measured in the barycentric and radial velocity corrected spectra. Their value is the center and $3\sigma$ error of the Gaussian fit to the line. (b) Integrated flux and FWHM of the Gaussian fit and $1\sigma$ error; upper limits of the integrated line flux are $3\sigma$ assuming a line width of 10 km s$^{-1}$; the error in the average is the maximum between the standard error of the mean and $\frac{\sqrt{\sum \sigma_{\sigma}^2}}{N}$. (c) Detection but Gaussian fit not possible; (d) No meaningful FWHM due to the low sigma detection; (e) Detected but severely affected by the atmospheric transmission.
## Appendix B: Flat disk model

The computation of a flat disk model starts by calculating the expected integrated line flux of an annulus at radius R. This is done by multiplying the intensity by the solid angle of the annulus projected by the inclination i:

$$F(R) = I(R) \Omega(R)_{\text{projected}}. \qquad (B.1)$$

Here R is the mid-point of two grid points in the radial grid

$$R = 0.5 (R_j + R_{j+1}), \qquad (B.2)$$

and

$$\Omega(R)_{\text{projected}} = \pi \cos(i) (R_{j+1}^2 - R_j^2)/D^2, \qquad (B.3)$$

where D is the distance to the source. The integrated line flux for each cell in the azimuthal direction is then determined by dividing the total integrated flux of the annulus by the number of points of the azimuthal grid ($N_{\theta}$):

$$F(R, \theta) = F(R)/N_{\theta}. \qquad (B.4)$$

We used 50 to 100 points in the radial direction and 1000 points in the azimuthal direction. The local line-profile $\phi$ of a grid point in R and $\theta$ is then obtained by convolving the integrated flux of the cell with a normalized Gaussian kernel, with a FWHM equal to the spectral resolution convolved with the turbulent and thermal broadening:

$$\phi(R, \theta, \nu) = F(R, \theta) * \phi_{\text{Gauss}}(\nu). \qquad (B.5)$$

The line-profile of each cell in the azimuthal direction $\theta$ is then velocity shifted to the expected local Keplerian velocity shift

$$\Delta V = \cos(\theta) \sin(i) \sqrt{\frac{GM_{\star}}{R}}, \qquad (B.6)$$

thus obtaining a $\phi(R, \theta, \nu)_{\text{shifted}}$ for each cell.

When no slit effects are taken into account, the 1D spectrum of the whole disk is obtained by summing the contributions of

each azimuthal cell in each annulus:

$$\phi(R, \nu) = \sum_{j=1}^{N_{\theta}} \phi(R, \theta, \nu)_{\text{shifted}}, \qquad (B.7)$$

and summing the spectra of all the annuli in the radial direction

$$\phi(\nu) = \sum_{j=1}^{N_R} \phi(R, \nu). \qquad (B.8)$$

To generate 3D channel maps the $\phi(R, \theta, \nu)_{\text{shifted}}$ of each cell is sampled in a Cartesian data cube with coordinates X, Y, $\nu$ where

$$X = R \cos(\theta), \qquad (B.9)$$

$$Y = R \sin(\theta) \cos(i). \qquad (B.10)$$

When the effect of the slit is taken into account, the 1D spectrum is extracted from the 3D (X, Y, $\nu$) channel map data cube generated by the model. First, the image in each velocity channel is convolved with a Gaussian beam of FWHM 206 mas to model the spatial resolution. Then, to simulate the effect of the slit, the 3D data-cube is rotated to account for the position angle of the disk on the sky and the slit orientation. Then a 2D spectrum is obtained from the 3D data by summing the pixels inside a 0.2" vertical aperture. This 2D disk model spectrum is scaled such that in the extracted 1D disk spectrum, the peak of the line is equal to the peak of the flux in the normalized observed composite 1D spectrum.

To calculate the spectroastrometric signature, a synthetic star+disk 2D spectrum is created by adding to the 2D disk spectrum a 2D star spectrum broadened by the PSF-FWHM. The 2D star spectrum is constructed such that the continuum in the extracted 1D spectrum is equal to 1. The model 2D star + disk spectrum is finally re-binned in the spatial and spectral directions such that its spatial and spectral pixel scales are the same as the CRIRES data. With this synthetic 2D star + disk spectrum, the theoretical spectroastrometric signature was measured using the formalism of Pontoppidan et al. (2011).

<figure><img src="image_27.png" /><figcaption>Fig. B.1. Line-profiles predicted for models around the best solution of the power-law intensity model: a) effect of varying the inner radius with the outer radius fixed; b) effect of varying the inner radius keeping $\alpha$ fixed. In the models, $\alpha$ or $R_{\text{out}}$ are adjusted such that $I(R_{\text{out}}) = 0.01 \times I(R_{\text{in}})$; c) line-profile for an intensity distribution with a sharp increase at 1.2 AU (in blue) and the line-profile of an intensity distribution that grows as a power-law from 0.01 AU to 1.2 AU (in red). In both models the intensity decreases with $\alpha = -1.8$ at $R > 1.2$ AU. Error bars on the composite spectrum are $1\sigma$.</figcaption></figure>
## Appendix C: Calculation of the Bayesian probability

In each model, a $\chi^2$ was calculated for each observational dataset (i.e., one $\chi^2$ for the $^{12}$CO P(9) line-profile, one $\chi^2$ for the $^{12}$CO rotational diagram, etc.) using

$$ \chi^2 = \frac{1}{N-1} \sum_i (Y_{\text{model } i} - Y_{\text{obs } i})^2 / \sigma_{Y_{\text{obs } i}}^2, \qquad (C.1) $$

where, Y corresponds to the line flux per velocity channel, in the case of the spectrum, and the Y axis of the plot, in the case of the rotational diagram. N is the number of channels in the spectrum or the number of data points in the rotational diagram. For the line-profiles we used the velocity channels from -15 to 15 km s$^{-1}$. This enabled us to cover the wings of the line and a small part of the continuum.

As there are six observational datasets, six $\chi^2$ values were calculated per model. Because the numerical value of $\chi^2$ can be very different for the rotational diagram and the line-profile, before calculating the combined $\chi^2$, the $\chi^2$ of each observational dataset (i.e. line-profile, rotational diagram) was normalized by the minimum values of $\chi^2$ of that observational dataset in the entire grid. The sum of the six normalized $\chi^2$ gave the final $\chi^2$ for a model. The Bayesian probability

$$ p = \exp(-\chi^2/2) \qquad (C.2) $$

was calculated for each model, and finally p was divided by the sum of all p. In this way, a normalized Bayesian probability p was obtained for each model. The 1D probability for each free parameter was calculated by summing the normalized p of all the models containing a particular value of the parameter in question. Similarly, 2D probability distributions were constructed by summing the normalized p of all the models containing the pair of values for the free parameters in the plot.

<figure><img src="image_28.png" /><figcaption>Fig. C.1. Bayesian probability plots for a sub-sample of the grid in which only the models with $\alpha_{N_{H,inner}} \le +1$ are considered (54 000 models).</figcaption></figure>
$$ \text{BEST-FIT GRID MODEL IF } \alpha_{N_{\text{H inner}}} \le +1.0: \quad \alpha_{N_{\text{H inner}}} = 0.0 \quad \delta_{\text{gas}} = 10^{-2} $$

<figure><img src="image_29.png" /><figcaption>**Fig. C.2.** Surface density, temperature, CO optical depth, flux density, cumulative line flux, rotational diagrams, line-profiles of the ¹²CO P(9), ¹³CO R(4), and C¹⁸O R(6) emission for best-fit grid model when $\alpha_{N_{\text{H inner}}} \le +1.0$. The model is plotted in red and the observations in black. Observed line-profiles are displayed in flux units after continuum subtraction with 3$\sigma$ error bars. The two branches seen in the rotational diagram correspond to the R and P branches of CO ro-vibrational emission. The rightmost panels compare the normalized theoretical line-profiles with the observed composite line-profile of each CO isotopolog with a 1$\sigma$ error bar.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. C.3.** Bayesian probability plots for the grid using only the ¹²CO and ¹³CO data (i.e. no C¹⁸O data). The models suggest a surface density drop of at least a factor 100 in the inner 6 AU, and an increasing surface density profile with radius (i.e. a power-law with a positive exponent).</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. C.4.** Bayesian probability plots for a sub-grid of models (28 800), in which we varied the radius of the gas density drop ($R_{gap max}$) between 4.0 and 6.0 AU. $R_{gap max}$ down to 5 AU are compatible with the data. The most likely value for the gas density drop is 6 AU, a radius similar to the dust density drop.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. C.5.** Flux density, cumulative line flux, rotational diagrams, and the line-profiles of the ¹²CO P(9), ¹³CO R(4), and C¹⁸O R(6) emission, for the model with the best combined fit to the rotational diagrams and line-profiles, extrapolating the surface density and temperature profile down to 0.1 AU. The model is shown in red and the observations in black. The observed line-profiles are displayed in flux units after continuum subtraction with 3σ error bars. The two branches seen in the rotational diagram correspond to the R and P branches of CO ro-vibrational emission. The rightmost panels compare the normalized theoretical line-profiles with the observed composite line-profile of each CO isotopolog with a 1σ error bar.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Fig. C.6.** Surface density and temperature profile, and predicted $v = 1 \rightarrow 0$ $^{12}$CO P(9), $^{13}$CO R(4), C$^{18}$O R(6), and $v = 2 \rightarrow 1$ $^{12}$CO P(3) and $^{12}$CO P(4) line-profiles for a model with $N_H$ at $R = 6$ AU three times larger than the best-fit grid model ($N_H$ ($R=6$ AU) = $3 \times 10^{23}$ cm$^{-2}$). The higher $N_H$ enables to describe the $v = 2 \rightarrow 1$ $^{12}$CO P(3) and $^{12}$CO P(4) line-profiles while keeping a good fit to the $v = 1 \rightarrow 0$ lines. The model has the same temperature structure and same surface density at $R < 6$ as the best-fit grid model (thus $\delta_{gas} = 3.3 \times 10^{-3}$). The cumulative flux plot shows that the $v = 2 \rightarrow 1$ lines are dominated by the contribution at $6 < R < 10$ AU. Errors in the plot are $3\sigma$, and the dashed horizontal line is the $5\sigma$ limit.</figcaption></figure>
<table><caption>**Table C.1.** Transition disks with a quantified gas density drop inside the dust cavity.</caption><thead><tr><th>Star</th><th>R<sub>cav</sub> dust [au]</th><th>δ<sub>dust</sub></th><th>λ</th><th>Ref.</th><th>R<sub>cav</sub> gas [au]</th><th>δ<sub>gas</sub></th><th>Σ<sub>gas</sub> inner [g cm<sup>-2</sup>]</th><th>g/d<sub>outer</sub></th><th>Gas tracers</th><th>Ref.</th><th>M&#773; [M<sub>☉</sub> yr<sup>-1</sup>]</th><th>Ref.</th></tr></thead><tbody><tr><td>HD 139614</td><td>6</td><td>10<sup>-4</sup></td><td>near-IR, mid-IR sub-mm</td><td>1</td><td>5-6</td><td>10<sup>-2</sup>-10<sup>-4</sup></td><td>7 × 10<sup>-5</sup>-2.4 × 10<sup>-3</sup></td><td>1-100</td><td>1<sup>2</sup>CO, <sup>13</sup>CO, <sup>C</sup><sup>18</sup>O ro-vib.</td><td>this work</td><td>1 × 10<sup>-8</sup></td><td>8</td></tr><tr><td>RXJ1615-3255</td><td>20</td><td>10<sup>-6</sup></td><td>sub-mm</td><td>2</td><td>20</td><td>&gt;10<sup>-4</sup></td><td>5 × 10<sup>0</sup>-1 × 10<sup>-1</sup></td><td>100</td><td><sup>12</sup>CO sub-mm</td><td>2</td><td>3.2 × 10<sup>-9</sup></td><td>9</td></tr><tr><td>SR21</td><td>25</td><td>10<sup>-6</sup></td><td>sub-mm</td><td>3</td><td>25</td><td>5 × 10<sup>-2</sup><br>≤10<sup>-5</sup></td><td>1 × 10<sup>1</sup>-1 × 10<sup>0</sup><br>1 × 10<sup>-1</sup>-1 × 10<sup>-2</sup></td><td>100</td><td><sup>12</sup>CO, <sup>13</sup>CO sub-mm</td><td>3</td><td>1.3 × 10<sup>-8</sup></td><td>9</td></tr><tr><td>DoAr 44</td><td>32</td><td>10<sup>-2</sup></td><td>sub-mm</td><td>3</td><td>32</td><td>10<sup>-2</sup><br>≤10<sup>-4</sup></td><td>5 × 10<sup>-1</sup>-1 × 10<sup>-1</sup><br>1 × 10<sup>-2</sup>-5 × 10<sup>-4</sup></td><td>100</td><td><sup>12</sup>CO, <sup>13</sup>CO sub-mm</td><td>3</td><td>6.3 × 10<sup>-9</sup></td><td>9</td></tr><tr><td>HD 135344 B</td><td>30, 45</td><td>10<sup>-2</sup></td><td>near-IR, sub-mm</td><td>4, 5, 6</td><td>45</td><td>10<sup>-1</sup></td><td>1 × 10<sup>-2</sup>-2 × 10<sup>-2</sup></td><td>4</td><td><sup>12</sup>CO ro-vib. &amp; sub-mm [OI] 63 µm</td><td>6<sup>a</sup></td><td>2 × 10<sup>-8</sup></td><td>10</td></tr><tr><td>HD 135344 B</td><td>40</td><td>2 × 10<sup>-4</sup></td><td>sub-mm</td><td>5, 3</td><td>30</td><td>2 × 10<sup>-4</sup></td><td>7 × 10<sup>-1</sup>-2 × 10<sup>-2</sup></td><td>80</td><td><sup>12</sup>CO, <sup>13</sup>CO sub-mm</td><td>3<sup>a</sup></td><td>2 × 10<sup>-8</sup></td><td>10</td></tr><tr><td>LkCa 15</td><td>45</td><td>10<sup>-5</sup></td><td>sub-mm</td><td>2</td><td>45</td><td>10<sup>-1</sup></td><td>1 × 10<sup>3</sup>-1 × 10<sup>1</sup></td><td>100</td><td><sup>12</sup>CO sub-mm</td><td>2</td><td>4 × 10<sup>-9</sup></td><td>9</td></tr><tr><td>IRS 48</td><td>60</td><td>10<sup>-3</sup></td><td>sub-mm</td><td>3</td><td>25</td><td>≤10<sup>-3</sup></td><td>1 × 10<sup>2</sup>-5 × 10<sup>-4</sup></td><td>12</td><td><sup>12</sup>CO, <sup>13</sup>CO sub-mm</td><td>3</td><td>3.2 × 10<sup>-9</sup></td><td>11</td></tr><tr><td>J1604-2130</td><td>70</td><td>10<sup>-1</sup></td><td>sub-mm</td><td>2</td><td>30</td><td>10<sup>-5</sup></td><td>1 × 10<sup>2</sup>-1 × 10<sup>-4</sup></td><td>100</td><td><sup>12</sup>CO sub-mm</td><td>2</td><td>&lt;10<sup>-11</sup></td><td>12</td></tr><tr><td>HD 142527<sup>b</sup></td><td>130</td><td>10<sup>-5</sup></td><td>sub-mm</td><td>7</td><td>90</td><td>2 × 10<sup>-2</sup></td><td>1 × 10<sup>2</sup>-4 × 10<sup>0</sup></td><td>100</td><td><sup>12</sup>CO, <sup>13</sup>CO, <sup>C</sup><sup>18</sup>O sub-mm</td><td>7</td><td>2 ± 1 × 10<sup>-7</sup></td><td>13</td></tr></tbody></table>

**Notes.** R<sub>cav</sub> dust is the dust cavity radius resolved by near-IR or sub-mm observations. When two values appear, the first corresponds to the small dust particles in the near-IR and the second to the large dust grains in the sub-mm. R<sub>cav</sub> gas is the reference radius for the drop in the gas surface density. An additional R<sub>cav</sub> gas is given if a second drop in the gas surface density is suggested in the literature. Σ<sub>gas</sub> inner is generally described as a power-law, the interval of values given correspond to Σ<sub>gas</sub> inner at the inner disk's radius R<sub>in</sub>, and at R<sub>cav</sub> gas. If two drops for the surface density are suggested in the literature, two intervals for the surface density in the inner disk are given. <sup>(a)</sup> Two models are available in the literature for HD 135344B. Although they have dust mass estimations similar within a factor 2 (1-2 × 10<sup>-4</sup> M<sub>☉</sub>) the gas mass and gas surface density in the outer disk differ. Carmona et al. (2014) find a gas-to-dust mass ratio of 4 and a gas mass of 8 × 10<sup>-4</sup> M<sub>☉</sub> for the outer disk, van der Marel et al. (2016) find a gas-to-dust mass ratio of 80 and a gas mass of 1.5 × 10<sup>-2</sup> M<sub>☉</sub>. Both models also used a different disk outer radius, 200 AU in Carmona et al. (2014) and 125 AU in van der Marel et al. (2016); <sup>(b)</sup> δ<sub>dust</sub> was calculated from the surface density profile provided by Perez et al. (2015), using a gas-to-dust mass of 100 at R > 130 AU, such that the dust mass between 0.2 and 6 AU is 10<sup>-9</sup> M<sub>☉</sub> as reported by the same authors.

**References.** (1) Matter et al. (2016); (2) van der Marel et al. (2015b); (3) van der Marel et al. (2016); (4) Andrews et al. (2011); (5) Muto et al. (2012); (6) Carmona et al. (2014); (7) Perez et al. (2015); (8) Garcia Lopez et al. (2006); (9) Manara et al. (2014); (10) Sitko et al. (2012); (11) Follette et al. (2015); (12) Mathews et al. (2012); (13) Mendigutía et al. (2014).