EDP Sciences
Free Access
Issue
A&A
Volume 598, February 2017
Article Number A118
Number of page(s) 29
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201628472
Published online 10 February 2017

© 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 disk1. 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).

Table 1

Stellar properties.

Table 2

Log of the science and calibrator observations.

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-8M/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 2. 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-4M 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, H2O, OH or CH+ in the 50200μ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.13 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.150.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 χ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 χ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 τ 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.

thumbnail Fig. 1

Composite normalized spectrum of the υ = 1 → 012CO, 13CO, and C18O lines. Error bars are 1σ in each spectrum.

Open with DEXTER

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)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-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 (RVbary = −1 × Vhelio). 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).

thumbnail Fig. 2

Examples of the υ = 1 → 012CO, 13CO, C18O, C17O and υ = 2 → 112CO 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 km s-1 are the theoretical wavelengths of each of the transitions. Note that the flux scale is larger for the υ = 1 → 012CO lines. Error bars are 3σ. Several υ = 3 → 212CO 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.

Open with DEXTER

3. Observational results

We have detected υ = 1 → 0 12CO, 13CO, C18O, C17O, and υ = 2 → 1 12CO emission lines. The υ = 3 → 2 12CO 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 → 012CO emission. To keep the notation short in the remaining of the paper, we mean by 12CO, 13CO, C18O emission υ = 1 → 012CO, 13CO, C18O 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 ~ 105) 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 12CO composite line-profile has a flat top and does not display evidence of asymmetries5. The composite 13CO and C18O lines are single peaked. Some asymmetric sub-structures are present in both lines but they are consistent with noise.

C18O emission is, at the 2σ level, 4 km s-1 narrower than 12CO emission. 13CO and 2 → 1 12CO lines are 1 km s-1 narrower than the 12CO line, at the 1σ level. To further test whether the 12CO, 13CO and C18O line-profiles are different, we ran a two-sample Kolmogorov-Smirnov (K-S) test6 on the composite 1D spectra in the ±15 km s-1 interval, after normalization by the line peaks. The K-S significance between the 12CO and C18O line-profiles is 8%, between the 12CO and 13CO line-profile is 30%, and between the C18O and 13CO line-profile is 97%. The K-S test indicates that the 12CO and C18O profiles are different (the C18O is narrower), and that statistically the 13CO profile resembles the C18O profile more than the 12CO 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 12CO 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 12CO emission extending to spatial scales larger than ~30 AU. More stringent limits are deduced in the next section.

Table 3

Types of models used to interpret the observations.

thumbnail Fig. 3

Observed composite of the normalized υ = 1 → 0 12CO 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 Rin = 1.2 AU and Rout = 15 AU (black cross in Fig. 4). Error bars on the composite line-profile are 1σ. Horizontal dotted lines are the average 1σ 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 ± 5 pc, the mean of PSF-FWHM is 27 AU.

Open with DEXTER

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 12CO 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 Rin to an outer radius Rout, where I0 is the intensity at Rin which is assumed initially to be 1. The exponent α is obtained for each pair of Rin and Rout such that I(Rout) = 0.01 × I0. 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 12CO 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.

thumbnail Fig. 4

contour plot for the grid of flat disk Keplerian models with a power-law intensity. The black cross displays the model with the lowest (0.35). The yellow and red curves show, for each Rin, the value of Rout that would generate a 1σ spectro-astrometric signal or a 1σPSF-FHWM broadening, respectively.

Open with DEXTER

We calculated a grid of disk models varying Rin between 0.1 and 70 AU and Rout between 0.2 and 100 AU. In Fig. 4 we present the contour plots of the reduced statistic (assuming three free parameters: Rin, Rout, and I0) of the model of the composite 12CO line-profile. Disk models with 0.9 <Rin< 1.8 AU, and 13 <Rout< 20 AU gave the best fit to the 12CO composite line-profile. The model that displays the smallest has Rin = 1.2 AU and Rout = 15 AU and an α exponent of the intensity −1.8.

Although disk models with Rout as large as 20 AU provide a good fit to the 12CO composite line-profile, star + disk models with a CO-emitting region with Rout> 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 Rin the value of Rout that would generate a displacement of the photocenter by 1σ). In a similar way, star + disk models with a CO-emitting region with Rout> 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 Rin the Rout 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 Rout to less than 15 AU. The model that best fits the 12CO 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 12CO composite line-profile have an inner radius around 1 AU. However, some models with smaller Rin are also compatible with the data. In Figs. B.1a,b we display the line-profiles expected for disks with Rin ranging from 0.1 to 1.2 AU. In panel (a) we show the results of the models with Rout fixed to 15 AU (α adjusted such that I(Rout) = 0.01 × I(Rin)). In panel (b) the line-profiles with α fixed to 1.8 (Rout set such that I(Rout) = 0.01 × I(Rin)). Depending on the value of α, models with Rin as low as 0.3 AU can be compatible with the observed 12CO 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 I0 at Rin. If instead we assume a soft inner edge, thus a smooth increase of the intensity from Rin up to the radius of the maximum intensity RImax = 1.2 AU, then Rin can be as small as 0.01 AU and the line profile would still be compatible with the data (see Fig. B.1c). The Rin 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 12CO line-profile shape.

thumbnail 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σ.

Open with DEXTER

4.1.2. A continuous or a gapped gas distribution?

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

The 12CO 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 12CO 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 12 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 C18O and C17O 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 12CO, 13CO and C18O line fluxes using a simple semi-infinite slab model in LTE with a single gas temperature and column density. We did not model the C17O 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 NH from 1018 to 1025 cm-2 with steps of 0.25 dex, and the gas temperature Tg 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 12CO abundance of 10-4 (N12CO = 1.0 × 10-4NH), a 12CO/13CO ratio of 100 and a 12CO/C18O ratio of 690 following Smith et al. (2009)7.

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 Rin = 0, Rout = 6.5 AU. We therefore used an average emitting region radius of 6.5 AU and a distance of 140 pc8 to determine the solid angle for the 1D slab model9.

For each NH and Tgas slab model, a rotational diagram for each CO isotopolog was calculated. We used the X and Y coordinates Y = ln(Ful/νguAul) and X = Eu. Here Ful is the line flux of the transition between the upper level u and the lower level l, gu is the degeneracy of the upper level (2J + 1), Aul is the Einstein coefficient of transition and Eu 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 (NH, Tgas, Ω), and N the number of data points (8 for 12CO, 7 for 13CO and 9 for C18O). In Fig. 6 we display the contour plots for each CO isotopolog and one combining the data of the three isotopologs.

thumbnail Fig. 6

contour plot of the modeled rotational diagrams for 12CO, 13CO, and C18O using a single temperature and surface density slab model. The combined contour plot is obtained by simultaneously using all the data of the three isotopologs. The cross indicates the location of the minimum in each panel. The numbers inside the contour indicate n-times the value of the minimum .

Open with DEXTER

thumbnail Fig. 7

Rotational diagrams and optical depths of the slab model (NH = 5 × 1021 cm-2, Tgas = 450 K) with the lowest combining all the 12CO, 13CO, and C18O 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σ.

Open with DEXTER

We find that the best fit to the rotational diagram is different for each CO isotopolog. 12CO is best described with Tgas ~ 450 K and NH> 1020 cm-2. 13CO is best reproduced by lower temperatures (Tgas ~ 380 K) and NH column densities of at least 5 × 1022 cm-2. The C18O is best fit by even lower temperatures ~350 K and NH column densities higher than 5 × 1022 cm-2. The errors on Tgas are 1020 K. The colder temperatures for 13CO and C18O emission indicate that they are produced at larger radii than the 12CO 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 NH = 5 × 1021 cm-2 and Tgas = 450 K. In this model the 12CO emission is optically thick and the 13CO and C18O optically thin (Fig. 7). The combined solution only satisfactorily describes the 12CO rotational diagram, however. The slopes of the 13CO and C18O rotational diagrams are not well reproduced. The curvature on the rotational diagrams suggests that there 13CO and C18O 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 NH, giving only lower limits for the column density.

4.2.1. Thermally excited or UV-pumped emission?

The detection of υ = 2 → 1 12CO 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 12CO/υ = 1 → 0 12CO 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 12CO emission and the upper limit of 0.04 on the υ = 3 → 2 12CO /υ = 1 → 0 12CO 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 12CO/υ = 1 → 0 12CO 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 12CO, 13CO, and C18O rotational diagrams and calculated the LTE υ = 2 → 1 12CO and υ = 3 → 2 12CO emission. We found that the observed υ = 2 → 112CO emission and the non-detection of the υ = 3 → 2 12CO emission can be well described by an LTE model, with Tgas = 460 K and NH = 1022 cm-2, in which the rotational temperature is equal to the vibrational temperature10. We show this model in Fig. 8, where we display the rotational diagrams of υ = 1 → 0 12CO and υ = 2 → 1 12CO emission of the model and the observations. This model generates υ = 3 → 2 12CO 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.

thumbnail Fig. 8

Rotational diagrams of the υ = 1 → 0 12CO (squares) and υ = 2 → 1 12CO (triangles) observed emission. Overplotted is a slab model with Trot = Tvib = 460 K and NH = 1022 cm-2, in red for υ = 1 → 0 12CO emission and in green for υ = 2 → 1 12CO emission.

Open with DEXTER

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 13CO and C18O 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 (NHRαNH and TgasRαTgas) and both distributions have no discontinuities (i.e. no gaps nor density drops or jumps in the temperature), and if the 12CO/13CO and 12CO/C18O abundance ratios are constant, then, a disk model able to reproduce the 12CO line-profile and 12CO rotational diagram would generate 13CO and C18O 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 13CO and C18O emission is optically thin up to relatively high column densities (NH ~ 1022 cm-2 and NH ~ 1023 cm-2 respectively), the 1D single power-law models generate 13CO and C18O lines that are too broad, too strong, and too warm because the column of gas at small radii is too large.

thumbnail Fig. 9

Examples of the predicted 13CO R(4) and C18O R(6) emission for a continuous power-law distribution of temperature and column density that describe the 12CO P(9) line-profile and the 12CO 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 12CO, 13CO, and C18O line-profiles and rotational diagrams.

Open with DEXTER

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 C18O and 13CO 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 C18O and 13CO 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 12CO 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 disk11. 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 Tgas and NH 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 resolution12. 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 NH(R) and Tgas(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 NH(R) and Tgas(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.

Table 4

Parameter space of the power-law NH and Tgas models.

A model has six free parameters: NH (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 12CO P(9) line at 4745.13 nm, 13CO R(4) at 4730.47 nm, and the C18O 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.

thumbnail 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 with decreasing (αinner< 0), flat (αinner = 0) and increasing (α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.

Open with DEXTER

thumbnail Fig. 11

Bayesian probability distributions of the grid of the NH and Tgas power-law models calculated (81 000 models). The χ2 statistic used to calculate the probability considers the rotational diagrams and line-profiles of the three isotopologs simultaneously.

Open with DEXTER

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 13CO and C18O lines that are too broad, strong and warm (too much warm 13CO and C18O emitted at small radii) to be consistent with the observations.

thumbnail 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σ. 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.

Open with DEXTER

thumbnail Fig. 13

Gas column density, temperature, CO optical depth, flux density, cumulative line flux, rotational diagrams, and line-profiles of the 12CO P(9), 13CO R(4), and C18O 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.

Open with DEXTER

The column density of gas traced by CO at R = 6 AU is well constrained to NH ~ 1023 cm-2. This gas column density is similar to the dust column density at R = 6 AU found by Matter et al. (2016)13. The gas column density at R = 6 AU is higher than the column density found in the single TgasNH slab model of the three CO isotopologs of Sect. 4.2. This is because higher column densities describe the C18O and 13CO rotational diagrams better. In fact, the emission of C18O and 13CO is optically thick in the 610 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 Tgas>Tdust. 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 NH = 3 × 1019 cm-2 and NH = 1021 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 Tgas>Tdust. 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 13CO and C18O 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 13CO and C18O 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 12CO P(9), 13CO R(4) and C18O R(6) line-profiles of the model in the grid that shows the best combined fit to the data. This model has an NH at R = 6 AU of 1023 cm-2, αNH inner(R< 6 AU) = + 2.0, αNH outer, δgas(R = 6 AU) = 10-2, T0(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 13CO and C18O 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 13CO and C18O line-profiles with line wings that are too strong. The effect is also seen in the 13CO and C18O 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 13CO and lines C18O disappear. The 12CO 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, 12CO becomes optically thin at R< 2 AU, the wings of the 12CO line become weaker, and the 12CO 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 12CO 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 12CO 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 NH inside the gap (at R = 6 AU) from 1021 to 1017 cm-2. The normalized theoretical 12CO P(9) spectrum was compared with the high S/N composite 12CO 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 NH inside the gap lower than 1018 cm-2 (yellow and red, δgas< 10-5) would have been seen in the 12CO composite spectrum as a line-profile with shoulders at ±15 km s-1.

Table 5

Parameters of the best-fit grid model.

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

12CO ro-vibrational emission requires relatively low column densities (NCO ~ 1015 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 NH = 5 × 1019 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 12CO 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 NCO = 5 × 1015 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 NCO > 1015 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 NH ≤ 1019 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.

thumbnail 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.

Open with DEXTER

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 physically14. 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 AU15.

thumbnail Fig. 15

12CO 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 12CO 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 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σ.

Open with DEXTER

C18O emission is detected at lower S/N than 12CO and 13CO emission. To further test the reliability of the gas density drop, we recalculated the Bayesian probability diagrams only taking into account the 12CO and 13CO 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 12CO/H2, 12CO/13CO and 12CO/C18O 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 13CO and C18O requires higher column densities than for 12CO. 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 magnitude16. 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 13CO and C18O lines (see Fig. 13) are optically thin17. Therefore, at R < 6 AU, the 13CO and C18O 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 13CO and C18O 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 13CO and C18O 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-plane18. The 12CO 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 12CO, 13CO and C18O 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 temperature19. The column of C18O (thus NH) 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 NH = 1021 cm-2. Higher NH, even with T0 (R = 1 AU) = 575 K (100 K lower than the best-fit grid model), would generate 13CO and C18O lines with high-velocity wings which would be too strong to be compatible with the CRIRES spectrum.

thumbnail Fig. 16

Predicted 12CO P(9) line-profile for a disk extending from 0.1 to 1.0 AU with a flat surface density NH = 5 × 1019 cm-2 and the extrapolated temperature profile from the best-fit grid model. Error bars are 3σ. The horizontal dashed line is the 5σ limit.

Open with DEXTER

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 12CO line-profile, and the fact that the average temperatures of 13CO and C18O 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 13CO and C18O 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.

thumbnail 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 αNH inner ≤ + 1.0. Note that the value of δ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.

Open with DEXTER

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 likely21.

The model grid was calculated assuming Rin = 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 12CO 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 12CO lines. We checked the predicted υ = 2 → 1 12CO 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 NH at 6 AU three times larger (NH (R = 6 AU) = 3 × 1023 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. C18O lines are a few km s-1 narrower than 12CO lines.

  3. 13CO and C18O emission are 50–100 K colder than 12CO 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 NH at 1 < R < 6 AU ranges between 3 × 1019 and 1021 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 NCO at R ≤ 1 is 5 × 1015 cm-2, which corresponds to a gas column density NH < 5 × 1019 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 hundred22. 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 12CO 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 12CO 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 disks23. This difference arises because in primordial disks the 12CO 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 12CO ro-vibrational emission becomes optically thick at relatively low column densities (NCO ~ 1015 cm-2), the lack of high-velocity wings in the 12CO 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 NH ~ 5 × 1019 cm-2 on the column of gas at R < 1 AU (assuming standard abundances) derived for HD 139614 is significantly lower than the NH = 1024−1026 cm-2 typical of primordial disks.

In HD 139614, the 13CO emission is narrower and colder than the 12CO 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 12CO and 13CO 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 (Tex = Trot = Tvib). 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 sources24. 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-8M/yr in HD 139614 (Garcia Lopez et al. 2006) is at the high end of the mass-loss rates (10-8−10-10M/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 (Rc,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 (105 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 × 1029 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-10M/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 56 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 observations25. 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 AU26; (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 (RH) of the planet. The Hill radius of a planet of mass Mp at a distance rp 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 MJ 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 MJ 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 MJ 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 MJ27.

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-simulations28.

Moreover, a 12 MJ 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 MJ 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 12CO 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 MJ, which is consistent with the upper limit of 2 MJ 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 MJ 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 aopt = 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 NH = 1019−1021 cm-2 (2.4 × 10-5−2.4 × 10-3 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 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-8M 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-9M 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 vR of the gas at radius R through disk = −2πRvR(R) Σ(R). If we now assume that disk ~ , then at a radius of 1 AU and for an accretion rate of 10-8M yr-1, we find that | vR | 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 (vK) 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 |vR| ~ αh2vK, 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 |vR| can reach at most ~10-5vK. 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-8M yr-1 would suggest an NH of 1025−1026 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 NH ~ 1025 cm-2 at R < 1 AU in a disk without dust in the cavity (best case for UV-photodissociation) around a Teff = 10 000 K central star (i.e., brighter in UV than HD 139614) NCO would be 1021 cm-2 at R < 1 AU (see Fig. 8 in Bruderer 2013), a value much higher than the 5σ limit of NCO = 5 × 1015 cm-2 from our CRIRES data. According to the Bruderer (2013) models (for a Teff = 10 000 K), to have NCO < 1015 cm-2, NH needs to below 1022 cm-2 at R < 1 AU for an inner disk without dust. For HD 139614, which has a Teff ~ 7800 K (and lower UV, accordingly) and some dust in the inner disk30, NH 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 12CO 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 12CO 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-8M 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 12CO, 13CO, C18O, C17O, and υ = 2 →1 12CO ro-vibrational emission. The lines observed are consistent with disk emission and thermal excitation. We find the following:

  1. The υ = 1 →0 12CO 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 13CO and C18O υ = 1 →0 emission are on average colder and emitted farther out in the disk (R> 6 AU) than the 12CO υ = 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 C18O and 13CO 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 NH = 3 × 1019−1021 cm-2 (7 × 10-5−2.4 × 10-3 g cm-2) assuming NCO = 10-4NH. We derived a 5σ upper limit on the CO column density at R < 1 AU NCO = 5 × 1015 cm-2, which suggests an NH < 5 × 1019 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 ALMA31).


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.

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.

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.

5

We note that the 1σ error in the flux in the 12CO 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.

6

KSTWO function in IDL.

7

The 12CO/13CO and 12CO/C18O 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 12C/13C from Smith et al. is nearly twice the expected interstellar medium (ISM) ratio, and the 12CO/C18O ratio is ×1.4 the ISM ratio. We used the 12CO/13CO and 12CO/C18O ratios of Smith et al. instead of the ISM ratios because we consider that they are more representative of the isotopologs ratios that would be expected in the inner disk of HD 139614.

8

Models were calculated before the recent Gaia distance measurement of 131 ± 5 pc, the difference in distance of 59 pc does not change the conclusions reached.

9

We note that the Rout 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.

10

The model with Tgas = 450 K and NH = 1022 cm-2 of Fig. 7 generated slightly weaker υ = 2 → 1 12CO emission.

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.

12

Including the instrument resolution in the convolution kernel of the thermal and turbulent broadening enable saving hundreds of convolutions per model, 107 convolutions in the grid, and it is equivalent to convolving the final data-cube by the instrument resolution.

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.

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.

15

We note, however, that models with flat surface density profiles produce a poorer fit to the observations than models with positive αNH inner because of the stronger 12CO high-velocity wings (Fig. 12).

16

For an application of selective photodissociation to the Solar Nebula see Lyons & Young (2005).

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.

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.

19

The best-fit grid model has a radial temperature that is a compromise between the 12CO that requires a higher temperature (emitted higher up) and 13CO and C18O emission that require lower temperatures (emitted farther down). The best-fit model combining 12CO, 13CO, and C18O has a gas T0 (R = 1 AU) = 675 K. This temperature is in between the temperature of the best-fit grid model for 12CO alone (T0 (R = 1 AU) = 725 K) and the temperature of the best-fit grid model for C18O emission alone (T0 (R = 1 AU) = 625 K).

20

That dust is seen at R ~ 6 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.

21

The change in the gas-density drop radius can be compensated for by changes in δgas or the surface density or temperature exponent.

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 × 10-17 W m-2, 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.

23

The 12CO ro-vibrational emission from primordial disks is two to five times broader than that of transition disks (Brown et al. 2013).

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.

25

We note, however, that the X-ray photoevaporation models described here were mostly developed for 1 M stars and may not be directly applicable to A-type stars.

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.

27

We note that the planet mass depends on assumptions such as the disk viscosity. Models of Matter et al. (2016) assumed α = 0.006.

28

For a study of the interplay of photoevaporation and planet formation, see, for example, Rosotti et al. (2013).

29

Here the hydrogen mass fraction of gas with solar composition is assumed (0.7) to calculate the gas surface density in g cm-2.

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.

31

At the time of writing, HD 139614 has not been observed 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

Table A.1

CO ro-vibrational line fluxes, upper limits, and average line ratios.

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).

thumbnail 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 α fixed. In the models, α or Rout are adjusted such that I(Rout) = 0.01 × I(Rin); 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 α = −1.8 at R> 1.2 AU. Error bars on the composite spectrum are 1σ.

Open with DEXTER

Appendix C: Calculation of the Bayesian probability

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

thumbnail Fig. C.1

Bayesian probability plots for a sub-sample of the grid in which only the models with αNH inner ≤ + 1 are considered (54 000 models).

Open with DEXTER

thumbnail Fig. C.2

Surface density, temperature, CO optical depth, flux density, cumulative line flux, rotational diagrams, line-profiles of the 12CO P(9), 13CO R(4), and C18O R(6) emission for best-fit grid model when αNH inner ≤ + 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σ 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.

Open with DEXTER

thumbnail Fig. C.3

Bayesian probability plots for the grid using only the 12CO and 13CO data (i.e. no C18O 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).

Open with DEXTER

thumbnail 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 (Rgap rmax) between 4.0 and 6.0 AU. Rgap rmax 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.

Open with DEXTER

thumbnail Fig. C.5

Flux density, cumulative line flux, rotational diagrams, and the line-profiles of the 12CO P(9), 13CO R(4), and C18O 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.

Open with DEXTER

thumbnail Fig. C.6

Surface density and temperature profile, and predicted υ = 1 → 012CO P(9), 13CO R(4), C18O R(6), and υ = 2 → 112CO P(3) and 12CO P(4) line-profiles for a model with NH at R = 6 AU three times larger than the best-fit grid model (NH (R = 6 AU) = 3 × 1023 cm-2). The higher NH enables to describe the υ = 2 → 112CO P(3) and 12CO P(4) line-profiles while keeping a good fit to the υ = 1 → 0 lines. The model has the same temperature structure and same surface density at R< 6 as the best-fit grid model (thus δgas = 3.3 × 10-3). The cumulative flux plot shows that the υ = 2 → 1 lines are dominated by the contribution at 6 <R< 10 AU. Errors in the plot are 3σ, and the dashed horizontal line is the 5σ limit.

Open with DEXTER

Table C.1

Transition disks with a quantified gas density drop inside the dust cavity.

All Tables

Table 1

Stellar properties.

Table 2

Log of the science and calibrator observations.

Table 3

Types of models used to interpret the observations.

Table 4

Parameter space of the power-law NH and Tgas models.

Table 5

Parameters of the best-fit grid model.

Table A.1

CO ro-vibrational line fluxes, upper limits, and average line ratios.

Table C.1

Transition disks with a quantified gas density drop inside the dust cavity.

All Figures

thumbnail Fig. 1

Composite normalized spectrum of the υ = 1 → 012CO, 13CO, and C18O lines. Error bars are 1σ in each spectrum.

Open with DEXTER
In the text
thumbnail Fig. 2

Examples of the υ = 1 → 012CO, 13CO, C18O, C17O and υ = 2 → 112CO 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 km s-1 are the theoretical wavelengths of each of the transitions. Note that the flux scale is larger for the υ = 1 → 012CO lines. Error bars are 3σ. Several υ = 3 → 212CO 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.

Open with DEXTER
In the text
thumbnail Fig. 3

Observed composite of the normalized υ = 1 → 0 12CO 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 Rin = 1.2 AU and Rout = 15 AU (black cross in Fig. 4). Error bars on the composite line-profile are 1σ. Horizontal dotted lines are the average 1σ 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 ± 5 pc, the mean of PSF-FWHM is 27 AU.

Open with DEXTER
In the text
thumbnail Fig. 4

contour plot for the grid of flat disk Keplerian models with a power-law intensity. The black cross displays the model with the lowest (0.35). The yellow and red curves show, for each Rin, the value of Rout that would generate a 1σ spectro-astrometric signal or a 1σPSF-FHWM broadening, respectively.

Open with DEXTER
In the text
thumbnail 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σ.

Open with DEXTER
In the text
thumbnail Fig. 6

contour plot of the modeled rotational diagrams for 12CO, 13CO, and C18O using a single temperature and surface density slab model. The combined contour plot is obtained by simultaneously using all the data of the three isotopologs. The cross indicates the location of the minimum in each panel. The numbers inside the contour indicate n-times the value of the minimum .

Open with DEXTER
In the text
thumbnail Fig. 7

Rotational diagrams and optical depths of the slab model (NH = 5 × 1021 cm-2, Tgas = 450 K) with the lowest combining all the 12CO, 13CO, and C18O 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σ.

Open with DEXTER
In the text
thumbnail Fig. 8

Rotational diagrams of the υ = 1 → 0 12CO (squares) and υ = 2 → 1 12CO (triangles) observed emission. Overplotted is a slab model with Trot = Tvib = 460 K and NH = 1022 cm-2, in red for υ = 1 → 0 12CO emission and in green for υ = 2 → 1 12CO emission.

Open with DEXTER
In the text
thumbnail Fig. 9

Examples of the predicted 13CO R(4) and C18O R(6) emission for a continuous power-law distribution of temperature and column density that describe the 12CO P(9) line-profile and the 12CO 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 12CO, 13CO, and C18O line-profiles and rotational diagrams.

Open with DEXTER
In the text
thumbnail 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 with decreasing (αinner< 0), flat (αinner = 0) and increasing (α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.

Open with DEXTER
In the text
thumbnail Fig. 11

Bayesian probability distributions of the grid of the NH and Tgas power-law models calculated (81 000 models). The χ2 statistic used to calculate the probability considers the rotational diagrams and line-profiles of the three isotopologs simultaneously.

Open with DEXTER
In the text
thumbnail 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σ. 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.

Open with DEXTER
In the text
thumbnail Fig. 13

Gas column density, temperature, CO optical depth, flux density, cumulative line flux, rotational diagrams, and line-profiles of the 12CO P(9), 13CO R(4), and C18O 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 15

12CO 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 12CO 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 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σ.

Open with DEXTER
In the text
thumbnail Fig. 16

Predicted 12CO P(9) line-profile for a disk extending from 0.1 to 1.0 AU with a flat surface density NH = 5 × 1019 cm-2 and the extrapolated temperature profile from the best-fit grid model. Error bars are 3σ. The horizontal dashed line is the 5σ limit.

Open with DEXTER
In the text
thumbnail 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 αNH inner ≤ + 1.0. Note that the value of δ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.

Open with DEXTER
In the text
thumbnail 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 α fixed. In the models, α or Rout are adjusted such that I(Rout) = 0.01 × I(Rin); 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 α = −1.8 at R> 1.2 AU. Error bars on the composite spectrum are 1σ.

Open with DEXTER
In the text
thumbnail Fig. C.1

Bayesian probability plots for a sub-sample of the grid in which only the models with αNH inner ≤ + 1 are considered (54 000 models).

Open with DEXTER
In the text
thumbnail Fig. C.2

Surface density, temperature, CO optical depth, flux density, cumulative line flux, rotational diagrams, line-profiles of the 12CO P(9), 13CO R(4), and C18O R(6) emission for best-fit grid model when αNH inner ≤ + 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σ 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.

Open with DEXTER
In the text
thumbnail Fig. C.3

Bayesian probability plots for the grid using only the 12CO and 13CO data (i.e. no C18O 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).

Open with DEXTER
In the text
thumbnail 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 (Rgap rmax) between 4.0 and 6.0 AU. Rgap rmax 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.

Open with DEXTER
In the text
thumbnail Fig. C.5

Flux density, cumulative line flux, rotational diagrams, and the line-profiles of the 12CO P(9), 13CO R(4), and C18O 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.

Open with DEXTER
In the text
thumbnail Fig. C.6

Surface density and temperature profile, and predicted υ = 1 → 012CO P(9), 13CO R(4), C18O R(6), and υ = 2 → 112CO P(3) and 12CO P(4) line-profiles for a model with NH at R = 6 AU three times larger than the best-fit grid model (NH (R = 6 AU) = 3 × 1023 cm-2). The higher NH enables to describe the υ = 2 → 112CO P(3) and 12CO P(4) line-profiles while keeping a good fit to the υ = 1 → 0 lines. The model has the same temperature structure and same surface density at R< 6 as the best-fit grid model (thus δgas = 3.3 × 10-3). The cumulative flux plot shows that the υ = 2 → 1 lines are dominated by the contribution at 6 <R< 10 AU. Errors in the plot are 3σ, and the dashed horizontal line is the 5σ limit.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.