© ESO, 2017

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¹. 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 horseshoe 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 (δ_gas) inside the dust cavity, with δ_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_⊙/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 ². and a dust density drop δ_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_⊙ 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 towards HD 139614 obtained with the ESO/VLT CRIRES instrument (Kaüfl 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^-6μm in the wavelength direction (0.14 km s^-1 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^-1.

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 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 (1)that gave the smallest χ² 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 χ² 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 τ was estimated using Here is median of the standard star flux at the wavelengths within 95% and 100% of the atmospheric transmission and σ 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)⁴. 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^-1, 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_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^-1 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^-1, and the 1D spectrum was continuum subtracted and normalized by the peak flux (median of the flux within ±2 km s^-1).

3. Observational results

We have detected υ = 1 → 0 ¹²CO, ¹³CO, C¹⁸O, C¹⁷O, and υ = 2 → 1 ¹²CO emission lines. The υ = 3 → 2 ¹²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 → 0 ¹²CO emission. To keep the notation short in the remaining of the paper, we mean by ¹²CO, ¹³CO, C¹⁸O emission υ = 1 → 0¹²CO, ¹³CO, C¹⁸O emission unless otherwise specified.

We reached a 3σ sensitivity of 2 × 10^-16 erg s^-1 cm^-2 for a line width of 3.3 km s^-1, and 1.2 × 10^-15 erg s^-1 cm^-2 for a line width of 20 km s^-1 (equivalent widths of 0.01 and 0.075 Å respectively). We achieved a spectral resolution of ~3.3 km s^-1 (R ~ 10⁵) as measured in an unresolved unsaturated sky-absorption line. The centers of the CO emission lines in the barycentric and radial-velocity-corrected spectra are located on average at v = 2 ± 1 km s^-1 (see Table A.1). As this value is close to zero and is lower than the uncertainty of ±2.3 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 ¹²CO composite line-profile has a flat top and does not display evidence of asymmetries⁵. The composite ¹³CO and C¹⁸O lines are single peaked. Some asymmetric sub-structures are present in both lines but they are consistent with noise.

C¹⁸O emission is, at the 2σ level, 4 km s^-1 narrower than ¹²CO emission. ¹³CO and 2 → 1 ¹²CO lines are 1 km s^-1 narrower than the ¹²CO line, at the 1σ level. To further test whether the ¹²CO, ¹³CO and C¹⁸O line-profiles are different, we ran a two-sample Kolmogorov-Smirnov (K-S) test⁶ on the composite 1D spectra in the ±15 km s^-1 interval, after normalization by the line peaks. The K-S significance between the ¹²CO and C¹⁸O line-profiles is 8%, between the ¹²CO and ¹³CO line-profile is 30%, and between the C¹⁸O and ¹³CO line-profile is 97%. The K-S test indicates that the ¹²CO and C¹⁸O profiles are different (the C¹⁸O is narrower), and that statistically the ¹³CO profile resembles the C¹⁸O profile more than the ¹²CO profile.

The 1σ 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σ error quoted is an average value. No displacement of the photocenter centroid is detected at the position of the ¹²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 ± 10 mas) and the telluric standard (172 ± 10) are consistent within the errors, which means that there is no evidence of extended continuum emission at 4.7 μm. We measured a stacked continuum PSF-FWHM of 2.40 ± 0.05 pixels (1σ), equivalent to 206 ± 4 mas or 27 ± 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 location of the line appears constant as a function of the wavelength. This directly indicates that there is no ¹²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 ¹²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 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: (4)extending from the an inner radius R_in to an outer radius R_out, where I₀ is the intensity at R_in which is assumed initially to be 1. The exponent α is obtained for each pair of R_in and R_out such that I(R_out) = 0.01 × I₀. In this model, all the physics of the excitation of the line is in the exponent α. 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 ¹²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_⊙, an inclination i = 20°, a PA = 292° (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 ± 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 reduced statistic (assuming three free parameters: R_in, R_out, and I₀) of the model of the composite ¹²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 ¹²CO composite line-profile. The model that displays the smallest has R_in = 1.2 AU and R_out = 15 AU and an α exponent of the intensity −1.8.

Although disk models with R_out as large as 20 AU provide a good fit to the ¹²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σ 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σ). 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σ 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σ 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 ¹²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 <v< 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 ≥ 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 ¹²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 (α adjusted such that I(R_out) = 0.01 × I(R_in)). In panel (b) the line-profiles with α fixed to −1.8 (R_out set such that I(R_out) = 0.01 × I(R_in)). Depending on the value of α, models with R_in as low as 0.3 AU can be compatible with the observed ¹²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₀ 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_Imax = 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 ¹²CO line-profile shape.

4.1.2. A continuous or a gapped gas distribution?

The ¹²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 ¹²CO line-profile could be described by an intensity distribution with a gap.

The ¹²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 ¹²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¹⁸O and C¹⁷O indicates that the emitting medium must be 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¹⁸ to 10²⁵ cm^-2 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^-1. We used a ¹²CO abundance of 10^-4 (N_¹²CO = 1.0 × 10^-4N_H), a ¹²CO/¹³CO ratio of 100 and a ¹²CO/C¹⁸O ratio of 690 following Smith et al. (2009)⁷.

The output of a slab model is in units of line flux per steradian 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., α = 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_in = 0, R_out = 6.5 AU. We therefore used an average emitting region radius of 6.5 AU and a distance of 140 pc⁸ to determine the solid angle for the 1D slab model⁹.

For each N_H and T_gas slab model, a rotational diagram for each CO isotopolog was calculated. We used the X and Y coordinates Y = ln(F_ul/νg_uA_ul) and X = E_u. Here F_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_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 . Here σ_{Yobs i} is the difference between Y calculated using the observed line flux and Y calculated using the observed line flux plus 3σ. The N−4 corresponds to three degrees of freedom (N_H, T_gas, Ω), and N the number of data points (8 for ¹²CO, 7 for ¹³CO and 9 for C¹⁸O). In Fig. 6 we display the contour plots for each CO isotopolog and one 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_gas ~ 450 K and N_H> 10²⁰ cm^-2. ¹³CO is best reproduced by lower temperatures (T_gas ~ 380 K) and N_H column densities of at least 5 × 10²² cm^-2. The C¹⁸O is best fit by even lower temperatures ~350 K and N_H column densities higher than 5 × 10²² cm^-2. The errors on T_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 × 10²¹ cm^-2 and T_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 contours, because solutions are degenerate with respect to N_H, giving only lower limits for the column density.

4.2.1. Thermally excited or UV-pumped emission?

The detection of υ = 2 → 1 ¹²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 υ = 2 → 1 ¹²CO/υ = 1 → 0 ¹²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 υ = 3 → 2 ¹²CO emission and the upper limit of 0.04 on the υ = 3 → 2 ¹²CO /υ = 1 → 0 ¹²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 υ = 3 → 2 ¹²CO/υ = 1 → 0 ¹²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 υ = 1 → 0 ¹²CO, ¹³CO, and C¹⁸O rotational diagrams and calculated the LTE υ = 2 → 1 ¹²CO and υ = 3 → 2 ¹²CO emission. We found that the observed υ = 2 → 1¹²CO emission and the non-detection of the υ = 3 → 2 ¹²CO emission can be well described by an LTE model, with T_gas = 460 K and N_H = 10²² cm^-2, in which the rotational temperature is equal to the vibrational temperature¹⁰. We show this model in Fig. 8, where we display the rotational diagrams of υ = 1 → 0 ¹²CO and υ = 2 → 1 ¹²CO emission of the model and the observations. This model generates υ = 3 → 2 ¹²CO emission lines with integrated fluxes on the order of 10^-18−10^-17 erg s^-1 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 isopotologs, suggest that the emission of each isotopolog is produced at different radial distances and/or vertical heights. The narrower and colder ¹³CO and C¹⁸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_H ∝ R^α_NH and T_gas ∝ R^α_Tgas) and both distributions have no discontinuities (i.e. no gaps nor density drops or jumps in the temperature), and if the ¹²CO/¹³CO and ¹²CO/C¹⁸O abundance ratios are constant, then, a disk model able to reproduce the ¹²CO line-profile and ¹²CO rotational diagram would generate ¹³CO and C¹⁸O lines that are too broad, too strong, and too warm to be consistent with the observations (see Fig. 9). In the next sub-section we provide the details of the model. Since ¹³CO and C¹⁸O emission is optically thin up to relatively high column densities (N_H ~ 10²² cm^-2 and N_H ~ 10²³ cm^-2 respectively), the 1D single power-law models generate ¹³CO and C¹⁸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 C¹⁸O and ¹³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 C¹⁸O and ¹³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: We allowed the models to have a reduced surface density at R< 6 AU by a factor δ_gas(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 ¹²CO line with a power-law intensity suggested that 1 AU is the radius of the maximum intensity.

We chose α_{NH inner} to range from −2.5 up to +3 to cover a wide range of possible surface density distributions in the inner disk¹¹. As Matter et al. (2016) found a dust depletion factor of 10^-4, we tested models with δ_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_gas and N_H using the CO slab model previously described, (8)The local broadening of the line is the convolution of the turbulent broadening (0.1 km s^-1), the local thermal broadening and the spectral resolution¹². 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_H(R) and T_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_H(R) and T_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_{H (R = 6 AU)},α_{NH inner},α_{NH outer},δ_gas,T_{(R = 1 AU)} and α_Tgas. Each model produces integrated line-fluxes and rotational diagrams for the three CO isotopologs and synthetic line-profiles at the CRIRES resolution for the ¹²CO P(9) line at 4745.13 nm, ¹³CO R(4) at 4730.47 nm, and the C¹⁸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 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 δ_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 ~ 10²³ cm^-2. This gas column density is similar to the dust column density at R = 6 AU found by Matter et al. (2016)¹³. 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 × 10¹⁹ cm^-2 and N_H = 10²¹ cm^-2 (7 × 10^-5−2.4 × 10^-3 g cm^-2). 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 δ_gas could be as large as 10^-4 depending on the total gas mass of the outer disk.

The Bayesian probability distributions indicate that the preferred values for the surface density exponent in the inner 6 AU are flat or positive, with α_{NH 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 α_{NH 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^-2, α_{NH inner}(R< 6 AU) = + 2.0, α_{NH outer}, δ_gas(R = 6 AU) = 10^-2, T₀(R = 1 AU) = 675 K, and α_Tgas = −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^-1, 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 line wings that are too strong. The effect is also seen in the ¹³CO and C¹⁸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 δ_gas = 10^-2 at R = 6 AU. As soon as α_{NH inner} = 0 is reached, the strong high-velocity line wings of the ¹³CO and lines C¹⁸O disappear. The ¹²CO line-profile hardly changes in all the models with α_{NH inner} –1.0 to +1.0 because it is optically thick. However, when α_{NH inner} is higher than +1.5, ¹²CO becomes optically thin at R< 2 AU, the wings of the ¹²CO line become weaker, and the ¹²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 ¹²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 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²¹ to 10¹⁷ cm^-2. The normalized theoretical ¹²CO P(9) spectrum was compared with the high S/N composite ¹²CO spectrum.

Figure 15 display 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¹⁸ cm^-2 (yellow and red, δ_gas< 10^-5) would have been seen in the ¹²CO composite spectrum as a line-profile with shoulders at ±15 km s^-1.

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

¹²CO ro-vibrational emission requires relatively low column densities (N_CO ~ 10¹⁵ cm ^-2) 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 × 10¹⁹ cm^-2 (1.2 × 10^-4 g cm^-2) would have produced line-profile wings (v> 15 km s^-1) 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 × 10¹⁵ cm^-2.

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¹⁵ cm^-2 (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 ≤ 10¹⁹ cm^-2 at R ≤ 1 AU, assuming standard abundances. This value is consistent 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 α_{NH inner} of +2 or +3. However, such high α_{NH inner} are hard to justify physically¹⁴. To test the robustness of the modeling conclusions, we recalculated the Bayesian probability plots for a sub-grid selecting only models with α_{NH inner} ≤ + 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 α_{NH inner} ≤ + 1 sub-grid. The model has δ_gas = 10^-2 and has a flat (α_{NH inner} = 0.0) surface density at R < 6 AU¹⁵.

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¹⁶. 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¹⁷. 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¹⁸. 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¹⁹. The column of C¹⁸O (thus N_H) could be underestimated because the 1D models use a higher average temperature. However, the column density in the inner 6 AU should be lower than N_H = 10²¹ cm^-2. Higher N_H, even with T_{0 (R = 1 AU)} = 575 K (100 K lower than the best-fit grid model), would generate ¹³CO and C¹⁸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 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 ¹²CO line-profile, and the fact that the average temperatures of ¹³CO and C¹⁸O lines (380 and 350 K respectively) are similar to that of the power-law temperature at 6 <R < 10 AU of the best-fit grid model (350–300 K, where most of the ¹³CO and C¹⁸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)²⁰. 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²¹.

The model grid was calculated assuming R_in = 1 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 AU (α_{NH inner}) between +1 and +3. Models with α_{NH inner} smaller than +1 gave too strong line wings for the ¹²CO emission. The best-fit grid model, which has α = + 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 AU derived in Sect. 4.5.

In the grid of models we did not fit the υ = 2 → 1 ¹²CO lines. We checked the predicted υ = 2 → 1 ¹²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 υ = 2 → 1 lines, but the model has weaker υ = 2 → 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 × 10²³ cm^-2), the same density profile at R < 6 AU (thus δ_gas = 3.3 × 10^-3), and the same temperature profile described the υ = 2 → 1 lines well while having a good fit to the υ = 1 → 0 lines. The model indicates that the υ = 2 → 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 AUto 15 AU in HD 139614, whichmeans that the dust gap observed in IR-interferometry datacontains molecular gas.

2. C¹⁸O lines are a few km s^-1 narrower than ¹²CO lines.

3. ¹³CO and C¹⁸O emission are 50–100 K colder than ¹²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 × 10¹⁹ and 10²¹ cm^-2 and has a distribution that most likely is flat or that increases with radius.

7. The 5σ upper limit on the CO column density N_CO at R ≤ 1 is 5 × 10¹⁵ cm^-2, which corresponds to a gas column density N_H < 5 × 10¹⁹ 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 α_{NH inner} ≤ 1.0.

We find that the α_inner = + 0.6 deduced for the dust is consistent with the range of α_{NH 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 (δ_gas and δ_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 × 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²². 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^-1 cm^-2) 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^-1) that are observed in primordial disks²³. 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 ~ 10¹⁵ cm^-2), 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σ upper limit of N_H ~ 5 × 10¹⁹ cm^-2 on the column of gas at R < 1 AU (assuming standard abundances) derived for HD 139614 is significantly lower than the N_H = 10²⁴−10²⁶ cm^-2 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^-1, while in HD 139614 (and other transition disks) it is only a few km s^-1. 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 υ = 1 → 0 CO and υ = 2 → 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 (logT = 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 υ = 2 → 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²⁴. The largest portion of transition disks shows stronger depletions in the dust than in the gas (δ_dust< δ_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_⊙/yr in HD 139614 (Garcia Lopez et al. 2006) is at the high end of the mass-loss rates (10^-8−10^-10 M_⊙/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_⊙ star (R_c,EUV ≃ 1.8 × M_∗/M_⊙ 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 probability given the EUV photoevaporation timescales (10⁵ 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 × 10^-14 erg s^-1 cm^-2, which translates into an X-ray luminosity of 1.2 × 10²⁹ 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 × 10^-10 M_⊙/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 δ_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²⁵. 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²⁶; (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 (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_⊙. 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²⁷.

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 δ_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²⁸.

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_⊙ 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 ¹²CO line-profile is symmetric and distortions are not detected at the 1σ 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 as the changes in opacity that are due to migrating solids reaching 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 features discovered in disks, it is not clear whether they are able to explain the gas density drop that the CO ro-vibrational emission reveals in HD 139614 (and other transition disks). The low gas column-density detected inside the dust cavities of transition 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 profile within 6 au of HD 139614 is indeed the result of a giant planet of mass between 1 and 2 M_J 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 trapping of the particles depends on how well they are coupled to the gas and therefore, it depends on their size and local gas surface density. Quantitatively, particles with size a_opt = 2Σ_gas/πρ_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 ρ_d is the internal density of the dust grains, typically 1–3 g cm^-3, and Σ_g is the gas surface density. Using N_H = 10¹⁹−10²¹ cm^-2 (2.4 × 10^-5−2.4 × 10^-3 g cm^-2)²⁹ 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 surface 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 (Ṁ_⋆) between 10^-9 and 10^-8 M_⊙ 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 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 Ṁ_⋆ can have very different gas surface densities in their inner regions. This is the case of IRS 48 and RXJ1615, both of which have Ṁ_⋆ ~ a few × 10^-9 M_⊙ 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, Ṁ_disk, is related to the surface density Σ and radial velocity v_R of the gas at radius R through Ṁ_disk = −2πRv_R(R) Σ(R). If we now assume that Ṁ_disk ~ Ṁ_⋆, then at a radius of 1 AU and for an accretion rate of 10^-8 M_⊙ yr^-1, we find that | v_R | should extend from about 70 km s^-1 down to 0.07 km s^-1 for Σ in the range [10^-3−1] 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 values expected in classical viscous disk models. These models tell us that |v_R| ~ αh²v_K, were h is the disk aspect ratio and α is the dimensionless turbulent viscosity of the disk. If the disk magnetic field is able to sustain vigorous turbulence in the cavities of transition disks, then |v_R| can reach at most ~10^-5v_K. This is at least two orders of magnitude lower than the above range of radial velocities obtained when assuming Ṁ_disk ~ Ṁ_⋆.

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 densities would produce. For HD 139614, the Ṁ_star = 10^-8 M_⊙ yr^-1 would suggest an N_H of 10²⁵−10²⁶ cm^-2 at R ≤ 1 AU for α ranging from 10^-2 to 10^-3. These gas columns are high enough for the CO to self-shield against UV photodissociation and produce strong CO ro-vibrational emission from R < 1 AU (not seen in our CRIRES spectrum). For example, if N_H ~ 10²⁵ cm^-2 at R < 1 AU in a disk without dust in the cavity (best case for UV-photodissociation) around a T_eff = 10 000 K central star (i.e., brighter in UV than HD 139614) N_CO would be 10²¹ cm^-2 at R < 1 AU (see Fig. 8 in Bruderer 2013), a value much higher than the 5σ limit of N_CO = 5 × 10¹⁵ cm^-2 from our CRIRES data. According to the Bruderer (2013) models (for a T_eff = 10 000 K), to have N_CO < 10¹⁵ cm^-2, N_H needs to below 10²² cm^-2 at R < 1 AU for an inner disk without dust. For HD 139614, which has a T_eff ~ 7800 K (and lower UV, accordingly) and some dust in the inner disk³⁰, N_H should be even lower to enable the efficient 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 Ṁ_disk ~ Ṁ_⋆. 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 turbulent 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 generate fast radial flows inside the cavity of its transition disk. Pontoppidan et al. (2011) detected in HD 142527 a ¹²CO ro-vibrational spectrum and spectro-astrometry signature displaying asymmetries, indicating non-Keplerian contributions to the emission. Additional resolved observations of the gas kinematics inside the cavities of transition disks will help assess the occurrence of this second scenario. In our case of HD 139614, the ¹²CO composite line-profile is symmetric, smooth, and consistent with Keplerian motion, which suggests the first scenario, namely that Ṁ_disk (1 AU) is most likely different from Ṁ_⋆.

6. Summary and conclusions

We have obtained VLT/CRIRES high-resolution spectra (R ~ 90 000) of CO ro-vibrational emission at 4.7 μm in HD 139614, an accreting (10^-8 M_⊙ 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 δ_dust of 10^-4 at R < 6 (Matter et al. 2016). We have detected υ = 1 →0 ¹²CO, ¹³CO, C¹⁸O, C¹⁷O, and υ = 2 →1 ¹²CO ro-vibrational emission. The lines observed are consistent with disk emission and thermal excitation. We find the following:

1. The υ = 1 →0 ¹²CO spectrum indicates that there is gas from1 AU up to 15 AU, and that there isno gap in the gas distribution. If a gap is present in the gas (i.e., aregion devoid of gas) then it should have a width smaller than2 AU.

2. The spectra of ¹³CO and C¹⁸O υ = 1 →0 emission are on average colder and emitted farther out in the disk (R> 6 AU) than the ¹²CO υ = 1 →0 emission. Keplerian flat-disk models clearly show that a drop in the gas density δ_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¹⁸O and ¹³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 δ_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_H = 3 × 10¹⁹−10²¹ cm^-2 (7 × 10^-5−2.4 × 10^-3 g cm^-2) assuming N_CO = 10^-4N_H. We derived a 5σ upper limit on the CO column density at R < 1 AU N_CO = 5 × 10¹⁵ cm^-2, which suggests an N_H < 5 × 10¹⁹ cm^-2 at R < 1 AU. The gas surface density in the disk of HD 139614 at R ≤ 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 α-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 (δ_dust < δ_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³¹).

---

Acknowledgments

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

Appendix A: Measured line centers, FWHM, integrated fluxes, and upper limits

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: (B.1)Here R is the mid-point of two grid points in the radial grid (B.2)and (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_θ): (B.4)We used 50 to 100 points in the radial direction and 1000 points in the azimuthal direction. The local line-profile φ of a grid point in R and θ 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: (B.5)The line-profile of each cell in the azimuthal direction θ is then velocity shifted to the expected local Keplerian velocity shift (B.6)thus obtaining a φ(R,θ,ν)_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: (B.7)and summing the spectra of all the annuli in the radial direction (B.8)To generate 3D channel maps the φ(R,θ,ν)_shifted of each cell is sampled in a Cartesian data cube with coordinates X,Y,ν where When the effect of the slit is taken into account, the 1D spectrum is extracted from the 3D (X,Y,ν) 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).

Appendix C: Calculation of the Bayesian probability

In each model, a χ² was calculated for each observational dataset (i.e., one χ² for the ¹²CO P(9) line-profile, one χ² for the ¹²CO rotational diagram, etc.) using (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 χ² values were calculated per model. Because the numerical value of χ² can be very different for the rotational diagram and the line-profile, before calculating the combined χ², the χ² of each observational dataset (i.e. line-profile, rotational diagram) was normalized by the minimum values of χ² of that observational dataset in the entire grid. The sum of the six normalized χ² gave the final χ² for a model. The Bayesian probability (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.

All Tables

All Figures