© ESO, 2016

1. Introduction

Viscous accretion, photoevaporation, and dynamical clearing are processes through which disks are thought to shape and dissipate most of their mass, which lies in the gas content (see Williams & Cieza 2011; Baruteau et al. 2014; Alexander et al. 2014,for a review). Dust also plays an important role since it dominates the disk opacity and provides the raw material for building the rocky planets and the giant-planet cores. Dust is affected by processes stemming from its coupling to the gas such as dust trapping, radial migration, and grain growth (e.g., Birnstiel et al. 2012; Laibe 2014), or from the stellar irradiation (e.g., Dominik & Dullemond 2011). Such evolution processes can produce specific signatures such as gaps, inner holes, and asymmetries (see, e.g., Crida & Morbidelli 2007; Owen et al. 2011; Meheut et al. 2012). Identifying them in the inner disk regions (~0.1–10 au) is essential since these regions are the expected cradle of telluric planets (Righter & OBrien 2011) and the location of the photoevaporation onset (Alexander et al. 2014). The emission deficit in the infrared (IR) SED of pre-transitional and transitional disks (e.g., Strom et al. 1989; Cieza et al. 2010; Merín et al. 2010) has been commonly interpreted as a clearing of their inner regions, which make these objects relevant laboratories for observing the signatures of disk evolution and dissipation processes. Infrared interferometry can specifically probe these inner disk regions (Dullemond & Monnier 2010; Carmona et al. 2014) and detect dust radial evolution (van Boekel et al. 2004), brigthness asymmetries (e.g., Kraus et al. 2013), and dust clearing (e.g., Benisty et al. 2010; Menu et al. 2014).

The Herbig star HD 139614 (see Table 1) is of particular interest here. It is associated with the Upper Centaurus Lupus (UCL) region of the Sco OB2 association located at 140 ± 27 pc (Preibisch & Mamajek 2008,and references therein). Its group-I SED (Meeus et al. 2001) presents pretransitional features (near-infrared (NIR) excess with mid-infrared (MIR) emission deficit at 6–7 μm) and weak MIR amorphous silicate features (Juhász et al. 2010), which suggests significant dust evolution in the inner regions. High-resolution imaging with T-RECS/Gemini only resolved the MIR continuum emission at 18 μm (17 ± 4 au, Mariñas et al. 2011), while MIR interferometry probed the inner N-band emitting region and revealed a narrow gap-like structure (located from ~2.5 to 6 au) depleted in warm μm-sized dust (Matter et al. 2014). Recent works showed or found hints of the presence of gaps in Group I disks (e.g., Maaskant et al. 2013; Menu et al. 2015), and HD 139614 is a rare, if not unique, case of object for which a small au-sized gap (~3 au) has been spatially resolved. Previous studies, mainly based on sub-mm interferometry (e.g., Andrews et al. 2011; Pérez et al. 2014) and NIR imaging (Quanz et al. 2013; Avenhaus et al. 2014), have focused on objects with large cavities (~10–100 au). Their origin is ambiguous and possibly combines, for example, photoevaporation (Rosotti et al. 2013), magnetorotational instability (Chiang & Murray-Clay 2007), and/or multiple unseen planets (Zhu et al. 2011). Single Jovian planets are expected to open au-sized gaps while inducing a gas and dust surface density decrease in the inner disk regions (e.g., Crida & Morbidelli 2007). Knowing that spectroscopic observations have proven the presence of gas (Panić & Hogerheijde 2009; Meeus et al. 2012) and gas tracers like PAHs (Acke et al. 2010) in the HD 139614 disk, this object constitutes a unique opportunity to characterize the early stages of inner disk dispersal and potentially witness planet-disk interaction. Spatially resolved IR observations are thus required. We report the first multiwavelength analysis of HD 139614, combining MIR VLTI/MIDI data with new NIR interferometric data obtained with VLTI/PIONIER (Le Bouquin et al. 2011) and VLTI/AMBER (Petrov et al. 2007), and far-IR photometry with Herschel¹/PACS (Pilbratt et al. 2010; Poglitsch et al. 2010). We aim to obtain new and robust constraints on the innermost ~0.1–1 au dust structure. We also aim to use radiative transfer to refine the constraints of the previous analytical modeling of the MIDI data (e.g., gap characteristics, outer disk properties). This is essential to determining the degree of differentiation of the inner and outer disk regions, on either side of the gap, and identify viable mechanisms for the inner disk clearing around HD 139614.

Section 2 summarizes the new observations that complement the previous VLTI/MIDI and optical/IR SED data used in Matter et al. (2014). Section 3 presents the analysis and geometrical modeling of the PIONIER and AMBER data. Section 4 describes the radiative transfer modeling of the broadband SED and the complete set of interferometric data. Section 5 discusses the modeling results against the mechanisms possibly responsible for the gap structure. This includes a comparative study with hydrodynamical simulations of gap opening by a planet. Finally, Sect. 6 summarizes our work and outlines the perspectives.

2. Observations and data processing

2.1. VLTI/PIONIER observations

We observed HD 139614 with the 1.8 m Auxiliary Telescopes (AT) in the frame of a large program on Herbig stars (ID: 190.C-0963) conducted with VLTI/PIONIER (see Table 2). It combines the light from four telescopes in H band, and provides six squared visibilities, noted as V², and four closure phases per ATs configuration. Our observations were obtained at low spectral resolution (R ~ 40) providing three spectral channels centered at 1.58 μm, 1.67 μm, and 1.76 μm. The calibration stars were selected with the SearchCal tool from the Jean-Marie Mariotti Center (JMMC). Each HD 139614 observation was bracketed by two calibrator observations, and the estimated transfer function was interpolated at the time of observation. The projected baseline lengths range from 10.3 m to 139.8 m (λ/ 2B ≃ 17 mas to 1 mas). The PIONIER field of view (FOV) is ~200 mas (~25 au at 140 pc). We performed a standard data reduction using the pipeline “pndrs”, described in Le Bouquin et al. (2011). Figure 1 shows the UV coverage and calibrated data. The final uncertainties include the statistical errors and the uncertainty on the transfer function.

2.2. VLTI/AMBER observations

HD 139614 was observed with AMBER (program 0.89.C-0456(A)) using the 8-m telescope triplet UT1-UT2-UT4. AMBER can combine three beams and provides spectrally dispersed V² and closure phases. Our observations were obtained in low spectral resolution (R ~ 35) in the H and K bands and included two interferometric calibrators (see Table 2) chosen using the SearchCal tool. The data consisted of two sets of eight and five exposures of 1000 frames. The data were reduced with the JMMC “amdlib” package (release 3.0.5). For each exposure, a frame selection was made to minimize the impact of the instrumental jitter and the non-optimal light injection into the optical fibers. Twenty percent of the frames with the highest fringe S/N provided the smallest errors on the resulting V². Despite this, the H-band data showed significant variability and low S/N (~1.5) probably because of the H = 7.3 magnitude of HD 139614 that was close to the AMBER limiting magnitude (H_corr = 7.5). Moreover, the low-resolution closure phases are affected by a strong dependency on the piston, which is not reduced by stacking frames (see AMBER manual, ESO doc. VLT-MAN-ESO-15830-3522). Therefore, from the AMBER observations, we only kept the K-band V² measurements. The instrumental transfer function was calibrated using the closest calibrator in time. The calibrated V² errors include the statistical error obtained when averaging the individual frames and the standard deviation of the transfer function over the calibrator observations. Our final dataset consists of six dispersed V² in the [2.0–2.5] μm range. The projected baseline lengths range from 52.0 m to 127.7 m (5 mas to 1.8 mas at λ = 2.2μm). The AMBER FOV is ~60 mas (~9 au at 140 pc). Figure 1 shows the UV coverage and calibrated data.

2.3. Photometric data

We complemented the SED used in Matter et al. (2014) with Herschel/PACS observations acquired on 7 March 2011 in the frame of the Key Program GASPS (Dent et al. 2013). Data reduction, flux extraction, and error estimation were performed as in Olofsson et al. (2013), and lead to (9.66 ± 0.96) × 10^-13 W.m^-2 at 70 μm, (4.93 ± 0.01) × 10^-13 W.m^-2 at 100 μm, and (2.43 ± 0.01) × 10^-13 W.m^-2 at 160 μm. We also included measurements at 800 μm and 1.1 mm taken from Sylvester et al. (1996). Except for the recent Herschel data, we assumed a 10% relative uncertainty on the SED points, which is conservative, especially for the 2MASS data with formal uncertainties ~5% (Skrutskie et al. 2006). The broadband SED is shown in Fig. 3.

2.4. Stellar parameters and spectrum

Table 1 shows the stellar parameters used for HD 139614. Following Preibisch & Mamajek (2008), we adopt a distance of 140 ± 27 pc hereafter. For the stellar flux, we use the same Kurucz spectrum (T_eff = 7750 K, log g = 4.0, Fe / H = −0.5) as in Matter et al. (2014). Given the temperature quoted in Table 1 (T_eff ~ 7850 K), this may be slightly too cold but remains consistent with other temperature estimates (e.g., Alecian et al. 2013).

3. Observational analysis and geometrical modeling

3.1. Observational results

As shown in Fig. 1, the PIONIER (top) and AMBER (bottom) V² shows circumstellar emission that is well resolved at a level of a few mas. The PIONIER V² data present an exponential-like profile with a steep decrease at low frequencies (≤10 m/μm) that may suggest a fully resolved emission. The PIONIER FOV (~25 au) partly encompasses the outer disk, which starts at ~6 au (Matter et al. 2014). NIR scattering by sub- and micron-sized grains at the outer disk’s surface may thus contribute to this steep decrease. At high frequencies (≥40 m/μm), the V² reach an asymptotic level between 0.5 at 1.58 μm and 0.35 at 1.76 μm, which translates to a visibility level between 0.7 and 0.6, respectively. This is close to the stellar-to-total flux ratio (STFR) evaluated to 0.7 at 1.65 μm using the stellar Kurucz spectrum of Sect. 2.5 and the NIR SED (2MASS measurements). This indicates that the circumstellar emission is fully resolved by PIONIER at high spatial frequencies. In the K band, the AMBER V² measurements range from 0.1 at 2 μm to 0.3 at 2.5 μm for the longest baselines (UT2-UT4 and UT1-UT4) and from 0.15 to 0.35 for the shortest one (UT1-UT2). The corresponding average V level is 0.4 (UT2-UT4 and UT1-UT4) and 0.5 (UT1-UT2) across the K band, which is close to the STFR evaluated to 0.4 at 2.2 μm. The inner disk is thus fully resolved in the K band, at least for the longest baselines. This suggests a NIR-emitting region that probably spreads over at least one au (~7 mas at 140 pc), as already suggested by Matter et al. (2014).

Our observations also show that V² depends on wavelength across both spectral bands. Such a chromatic effect has been studied, for instance, in the frame of IR image reconstruction (Kluska et al. 2014). At high spatial frequencies, where the disk is fully resolved, the chromatic variation in the asymptotic level of visibility is ~0.1 (expressed in V) across the H band and ~0.2 across the K band. This reveals the decrease in the unresolved stellar contribution combined with a positive chromatic slope due to the emission from the hottest dust grains. Assuming a blackbody law for this hottest component, we could reproduce the NIR SED (from λ = 1.25μm to λ = 3.4μm) and the STFR chromaticity across the H band (~0.1) and the K band (~0.2), with dust grain temperatures from 1300 to 1500 K, which is close to the sublimation temperature for silicates. We note that the hottest emission at 1500 K induces a slightly lower STFR (~0.65) at 1.65 μm than the ratio induced by the blackbody emission at 1300 K (~0.7). The latter is closer to the STFR derived from the 2MASS measurement at 1.65 μm, namely 0.7.

The PIONIER closure phases do not show a clear departure from zero, hence noticeable signatures of brightness asymmetries. We do note a dispersion of ≲2° for the yellow and red measurements and ≲5° for the noisier blue and green ones (see Fig. 1). Since the closure phase produced by a binary system is, at the first order, proportional to the flux ratio between the components (Vannier et al. 2006), a closure phase dispersion of ±2–5° would translate to an upper limit of ~3–8 × 10^-2 on the flux ratio. A flux ratio of ~10^-2 would imply an upper mass limit for the hypothetical companion of about 0.11 M_⊙, using the recent BT-SETTL atmospheric models for very low-mass stars, brown dwarfs, and exoplanets of Allard et al. (2012). Based on the available PIONIER data, it is thus unlikely that HD 139614 is actually a tight binary system hosting a companion that is more massive than 0.11 M_⊙ in its nearby environment (~au).

3.2. Geometrical modeling

We used a geometrical approach to derive the basic characteristics of the NIR-emitting region. Given the exponential-like profile of PIONIER V² measurements, we considered a Lorentzian-like brightness distribution to represent the NIR emission, independently in H and K bands. This centrally peaked brightness profile has broader tails than the usual Gaussian profile. It can be used to estimate a characteristic size for a spatially resolved emitting region that presents a gradually decreasing brightness profile with smooth outer limits. This is relevant for HD 139614, given its expected spatially extended NIR-emitting region. The squared visibility of a Lorentzian-like circumstellar emission can be written as $\begin{array}{ccc}{\mathit{V}}_{\mathrm{circ}}^{\mathrm{2}}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}{\left|\exp \left(\mathrm{-}\mathrm{2}\mathit{\pi }\frac{{\mathit{r}}_{\mathrm{50}\mathrm{\%}}}{\sqrt{\mathrm{3}}}\frac{{\mathit{B}}_{\mathrm{eff}}\mathrm{\left(}\mathit{i,}\mathrm{PA}\mathrm{\right)}}{\mathit{\lambda }}\right)\right|}^{\mathrm{2}}\mathit{,}& & \end{array}$(1)where r_50% is defined as the angular radius of half-integrated flux (containing 50% of the total flux), and B_eff(i,PA) is the effective baseline (see, e.g., Matter et al. 2014) with i and PA the inclination and position angle of the circumstellar component. Then, the total visibility ${\mathit{V}}_{\mathrm{tot}}^{\mathrm{2}}$ is ${\mathit{V}}_{\mathrm{tot}}^{\mathrm{2}}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}\mathrm{\left|}{\mathit{f}}_{\mathrm{\ast }}\mathrm{\right(}\mathit{\lambda }\mathrm{\left)}{\mathit{V}}_{\mathrm{\ast }}\mathrm{\right(}\mathit{u,v}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{\left)}{\mathit{V}}_{\mathrm{circ}}\mathrm{\right(}\mathit{u,v}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}$, with V_∗ = 1 the visibility of the unresolved star. Here, f_∗(λ) is the star-to-total flux ratio that is estimated at each wavelength using the Kurucz spectrum of Sect. 2.5 for the star and a single-temperature blackbody for the circumstellar contribution. To limit the number of free parameters, we only consider two STFRs f_{∗ ,1300} and f_{∗ ,1500}, calculated with a 1300 K and a 1500 K blackbody emission (see Sect. 3.1). The free parameters of the model are r_50% (in mas), i, and PA.

3.2.1. PIONIER

Considering the complete set of V² data, we computed a grid of models by scanning the parameters space in 70 steps and, from the χ² calculated from the measured and modeled V², derived the Bayesian probability (exp[−χ²/ 2]) for each of the parameters. These marginal probability distributions are the projection of exp[−χ²/ 2] along the three dimensions of the parameters space. The range of explored values is [1–20] mas (0.1–2.8 au at 140 pc) for r_50%, [0°–70°] for i, and [0°–175°] for PA. It appears that low i values (≲35°) are favored, while no clear constraint is obtained on PA. Indeed, with i ≲ 35°, all PA values between 0° and 150° fit the data almost equally well (${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\simeq }\mathrm{2.3}$–2.6). Since our available PIONIER dataset does not seem to highlight a significant difference in position angle and inclination relative to the outer disk, we adopted the values derived by Matter et al. (2014) for the outer disk: 112 ± 9° and 20 ± 2° (with 1σ uncertainties), respectively. We again explored the same range of values for r_50% for the two STFRs. The best-fit solution is represented by r_50% = 3.9 ± 0.1 mas (0.55 ± 0.01 au at 140 pc) for f_{∗ ,1300}, and by r_50% = 2.7 ± 0.2 mas (0.4 ± 0.02 au at 140 pc) for f_{∗ ,1500}. The 1σ uncertainties correspond to the 68% confidence interval derived directly from the marginal probability distribution. Using the f_{∗ ,1500} ratio (f_{∗ ,1500}<f_{∗ ,1300} in H band) leads to a better fit (${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{=}\mathrm{2.6}$) than the f_{∗ ,1300} ratio (${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{=}\mathrm{6.6}$). As shown in Fig. 2, the V² decrease at low spatial frequencies is reproduced in both cases, while the modeled V² for f_{∗ ,1300} significantly overestimates the V² measurements at the highest spatial frequencies. This suggests a discrepancy between the STFR predicted by the asymptotic V² level and the one estimated from the 2MASS measurements and the Kurucz spectrum. Possible causes are 1) the uncertainties on the 2MASS measurements (~5% Skrutskie et al. 2006) and on the stellar parameters of the Kurucz model; 2) a change in the STFR between the time of the 2MASS and PIONIER observations, which however, appears unlikely given the absence of significant visible or MIR variability (Meeus et al. 1998; Kóspál et al. 2012) and the face-on orientation of the disk (see Matter et al. 2014,and references therein); and 3) the degradation of V² measurements at long baselines, where the coherent flux is lower. HD 139614 has a H mag = 7.3 that is close to the PIONIER sensitivity limit (H ~ 7, Le Bouquin et al. 2011).

Finally, no fully resolved emission was needed to reproduce the steep V² decrease at low spatial frequencies.

3.2.2. AMBER

Considering the complete set of AMBER dispersed V² data, we followed the same procedure as for PIONIER (same parameter space and same derivation of the formal errors). It quickly appeared that our AMBER interferometric dataset is too sparse in UV coverage (one baseline triplet covering only projected baselines lengths longer than 50 m) to constrain the inclination and the position angle of the K-band emitting region. Therefore, assuming again that the inner component is coplanar with the outer disk, we explored a broad range of r_50% values (1–20 mas) and found a best-fit solution represented by r_50% = 4.3 ± 0.1 mas (0.60 ± 0.02 au at 140 pc) for f_{∗ ,1300}, and by r_50% = 4.2 ± 0.1 mas (0.58 ± 0.02 au at 140 pc) for f_{∗ ,1500}. As shown in Fig. 2, both solutions agree well with the measured V² (${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{~}\mathrm{1.6}$), since f_{∗ ,1300}(λ) ≃ f_{∗ ,1500}(λ) across the K-band.

Our new NIR data spatially resolved the innermost region of HD 139614 and suggest that hot dust material is located around the expected dust sublimation radius of micron-sized silicate dust (≃0.2 au). We do not rule out the presence of refractory dust material within this radius. Moreover, the inner dust component probably extends out to 1 au or 2 au. Indeed, our derived r_50% values imply that, for instance, 80% of the flux in such Lorentzian profiles would then be contained within ~1.1 au in the H band, and ~1.7 au in the K band. Considering the uncertainty on the stellar distance (±27 pc), the r_50% values would vary from 0.35 au to 0.45 au in the H band and from 0.5 au to 0.7 au in the K band. These values are still consistent with our conclusions. The slight difference in size estimate between the H and K bands suggests a disk chromaticity that is probably related to a temperature gradient. Finally, the data do not allow us to constrain a difference in inclination and position angle between the inner and outer components and are consistent with a coplanarity.

4. Radiative transfer modeling

4.1. Disk model

Based on Matter et al. (2014), we consider a two-component model composed of an inner and an outer dust component spatially separated by a dust-depleted region. Although we first assume an empty gap when performing the model’s grid computation and χ² minimization, an upper limit on the dust mass in the gap will then be estimated. We mention that the results of Matter et al. (2014), which showed an incompatibility of the IR SED and the mid-IR visibilities with a continuous disk structure, confirmed previous results based on SED modeling. Notably, the possibility of a partially self-shadowed continuous disk was explored by Dominik et al. (2003). Using a passive disk model with a puffed-up inner rim, they could not reproduce both the NIR excess and the rising MIR spectrum of HD 139614 (see their Fig.4). The artificial increase in the inner rim scale height, which is required to match the NIR excess, implied a strong shadowing of the outer disk. This led to decreasing and too low MIR emission. We performed radiative transfer modeling tests with a continuous disk, including a puffed-up and optically thick inner rim to induce self-shadowing over the first few au. However, this model systematically produced too much flux in the 5–8 μm region, and decreasing and too weak emission at λ> 8μm. This implied MIR visibilities without sine-like modulation, in contrast with the MIDI data. Moreover, the puffed-up inner rim systematically induced a too spatially confined NIR-emitting region that is incompatible with the NIR interferometric data.

4.1.1. RADMC3D

We used the radiative transfer code RADMC3D to produce disk images and SED (Dullemond & Dominik 2004). Its robustness and accuracy were validated through benchmark studies (e.g., Pinte et al. 2009). This code can compute the dust temperature distribution using the Monte-Carlo method of Bjorkman & Wood (2001) with improvements like the continuous absorption method of Lucy (1999). We considered an axially symmetric two-dimensional disk in a polar coordinate system (r, θ) with a logarithmic grid spacing in r and θ. An additional grid refinement in r was applied to the inner edge of both components to ensure that the first grid cell is optically thin. The radiation field and temperature structure computed by the Monte Carlo runs were used to produce SEDs and images by integrating the radiative transfer equation along rays (ray-tracing method). Isotropic scattering was included in the modeling. While the thermal source function is known from a first Monte Carlo run, the scattering source function is computed at each wavelength through an additional Monte Carlo run prior the ray-tracing.

4.1.2. Disk structure

Each dust component is represented by a parameterized model of passive disk with a mass M_dust and inner and outer radii, r_in and r_out. Assuming it is similar to the gas density distribution in hydrostatic equilibrium, the dust density distribution is given by $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{z}}{\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right)}^{\mathrm{2}}\right]\mathit{,}& & \end{array}$(2)where z ≃ (π/ 2−θ)r is the vertical distance from the midplane, in the case of a geometrically thin disk (z ≫ r). The dust surface density Σ(r) and scale height H(r) are parameterized as Σ(r) = Σ_out⁽r/r_out^)p and H(r) = H_out⁽r/r_out^{)1 + β}, where β is the flaring exponent, Σ_out and H_out are the dust surface density and scale height at r_out. To enable a smooth decrease in density after r_out, we apply another p exponent to the surface density profile. We also include the possibility of rounding off or “puffing up” the inner rim of each component. For that, we artificially reduce or increase the dust scale height H(r) to reach a chosen value Ĥ_in at r_in. The width ϵr_in over which this reduction or increase is done – from no change at r₁ = (1 + ϵ)r_in to the dust scale height Ĥ_in – sets the “sharpness” of the rim. Following Ratzka et al. (2007), the modified dust scale height Ĥ(r), between r_in and r₁, writes as $\begin{array}{ccc}\mathit{Ĥ}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\left(\mathrm{1}\mathrm{-}\frac{\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{in}}}{{\mathit{r}}_{\mathrm{1}}\mathrm{-}{\mathit{r}}_{\mathrm{in}}}\right)\mathit{Ĥ}\mathrm{in}\mathrm{+}\left(\frac{\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{in}}}{{\mathit{r}}_{\mathrm{1}}\mathrm{-}{\mathit{r}}_{\mathrm{in}}}\right)\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}& & \end{array}$(3)

4.1.3. Dust properties

The dust grain properties are described by their optical constants taken from the Jena database². The mass absorption and scattering coefficients (in cm² g^-1), κ_abs(λ), and κ_sca(λ) are then computed using the Mie theory. The spherical shape approximation is usually safe for amorphous or featureless dust species for which the extinction properties are much less sensitive to grain shape effects than to crystalline material (e.g., van Boekel et al. 2005). Assuming a grain size distribution n(a) ∝ a^-3.5 (Mathis et al. 1977) with minimum a_min and maximum a_max grain sizes, the size-averaged mass absorption/scattering coefficient is obtained by adding the mass absorption/scattering coefficients of each grain size times their mass fraction. Then, the global κ_abs(λ) and κ_sca(λ) are derived by adding the size-averaged mass absorption/scattering coefficients of each dust species times their abundance to form a single “composite dust grain” that is the mix of the constituents. This averaged approach has already been justified in other disk studies (Mulders et al. 2011) and implicitly assumes a thermal coupling between the grains. This is expected since dust grains in disks are likely to be in the form of mixed aggregates in thermal contact. Each disk component has a homogeneous composition in the radial and vertical directions. No dust settling or radial segregation is considered here.

4.2. Modeling approach

We split our modeling approach into several steps, considering the mutual radiative influence between the different disk components and the associated dataset. For the star emission, we use the synthetic spectrum detailed in Sect. 2.5 and Table 1.

4.2.1. Inner disk

We address first the inner disk to reproduce the NIR SED and the PIONIER and AMBER V². Following Sect. 3.2, we set the inner disk’s inclination and position angle to i_disk = 20° and PA_disk = 112°. The inner rim is modeled as a vertical wall.

For the radial structure, r_in is first set to 0.2 au (Matter et al. 2014). We vary the surface density profile exponent p to explore different dust radial distributions. We also prevent Σ(r) from decreasing too abruptly after r_out by setting p = −10. Greater p values (e.g., ~−5) would induce too smooth a decrease, keeping r_out from representing the inner disk’s size.

For the vertical structure, the inner disk scale height at r_out is assigned a broad range of values, namely H_out/r_out = [ 0.03,0.05,0.1.0.2,0.3 ]. This encompasses the dust scale height values that are typically inferred or considered (0.1–0.15) for disks (e.g., Brown et al. 2007; Benisty et al. 2010).

We consider a dust mix of silicates and graphite. Including graphite is justified by its being a major constituent of the interstellar grains (Draine 2003) and by recent disks modeling that shows the importance of featureless and refractory species like graphite to explain the observations (Meeus et al. 2002; Carmona et al. 2014). We considered three graphite mass fractions (see Table 3) to evaluate its impact on the modeling. We assumed the same graphite composition as in Schegerer et al. (2008) and considered grains with a_min = 0.05μm and a_max = 0.2μm to maximize the IR opacity. Since graphite is a featureless species, this choice is not critical. For the silicate content, we assumed a pure iron-free olivine composition (see Juhász et al. 2010) with a_min = 0.1μm and 5 μm to keep or suppress the silicate feature. We fixed a_max = 20μm since the optical and NIR opacity contribution from mm grains is negligible.

The inner disk free parameters are p, r_out, and M_dust. The last is a lower limit since we only include μm-sized and sub-μm-sized grains. Then, the modeling steps are:

• computing a grid of 10 × 10 × 10 models on M_dust, p, r_out for each model setup (combination of a dust setup and a H_out/r_out value; see Table 4). We then calculate a reduced ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ separately for each dataset: the NIR SED (5 points from 1.25 μm to 4.5 μm), the dispersed PIONIER V², and dispersed AMBER V². The synthetic V² are computed from the RADMC3D images multiplied by a 2D Gaussian with a FWHM equal to the instrument FOV. The best-fit model of each model setup is found by minimizing the sum of the reduced ${\mathit{\chi }}_{\mathit{r}}^{\mathrm{2}}$ values;

• for each model setup, computing a grid of 8 × 8 × 8 models around the global minimum ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ and calculating the marginal probability distribution (exp [ −χ²/ 2 ]) to determine the best-fit value and 1σ uncertainty (68% confidence interval) on M_dust, p, and r_out. The best-fit model with the lowest reduced ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ value is chosen as the global best-fit solution;

• once the global best-fit solution is found, slightly varying r_in around its initial value (0.2 au) to further improve the fit to the NIR visibilities at high spatial frequencies.

4.2.2. Outer disk

We then address the outer disk to find a solution that is consistent with the Herschel/PACS data and the sub-mm SED. Based on the model of Matter et al. (2014), we set r_in = 5.6 au and r_out = 150 au and assume a pure-olivine composition for the silicate content, with a grain size distribution from 0.1 μm to 3 mm to account for the weak silicate emission feature in the N band and the low sub-mm spectral index (Sylvester et al. 1996). We set the surface density profile exponent to p = −1, as is typically observed in the outer regions of disks (e.g., Andrews & Williams 2007; Andrews et al. 2010). Then, two typical flaring profiles for irradiated disks are explored, i.e., β = 1 / 7 and β = 2 / 7 (e.g., Chiang & Goldreich 1997), and the dust scale height H_out/r_out is varied between 0.1 and 0.15, which is typical of disks (e.g., Benisty et al. 2010). We also vary the graphite abundance (from 0 to 20%) and the dust mass between 10^-6M_⊙ and 10^-4M_⊙.

4.2.3. Outer disk inner edge

We then focus on the inner edge to reproduce the MIDI data and improve the fit to the MIR SED. Both trace the geometry and flux fraction intercepted by the outer disk’s inner edge, hence its scale height and radial position. Since the SED measurements at 4.8 μm, 12 μm, and 25 μm are consistent with the Spitzer spectrum, modeling and fitting them will be sufficient and lead to a decent fit of the Spitzer spectrum. Reproducing the full Spitzer spectrum in detail is beyond the scope of this paper. Using the global best-fit solution of the inner disk, we compute a grid of 8 × 8 × 8 models on the outer disk inner radius r_in, the modified scale height Ĥ_in at r_in, and the rounding-off parameter ϵ (see Table 4). We then calculate a reduced ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ separately for each dataset, namely the dispersed MIDI visibilities (47 data points between 8 μm and 13 μm) and the MIR SED (at 4.8 μm, 12 μm, and 25 μm). The best-fit model is found by minimizing the sum of the reduced ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ values. The parameters’ best-fit values and 1σ uncertainties are computed in the same way as for the inner disk.

4.3. Results

4.3.1. A tenuous and extended inner disk

Table 5 shows the results of the inner disk modeling. Several models with different r_out values for the inner disk fit the NIR data equally well. Since our MIDI data partly resolve the inner disk (Matter et al. 2014), we can use them to break out the degeneracy in r_out and identify the global best-fit solution. For each inner disk solution, we computed a grid of 8 × 8 × 8 models on the outer disk’s inner edge parameters, and kept the model giving the lowest reduced ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ on the MIR SED and visibilities.

Size and location.

Only the inner disk models with r_out> 2 au are compatible with the MIR visibilities, especially from 8 to 9.5 μm. We constrained the inner disk’s outer radius to ${\mathit{r}}_{\mathrm{out}}\mathrm{=}\mathrm{2.5}{\mathrm{5}}_{-0.13}^{\mathrm{+}\mathrm{0.13}}$ au. With an uncertainty of (±27 pc) on the stellar distance, r_out varies from 2.0 to 2.9 au, which is still consistent with our conclusions. This small a variation in r_out would not significantly modify the conditions of irradiation of the inner disk and thus its temperature and brightness profile. From the best-fit solution, we varied r_in by steps of 0.05 au around r_in = 0.2 au to optimize the fit to the PIONIER and AMBER V². The best agreement is nevertheless found for r_in = 0.2 au. The r_in< 0.2 au values slightly increased the amount of unresolved inner disk emission and induced too high a K-band V² at low spatial frequencies. The r_in> 0.2 au values induced too prominent a lobe in the H-band V² profile at high spatial frequencies (see Appendix). As shown in Figs. 3 and 4, our model agrees well with the NIR SED and the PIONIER and AMBER V². Nevertheless, the best-fit H band V² are still slightly too high at high spatial frequencies, which probably suggests that the inner rim has a smoother shape than the assumed vertical wall. The modeled closure phases in the H band appears consistent with the measured ones. Figure 5 shows the best-fit model images at λ = 1.6μm and λ = 2μm.

Dust composition.

Our modeling favored the dust setup “2” (80% olivine +20% graphite, see Table 3). This implies a temperature for the hottest composite grains of T ≃ 1660 K at r_in = 0.2 au, which is on the warm side of sublimation temperatures for μm-sized silicates (Pollack et al. 1994). We also tested the effect of computing the temperature distribution separately for the two dust species. This led to olivine grains that were too cold (~1150 K) to be responsible for the NIR excess and to very hot graphite grains (~1900 K). The latter would induce too much flux at λ< 2μm and would disagree with our conclusions in Sect. 3.1. Therefore small graphite grains (or any refractory, featureless and efficient absorber/emitter) in thermal contact with the silicate grains seems required to induce enough heating and emission in the inner disk’s outer parts (>1 au) and to produce a spatially extended NIR-emitting region. The solution with a 100% graphite composition and ${\mathit{r}}_{\mathrm{out}}\mathrm{=}\mathrm{2.2}{\mathrm{0}}_{-0.10}^{\mathrm{+}\mathrm{0.10}}$ au is also compatible with the MIR visibilities. The dust compositions including smaller (sub-μm-sized) olivine grains (dust setups “1” and “3”) induce high inner rim temperatures (~1550 K), for which these small grains may sublimate, and overpredict the strength of the 10 μm feature in the Spitzer spectrum. This suggests that the inner disk does not contain a detectable amount of sub-μm-sized silicate grains and that the weak silicate feature originates in the outer disk.

Surface density profile.

We constrained the surface density profile to $\mathit{p}\mathrm{=}\mathrm{0.}{\mathrm{6}}_{-0.3}^{\mathrm{+}\mathrm{0.3}}$ (see Fig. 7). This positive profile suggests a radially increasing distribution of dust grains, and was required to create the extended NIR-emitting region predicted by the interferometric data. Such a radial distribution induces less shadowing from the inner rim and sufficient heating and re-emission from the outer parts of the inner disk. The radially integrated midplane optical depth (at λ = 0.55μm) reaches τ ≃ 1 only at r ≃ 0.6 au, as shown in Fig. 6. All the models with p< 0 induced too strong a flux contribution from the inner disk rim in the SED at λ ≲ 2μm and therefore a NIR-emitting region that is too spatially confined. This implied H-band and K-band V² that are too high at low spatial frequencies, and too low at high spatial frequencies since the STFR is underestimated (see Appendix).

Dust mass and surface density level.

We estimated a mass of ${\mathit{M}}_{\mathrm{dust}}\mathrm{=}\mathrm{8.}{\mathrm{8}}_{-3.7}^{\mathrm{+}\mathrm{4.2}}\mathrm{\times }{\mathrm{10}}^{-11}{\mathit{M}}_{\mathrm{\odot }}$ in small dust grains, which dominate the optical/IR extinction efficiency, and set the NIR emission level. The NIR absorption efficiency is dominated by the small graphite grains with a mass absorption coefficient 100 times larger than for the μm-sized olivine grains. Therefore the total dust mass estimate largely depends on the graphite fraction, and a significant amount of silicate grains could be hidden in the inner disk.

To estimate an upper limit on the dust mass, we extended the olivine grain size distribution up to 3000 μm (as in the outer disk) and kept the same absolute mass in graphite grains of our best-fit model. We then increased the mass in olivine grains (by steps of 10^-10M_⊙) until the modeled NIR emission (at 2.2 μm, 3.4 μm, and 4.8 μm) and/or NIR visibilities deviated by more than 3σ from the observed ones. As a result, we found an upper limit of about 6 × 10^-10M_⊙. Higher mass values implied a deviation of more than 3σ with respect to several AMBER V² measurements (all across the K band) and to the observed SED at 4.8 μm. We show in Fig. 7 the best-fit dust surface density and overplot its upper limit (dotted line).

Relative to the outer disk, the inner disk seems strongly depleted by at least ~10³. However, this dust depletion may be biased by a difference in dust composition between the inner disk that contains graphite grains and the outer disk that only contains olivine grains (see Table 3). However, the NIR mass absorption coefficient of the inner disk’s dust mixture is larger, by a factor 5, than that of the outer disk. In the MIR, the two coefficients are equal. The difference in dust composition thus cannot account for the dust surface density ratio of ~10³ between the two components.

Dust scale height.

A dust scale height of H_out ≃ 0.25 au (H_out/r_out = 0.1) at r_out = 2.55 au is favored. This translates to H_in ≃ 0.01 au (H_in/r_in = 0.044) at the inner rim. We recall that H(r) is the height from the midplane at which the dust density has decreased by a factor e^-0.5 (see Eq. (2)). Assuming a perfect dust-gas coupling, this would equal the pressure scale height of the gas in hydrostatic equilibrium. With a midplane dust temperature of 1660 K at r_in, the gas pressure scale height would be smaller than the dust scale height (H_in,gas = 0.006 au, i.e. H_in,gas/r_in ≃ 0.028). If dust grains are coupled to the gas, this suggests that 1) gas is actually hotter (with a required temperature of about 5000 K), as expected in the inner disk regions from thermo-chemical modeling (e.g., Woitke et al. 2009); and/or that 2) gas is vertically supported by additional sources, such as magnetic forces arising from the inner disk accretion driven by Magneto-rotational turbulence (Turner et al. 2014). With this dust scale height, the inner disk does not cast a significant shadow on the inner rim of the outer disk. Indeed, while the total integrated optical depth (at λ = 0.55μm) in the inner disk midplane is τ ≃ 4, it decreases quickly to τ< 1, above the midplane, so that the outer disk’s inner rim can intercept a large fraction of the stellar irradiation (see Fig. 6). The impact can be clearly seen in the emission increase at λ ≤ 8μm in the best-fit model SED (see Fig. 3) and in the ring-like morphology of the outer disk’s rim in the synthetic image at 10 μm (see Fig. 5).

4.3.2. An au-sized gap

As shown in Table 6, the outer disk starts at 5.7 au, and the gap, depleted in warm small grains, is au-sized (~3.2 au), which confirms the results of Matter et al. (2014). Considering the stellar distance uncertainty (±27 pc), the outer disk’s inner radius varies between 4.7 au and 6.7 au. The corresponding temperature variation is less than 5 K in the disk midplane and ~10 K in the higher disk surface. The gap width would range from 2.7 au to 3.8 au, which is still in the au-sized range. Therefore, we do not expect significant changes in the emission profile of the outer disk’s inner edge, hence no impact on the modeled MIR SED and visibilities.

The best-fit N-band visibilities are shown in Fig. 4. Both the level and the shape (sine-like modulation) in the measured visibilities are reproduced by our model. Figure 5 shows the best-fit model image at λ = 10μm (right panel), which shows the extended inner disk and the ring-like morphology of the outer disk’s inner rim. Then, we estimate an upper limit on the dust mass that could lie within the gap. Gaps opened by substellar companions can in principle filter dust grains partially decoupled from the gas (e.g., Rice et al. 2006). Therefore, the dust in the gap and the inner disk likely differs from the outer disk.

For the gap, we thus assumed the same dust composition (size distribution from 5 μm to 3000 μm for olivine grains and graphite fraction of 3%) and scale height profile as for the inner disk and considered the simplest case of a constant dust surface density profile (p = 0). We then varied the amount of dust in the gap until the modeled MIR SED (at 4.8 μm, 12 μm and 25 μm) and/or the MIR visibilities start deviating by more than 3σ relative to the observations. As a result, we estimated an upper limit of ~2.5 × 10^-10M_⊙. Above this value, the modeled SED overestimated the observed SED at 4.8 μm by more than 3σ, while two of the modeled MIR visibilities became lower, by more than 1σ, than the MIDI visibilities between 11.5 and 13 μm. As the gap is being filled in, the amplitude of the sine-like modulation also decreases in the modeled MIR visibilities. Interestingly, the mass estimate inside the gap is comparable to the upper limit found for the inner disk. We recall that other mass reservoirs could lie within the gap, such as cold mm-sized grains, pebbles, or minor bodies, which do not significantly contribute to the disk IR emission.

4.3.3. An outer disk with a “smooth” inner edge

A 1.3 × 10^-4M_⊙ outer dust disk with a flaring index β = 1 / 7 and a scale height H_out/r_out = 0.135 at r_out consistently reproduces the Herschel/PACS and sub-mm photometry measurements. A 100% olivine dust composition with 0.1 μm <a< 3 mm was sufficient to account for the Spitzer/IRS measurements and its weak 10 μm silicate feature, along with the sub-mm spectral slope. Adding graphite in our model was not required. The best-fit broadband SED is presented in Fig. 3.

Our MIDI data favored a slightly rounded shape for the outer disk’s inner rim over a purely vertical wall. A reduction of the scale height of the dusty disk from 1.2 × r_in = 7.4 au to r_in = 5.7 au (ϵ = 0.3), leading to H_in/r_in = 0.04 (instead of H_in/r_in = 0.085), was needed to reproduce the modulation in the measured MIDI visibilities. A purely vertical wall induced too strong a modulation with the apparition of a pronounced lobe for the third visibility measurement. The rim shape could be explored further in a future study using hydrodynamical simulations to infer the mass of an hypothetical substellar companion, as previously done for HD 100546 (Mulders et al. 2013).

4.4. Caveats and limitations

First, we only explored a limited set of dust compositions. Considering crystalline silicates (Juhász et al. 2010) or other featureless species (e.g., iron) may influence the dust mass and change the disk temperature and brightness profile. However, the available data do not allow us to disentangle more complex dust compositions.

Our best-fit solution assumed a radially and vertically homogeneous composition. Considering dust radial segregation or settling will affect the disk opacity profile and may lead to different scale heights and surface density profiles that are compatible with the data. However, the available dataset on HD 139614 currently lacks resolved mm observations that probe larger grains and cannot allow us to investigate dust segregation.

With a disk orientation close to face-on, the assumption of isotropic scattering may have induced too much NIR flux being scattered in our line of sight and biased our estimation of the inner disk dust mass. However, this should be limited since the NIR scattering opacity is strongly dominated by the sub-μm-sized grains of our disk model. These grains are still in the Rayleigh regime, in the IR, and scatter almost isotropically.

The 2D axisymmetric geometry we assumed is relevant for identifying the main disk structural aspects but not possible asymmetries. However, this is not problematic here since the available closure phases do not suggest strong brightness asymmetries. We also considered a parameterized dust density that follows the prescription of a gaseous disk in hydrostatic equilibrium. A self-consistent computation of the disk vertical structure, as in Menu et al. (2014), may have modified our derived dust scale height. However, this approach still relies on a thermal coupling between dust and gas and does not necessarly allow exploring vertical structures departing from the “hydrostatic” prescription. Finally, the power-law surface density we used is an approximation that allows trends in the radial dust distribution to be identified. It remains very simplified compared to the dust density profiles derived from hydrodynamical simulations coupling gas and dust (e.g., Fouchet et al. 2010) even though the latter can approach power-law profiles in some cases (e.g., Birnstiel et al. 2012).

5. Discussion

5.1. A group I Herbig star with a pre-transitional disk

Our new results on the gapped disk of HD 139614 support the idea that most of the group-I Herbig objects are in the disk-clearing stage, as suggested by (Maaskant et al. 2013) and confirmed recently by Menu et al. (2015) from a global analysis of all the MIDI data obtained on Herbig stars. Interestingly, Menu et al. (2015) find that Group I objects known to harbor very large gaps of tens up to a hundred au, such as HD 142527 (e.g., Fukagawa et al. 2006), would also have a small gap in the inner region (≤10 au). Moreover, several supposedly gapless and settled Group II disks show hints of a very narrow gap in their inner (~au) dust distribution. HD 139614 stands out here since it is a rare, if not unique, case of Herbig star’s disk for which an au-sized gap (~3 au) has been clearly spatially resolved in the dust distribution (see Espaillat et al. 2014, for a review on gap sizes). MIR gaps correspond to local depletion in small sub-μm to μm-sized grains, which are good tracers of the gas (e.g., Barrière-Fouchet et al. 2005). As a result, HD 139614 is a relevant laboratory for investigating gap formation and disk dispersal in the inner regions, which are only reachable by IR interferometry.

As shown in Sect. 4, self-shadowing by the inner disk does not convincingly reproduce the pretransitional features of the HD 139614’s SED and the interferometric data. Among the main disk clearing mechanisms, two could explain a narrow gapped structure: photoevaporation and planet-disk interaction (Espaillat et al. 2014). The photoevaporation scenario may be relevant for HD 139614 given its accretion rate (≲10^-8M_⊙/yr, Garcia Lopez et al. 2006), which is in the range of predicted photoevaporative mass-loss rates 10^-10–10^-8M_⊙/yr (Alexander et al. 2014). Moreover, the photoevaporative wind is expected to first open a gap at the critical radius where the mass-loss rate is at its maximum (Alexander et al. 2014). In the extreme-UV regime, this critical radius was estimated to be R_c,EUV ≃ 1.8 M_∗/M_⊙ au. For HD 139614 (M_∗ = 1.7 M_⊙), R_c,EUV ≃ 3.1 au, which is consistent with the gap location. However, several aspects of the HD 139614 disk hardly appears consistent with this scenario:

• the detection of warm molecular gas in the 0.3–15 au region from observations of therovibrational emission of ¹²CO and its isotopologue lines with CRIRES (Carmona et al., in prep.). Preliminary results suggest at least 0.4 M_⊕ and up to 2 M_Jup of gas within 6.5 au. Moreover, the CO lines present a clear double-peaked profile that is indicative of a gaseous disk in Keplerian rotation, which is not consistent with the presence of a photoevaporative molecular disk wind;

• the non-detection of the [OI] line at 6300 Å, which suggests the absence of a disk wind (Acke et al. 2005);

• the short timescale (<10⁵ yrs) expected for the inner disk dissipation (in gas and dust) after the photoevaporative wind has opened a gap (Alexander et al. 2014).

For these reasons, we then focus on the planet-induced gap scenario. Single massive planets are expected to carve a few au wide gap in the gas; the precise width and depth depending on the planet mass and disk viscosity (e.g., Baruteau et al. 2014). A similar clearing is expected for the gas-coupled small dust grains. Moreover, the tenuous HD 139614 inner disk suggests a drastic evolution and differentiation of the inner regions, which is possibly reminiscent of dust filtration by a planet on the gap’s outer edge (e.g., Rice et al. 2006). Interestingly, Carmona et al. (2014) find similar differences for the HD 135344 B transition disk. To determine whether our observational constraints on the dust can sustain the planet-induced gap scenario, we performed a comparative study with hydrodynamical simulations of gap opening by a giant planet, adapted to the case of HD 139614.

5.2. Investigation of disk-planet interaction

5.2.1. Hydrodynamical setup

We use the code FARGO-2D1D³, which is designed to study the global evolution of a gaseous disk perturbed by a planet (Crida et al. 2007).

Our planet is on a circular orbit of radius a_pl and does not migrate. We consider planet masses of 1, 2, and 4 × 10^-3M_∗, with M_∗ the mass of the star. For HD 139614, this corresponds to 1.7, 3.4, and 6.8 M_Jup. The 2D grid extends from r = 0.35 a_pl to 2.5 a_pl with N_r = 215 rings of N_s = 628 cells arithmetically spaced so that the resolution is dφ = dr/r = 0.01 at the planet location, where φ is the azimuth. The 1D grid extends from 0.04 a_pl to 20 a_pl and has open boundary conditions to allow for disk spreading. The simulation uses a locally isothermal equation of state: P = c_S²Σ with P the pressure, Σ the gas surface density, and c_s = HΩ the sound speed. Here, c_s is fixed such that H/r is constant, with H the gas pressure scale height and Ω the Keplerian angular velocity. We set H/r = 0.04 to be consistent with the H/r value inferred for the dust at the outer disk’s inner edge. The disk scale height influences the gap’s depth but is not a critical parameter here, given our limited constraint on the dust gap depth. The initial density profile is Σ = Σ₀(r/a_pl)^-1. The gas viscosity is $\mathit{\nu }\mathrm{=}{\mathrm{10}}^{-5}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}{\mathit{a}}_{\mathrm{pl}}}$, with G the gravitational constant; at the planet location, the Reynolds number is R = r²Ω /ν = 10⁵, which gives α = 6.25 × 10^-3 in the Shakura & Sunyaev (1973) prescription, but here the viscosity is assumed to be independent of time and space. With this equation of state and the planet on a fixed orbit, the units are arbitrary: regardless of the physical value of a_pl and Σ₀, the results scale accordingly. For easier comparison with the case of HD 139614, we have set a_pl = 4.5 au, and Σ = 15 g cm^-2 at 10 au as extrapolated from the dust surface density at 10 au (Σ_dust = 0.15 g cm^-2), assuming a gas-to-dust ratio of 100. Figure 7 shows the gas density profiles.

5.2.2. Simulated gas density profile

With the considered disk viscosity and scale height, a ~3 M_Jup planet at 4.5 au produces a gas gap between ~3 and 6 au. This is consistent with the dust gap width, the inner disk’s outer radius (r_out = 2.55 ± 0.13 au), and the outer disk’s inner radius (r_in = 5.7 au) of our best-fit dust model. The gas gap appears shallower than in the dust, which is expected since gas pressure minima are strongly depleted in dust (e.g., Rice et al. 2006). The 1.7 and 6.8 M_Jup planets produced a gap that was either too narrow (≲2 au) or too wide (≥4 au).

Another noticeable aspect is the surface density profile. From the radiative transfer, we highlighted a radially increasing dust surface density profile in the inner disk, while the usual profile in r^-1 was kept for the outer disk. The gas surface density profile appears very consistent with this. In the inner disk, the gas surface density increases from 0 at the inner edge until it reconnects to the original profile in r^-1 at r> 5 au. Moreover, if fitted by a power law, the gas profile’s exponent in the inner disk (~0.6) is similar to the dust profile exponent (0.64). This radially increasing gas surface density, which we highlighted for the small dust grains, too, is intrinsic to the disk viscous evolution. As shown in Lynden-Bell & Pringle (1974) and Crida & Morbidelli (2007), such a profile naturally appears in the inner region of a viscously spreading gaseous disk, which has a finite inner edge that is not too close to the star (possibly 4% of the planet’s orbital radius for HD 139614). Interestingly, a radially increasing gas density profile was required to describe the CO emission within the cavity of the HD 135344B disk (Carmona et al. 2014).

After a gap has been opened, the gas surface density ratio between the inner and outer components is mainly ruled by the amount of gas that can cross the gap. Moreover, if the gap is opened close to the disk’s innermost edge, where the density profile is radially increasing, the inner disk’s surface density will be lower than in the outer disk (see Fig. 6 of Crida & Morbidelli 2007). Figure 7 shows an inner disk surface density deficit of a factor 10 in gas, relative to the outer disk. This does not vary much as the whole disk viscously spreads and the gas density decreases as a whole. Even after 1 Myr, the inner disk is much less depleted in gas (~10) than in dust (~10³). This suggests that a planet opening a gap consistent with our observations cannot affect the accretion flow enough to the inner regions. With the 6.8 M_Jup planet or more viscous disks, we could not deplete the inner disk any more. Having the planet on a fixed orbit actually helps the inner disk’s depletion since a planet migrating in type-II migration follows the disk viscous spreading without hampering it. Also, when the disk is less massive than the planet, no migration occurs, and the situation is similar to our simulation (see Eq. (11) of Crida & Morbidelli 2007). However, during the disk evolution, dust grains grow and decouple from the gas, or fragment and stay coupled. The inner disk depletion may thus differ for the small dust grains and the gas, as shown hereafter.

5.2.3. Grain growth and fragmentation in the inner disk

The outer edge of a gap opened by a planet is a pressure maximum and can trap the partially gas-decoupled dust grains. In the inner disk, these grains drift inward and are eventually lost to the star. The radial drift velocity is highest when the Stokes number St = Ωρ_da/ρ_gc_s ~ 1, where Ω is the Keplerian angular velocity, ρ_d the intrinsic dust density, a the grain size, ρ_g the gas density, and c_s its sound speed. The optimal size for radial migration is ${\mathit{a}}_{\mathrm{opt}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{/}\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\rho }}_{\mathrm{d}}$ in the midplane (Fouchet et al. 2010). Since the gap acts as a dust filter, the dust-to-gas ratio in the outer disk is likely to be close to the interstellar value of 0.01, leading to Σ_g ~ 10 g cm^-3 at the gap’s outer edge (see Fig. 7).

Our simulations indicate a gas surface density that is ten times lower in the inner disk, which implies a_opt ≃ 1 mm for typical values of ρ_d (1–3 g cm^-3). Grain growth is expected to be efficient in the inner disk regions (Brauer et al. 2008; Laibe et al. 2008) and to bring small grains to mm or cm sizes. Once grains have grown up to sizes ~a_opt in the inner disk, they migrate inwards quickly and sublimate when approaching the star. This could occur on a timescale of 10² yr for cm-sized bodies at 1 au in a typical accreting disk (Brauer et al. 2008). Only the small grain (μm-sized or smaller) reservoir, partly replenished by fragmentation, would stay coupled to the gas and remain in the inner disk longer. This depletion of the small grains reservoir will lead to a tenuous inner disk with reduced optical depth. This is supported by hydrodynamical simulations of a disk of gas and dust containing a planet, where dust grain growth, fragmentation, and dynamics are treated self-consistently (Gonzalez et al. 2015). Although the disk and planet parameters are different from ours, they show an inner disk that depletes more efficiently in dust than in gas. For most of the considered fragmentation velocities (<20 m s^-1), growth is fast for small grains but more difficult as they keep growing. Fragmentation either prevents grains from growing above a_opt, promoting their fast migration, or breaks them down again to smaller sizes, thereby preventing their fast migration and partly replenishing the small grain reservoir. The inner dust disk thus drains out on timescales shorter than the viscous disk evolution and only keeps the remaining small grains and the larger bodies that could have overcome the radial-drift and fragmentation barriers. After 10⁵ yr, the inner-to-outer disk dust density ratio has dropped to ~10^-3 and leads to a tenuous inner disk with a reduced optical depth, which agrees with our observational results. Quantitatively reproducing our observed dust disk with similar simulations is beyond the scope of the paper. Nevertheless, the qualitative agreement supports the scenario of a planet-induced gap in the HD 139614 dust disk.

6. Conclusion

This first multiwavelength modeling of the dust disk around HD 139614 provided the following results:

• We confirmed a gap structure, between 2.5 and 5.7 au, depleted (but not necessarily empty) in warm sub-μm- and μm-sized grains. HD 139614 appears to be a rare case of a Herbig star with a narrow au-sized gap in its dust disk.

• The NIR-emitting region was found to be spatially extended from 0.2 to 2.5 au, with a radially increasing surface density profile (p> 0), a dust scale height of H ~ 0.01 au at r = 0.2 au, and a strong dust depletion of at least ~10³ (in surface density) relative to the outer disk. This suggests a drastic evolution of the inner regions.

• With a dispersion of ≲2–5° around zero, the PIONIER closure phases does not suggest significant brigthness asymmetries. Notably, HD 139614 is unlikely to be a tight binary system with a companion more massive than about 0.11 M_⊙.

Currently available data and modeling do not support self-shadowing from the inner disk or photoevaporation as the origin of the dust gap in the HD 139614 disk. We thus tested whether the mechanism of disk/planet interaction could be viable against our constraints obtained on the dust from the radiative transfer. Assuming that small dust grains, probed by IR interferometry, are coupled to the gas, we performed hydrodynamical simulations of gap formation by a planet in a gaseous disk using the code FARGO-2D1D. It appears that:

• For typical values of disk viscosity and scale height, a 3 M_jup planet, fixed at ~4.5 au, could produce a gap in the gas consistent with the dust gap, although a bit shallower.

• The radially increasing dust surface density in the regions interior to the gap is reproduced in the gas. The original gas density profile (∝r^-1) is maintained in the outer disk for a long time. The radially increasing profile is predicted analytically and represents the accretion profile of the gas (and coupled dust grains) close to the inner disk rim. To our knowledge, it is the first time that this behavior is observationally identified in the innermost dust distribution of a disk.

• The inner disk dust depletion (~10^-3) is not reproduced in gas (~0.1). However, this can be explained by the radial drift and growth and fragmentation processes affecting dust and is predicted by hydro simulations coupling gas and dust.

All of this supports the hypothesis of a planet-induced gap in the HD 139614’s dust disk and reinforces the idea that disks around Group-I Herbig stars are already in the disk-clearing transient stage. Confirming the planet-induced gap scenario will require further observations such as upcoming MIR imaging of the gap with the second-generation VLTI instrument MATISSE (Lopez et al. 2014). Moreover, ALMA will be able to probe the mm grains and obtain a robust estimate of their mass at r< 10 au, while mapping the gas density distribution and dynamics inside and outside the gap. Although the HD 139614 gap seems out of reach for VLT/SPHERE, direct imaging of the outer regions could be attempted to detect other signatures of disk-planet interaction. HD 139614 will also constitute an exciting target for the future E-ELT/MICADO, PCS, and METIS instruments.

---

Acknowledgments

We thank the anonymous referee for the comments that helped to improve this manuscript significantly. A. Matter acknowledges financial support from the Centre National d’Études Spatiales (CNES). J.-F. Gonzalez is grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the program “Investissements d’Avenir” (ANR-11-IDEX-0007) of the French government operated by the ANR. This work is supported by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR-10-BLAN-0511).

References

Appendix A: Parameter effects in radiative transfer modeling

In this section, we illustrate the individual effect of the main disk parameters on our radiative transfer modeling. From our best-fit model shown in Table 6 and Figs. 3 and 4, we vary each parameter separately around its best-fit value, and show the impact on the modeled observables (SED, NIR V², and MIR visibilities). In this way, we illustrate the effects of the main disk parameters and how much they can be constrained by our available dataset.

Appendix A.1: Inner disk’s inner radius

We show here the effect of varying the inner radius of the inner disk on the SED and the H-band and K-band V². As mentioned in Sect. 4.3.1, r_in< 0.2 au values slightly increase the amount of unresolved inner disk emission. As shown in Fig. A.1, this induces higher modeled H-band V² combined with a more gradual V² decrease as a function of spatial frequency. However, the modeled K-band V² become overestimated at the lowest spatial frequencies (smallest baselines). The r_in> 0.2 au values induce too prominent a lobe in the PIONIER V² profile at high spatial frequencies. The modulation amplitude thus increases too much as we increase r_in above 0.2 au.

Appendix A.2: Surface density p exponent

We examine here the effect of varying the power-law exponent p of the inner disk’s dust surface density on the NIR SED and V². As mentioned in Sect. 4.3.1 and shown in Fig. A.1, the models with p ≤ 0.3 start inducing too strong a flux contribution from the inner disk rim and an overestimated NIR emission at λ ≲ 2μ in the modeled SED. This too spatially confined NIR-emitting region implies H-band and K-band V² that is too high at low spatial frequencies and V² that is too low at high spatial frequencies since the stellar-to-total flux ratio is lower than in the case of our best-fit model, especially at λ ≲ 2μ.

Appendix A.3: Inner disk’s outer radius

We examine here the effect of varying the outer radius r_out of the inner disk. In addition to the NIR SED and V², we also consider that the N-band visibilities illustrate the constraint provided by the MIDI data on the extension of the inner disk. Figure A.2 shows that decreasing or increasing the size of the inner disk by at least 0.5 au, relative to our best-fit model, has a noticeable effect on all the observables. A smaller inner disk will induce too high N-band visibilities especially between 8 and 9.5 μm, while larger and over-resolved inner disks (r_out ≥ 3 au) will produce N-band visibilities that are too low and flatter. The effect on the H-band and K-band V² is directly correlated to the change in the NIR emission level in the modeled SED. Indeed, a smaller inner disk (r_out ≲ 2.0 au) is more optically thick (for the same dust mass) and will produce more NIR emission. This implies an overestimation of the disk contribution and consequently a lower stellar-to-total flux ratio, which will then induce lower NIR V² (see Fig. A.2), especially at long baselines where the inner disk is fully resolved. For a larger inner disk (r_out ≥ 3 au), the effect on the observables is the opposite.

Appendix A.4: Inner disk scale height

We show in Fig. A.3 the effect of varying the dust scale height H_out of the inner disk at its outer radius r_out. A slight departure from the best-fit value H_out/r_out = 0.1 has a noticeable impact on all the NIR observables. Indeed, increasing the scale height increases the amount of stellar flux captured and reprocessed by the inner disk and thus the NIR emission level. This implies a decrease in the stellar-to-total flux ratio and therefore lower NIR V², especially at high spatial frequencies. We thus quickly underestimate the V² level in H and K bands. On the other hand, slightly reducing the dust scale height to H_out/r_out ≃ 0.07 produces modeled SED and NIR V² that are still consistent with the observations.

Appendix A.5: Outer disk’s inner radius

Finally, we show in Fig. A.3 the effect of varying the inner radius of the outer disk r_in on the MIR visibilities, which are directly probing this region of the HD 139614 disk. Indeed, it appears that a deviation by more than ~0.3 au, which is the 1σ uncertainty on r_in (see Table 6), induces a noticeable shift in the sine-like modulation in the modeled MIR visibilities. The latter are no longer consistent with the measured MIDI visibilities, although the visibility level remains similar. No significant effect can be seen in the MIR SED.

All Tables

All Figures