The GRAVITY Young Stellar Object survey -- I. Probing the disks of Herbig Ae/Be stars in terrestrial orbits

The formation and the evolution of protoplanetary disks are important stages in the lifetime of stars. The processes of disk evolution and planet formation are intrinsically linked. We spatially resolve with GRAVITY/VLTI in the K-band the sub au-scale region of 27 stars to gain statistical understanding of their properties. We look for correlations with stellar parameters, such as luminosity, mass, temperature and age. Our sample also cover a range of various properties in terms of reprocessed flux, flared or flat morphology, and gaps. We developed semi-physical geometrical models to fit our interferometric data. Our best models correspond to smooth and wide rings, implying that wedge-shaped rims at the dust sublimation edge are favored, as found in the H-band. The closure phases are generally non-null with a median value of ~10 deg, indicating spatial asymmetries of the intensity distributions. Multi-size grain populations could explain the closure phase ranges below 20-25 deg but other scenarios should be invoked to explain the largest ones. Our measurements extend the Radius-Luminosity relation to ~1e4 Lsun and confirm the significant spread around the mean relation observed in the H-band. Gapped sources exhibit a large N-to-K band size ratio and large values of this ratio are only observed for the members of our sample that would be older than 1 Ma, less massive, and with lower luminosity. In the 2 Ms mass range, we observe a correlation in the increase of the relative age with the transition from group II to group I, and an increase of the N-to-K size ratio. However, the size of the current sample does not yet permit us to invoke a clear universal evolution mechanism across the HAeBe mass range. The measured locations of the K-band emission suggest that these disks might be structured by forming young planets, rather than by depletion due to EUV, FUV, and X-ray photo-evaporation.


Introduction
An understanding how disks of gas and dust around young stars evolve and dissipate is essential to gaining a better understanding of planet formation theories.The morphology of protoplanetary disks provides information on their evolution and may relay evidence relevant to planet formation.In recent years, outer disk features have been probed in detail through scattered light imaging with instruments like SPHERE in the optical and near-infrared (Beuzit et al. 2019) or through imaging at (sub-Fig.1. GRAVITY (squares, diamonds, triangles) and Gaia (circles) observations (Vioque et al. 2018) of Herbig Ae/Be stars as put in the Hertzsprung-Russell diagram.The symbols code the Meeus groups (Meeus et al. 2001).The dashed line corresponds to the 2.5 Ma isochrone.The numbering refers to Table 1. and broad diversity of the morphology of the outer disks (Garufi et al. 2018).
Currently, the level of detailed understanding with regard to the disk structure that is accessible with SPHERE or ALMA for the disks' outer region (i.e., ∼ 20 to ∼ 500 au) is lacking when we are looking at the innermost regions of a few au in size.However, the dynamical processes at work in these first central au might be just as complex and as diverse as what has been found in the outer regions.They most likely drive disk evolution in a key region where terrestrial planets form or migrate.The inner rim is the place where dust is processed thermally and whence it can be redistributed to the outer disk.There is evidence that the magneto-rotational instability activation in the innermost disk and its suppression in the dead zone are at the basis of the inside-out planet formation scenario (Mohanty et al. 2018).In this context, dead zone and dust traps, resulting in the possible formation of local asymmetries, play a crucial role in the vortex generation and planetesimal formation (Flock et al. 2016(Flock et al. , 2017)).The interplay between these different aspects of the physics of disks is one of the important open question in the field, and a direct comparison between inner/outer disk phenomena is a possible way forward with regard to sharpening our view on young disk evolution.
The spectral energy distribution (SED) of the disk is too degenerate to allow the inner rim to be fully constrained, whereas zooming into the innermost regions of disks in the nearest starforming regions (d∼ 120-130 pc) requires observations at milliarcsecond (mas) resolution to reach the desired (sub-)au scales.As of a few recent decades, long-baseline infrared interferometry (LBI) has demonstrated such a unique ability, as, for example, with the IOTA interferometer (Millan-Gabet et al. 2001) and the four-telescope H-band instrument PIONIER (Le Bouquin et al. 2011) at the VLTI.Despite its poorer coverage of the spatial frequency plane in comparison to ALMA, LBI was able to bring new insights in the fundamental morphological properties of the inner disks of Young Stellar Objects (YSO; see Dullemond & Monnier (2010) and Kraus (2015) for a review).The size-luminosity relation adequately correlates the measured inner radii of disks with the luminosities of the central stars (Monnier et al. 2005), which is in strong agreement with the presence of a directly illuminated rim at the dust sublimation radius (Natta et al. 2001).To explain both photometric and interferometric measurements, and, in particular, the near-infrared (NIR) contribution of the rim to the Herbig NIR excess, several disk models have been proposed as vertically extended geometry (Natta et al. 2001;Flock et al. 2017), disk wind components (Vinković & Jurkić 2007;Bans & Königl 2012), or a tenuous dusty halo around the disk inner regions (Vinković et al. 2006;Maaskant et al. 2013).Regarding the inner rim structure, a wide range of shapes have been proposed: a puffed-up rim due to direct heating from the star, with a curved bright side of the rim linked to the gradient of evaporation temperature (Natta et al. 2001;Isella & Natta 2005); a wedge-shaped rim attributed to differences in cooling properties between grains of different sizes (Tannirkulam et al. 2007); a vertical structure of the rim when considering different grain sizes and composition was self-consistently calculated by Kama et al. (2009); see also Klarmann et al. subm.).The structure and the composition of these regions are still a matter of debate.
With a particular focus on Herbig Ae/Be intermediate-mass objects, the H-band PIONIER Large Program (Lazareff et al. (2017); L17) provides a detailed view on the inner rim region for a large number of stars (∼ 50) and it has shown that the inner rim, where dust is sublimated, appears smooth and wide, which is in agreement with pioneering CHARA observations (Tannirkulam et al. 2008) and disk models with multi-grain populations.Because the disk/star contrast increases in the K-band, and the disk morphology in the continuum can be better compared to  Vioque et al. (2018): distance d, effective temperature T eff , luminosity L, mass M, age, extinction A v , and binarity flag Bin.The last column gives the group of the Meeus classification (Meeus et al. 2001;Juhász et al. 2010) when available.Spectral type is from CDS-Simbad.Symbol indicates the low-quality Herbig Ae/Be sample in Vioque et al. (2018) et al. 2017) and which includes about a hundred young Herbig Ae/Be and T Tauri stars observed in the K-band continuum and across spectral lines such as the Hydrogen Brγ and the CO band heads.Our targets span a wide range of masses (from 1 to 10 solar masses), ages (from 10 4 to 10 7 a), accretion rates (from 10 −6 to 10 −9 M /yr −1 ), and disk morphologies (i.e., full and transitional disks), as we aim to explore the dependence of our results on the properties of the central star (spectral type, mass, luminosity, age) and those of the disk (reprocessed flux, presence of gaps, flared/flat morphology).
This paper reports on new K-band continuum observations with GRAVITY for 27 Herbig Ae/Be stars observed during the first two years of the GTO YSO program, with a focus on possible trends among different properties identified within the systems.The paper is organized as follows: GRAVITY observations are described in Sect.2; we present our dataset in Sect.3; our data analysis is detailed in Sect.4; the results are given in Sect. 5 and discussed in Sect.6.

Sample
We observed a sample of 27 Herbig Ae/Be stars covering a broad range in luminosity (1-10 4 L ), effective temperature (6375-17500 K), mass (1.4-12.2M ), and age (0.04-14.5 Ma), as detailed in Table 1 and Fig  ONIER Large Program since they correspond to the brightest candidates, which will allow us to compare the disk morphology between the H-and K-bands.Since flat, flared, gapped, and ungapped disks are expected to have different signatures when observed at high angular resolution, we selected our targets to exhibit IR-excess and to include 10 group I objects, 12 group II objects, and 5 unclassified objects following the Meeus et al. (2001) classification (see Sect. 6.3 for details).

Observations
All our targets were observed with the GRAVITY instrument, using the four 1.8-m Auxiliary Telescopes, with the exception of HD 100546, which was observed using the four 8-m Unit Telescopes.The interferometric signals were recorded on six baselines simultaneously on the fringe tracker (FT) and the science channel (SC) detectors.The SC records the interferometric observables at high spectral resolution (R ∼ 4000) over the whole K-band with individual integration times of, typically, 10 s to 30 s.The FT records at low spectral resolution (five spectral channels over the K-band) at frame rates ranging from ≈300 to ≈ 900 Hz (Lacour et al. 2019).Working at such a high speed allows for the atmospheric effects to be frozen.Each observation file corresponds to five minutes on the object.To calibrate the atmospheric transfer function, we interleaved our target observations with observations of interferometric calibrators that are established as single stars, bright, small enough, and close to the target.A detailed log of the observations is presented in Table A.1.

Data reduction
All the data were reduced using the GRAVITY data reduction pipeline (Lapeyrere et al. 2014).In this paper, we focus on the FT observations for probing the inner dust rim in the K-band continuum as the turbulence effects are much lower at the speed of the FT.The SC observations at a spectral resolution of about 4000 will be used to constrain the Brγ line emitting regions (Garcia-Lopez et al., in prep).For each file, we obtained six squared visibilities and four closure phases for six spectral channels.In the following, we discard the first spectral channel that might be af-fected by the metrology laser working at 1.908 µm and by the strong absorption lines of the atmospheric transmission.

Dataset
All the calibrated data are displayed in Appendix B. Our GRAV-ITY observations span a spatial frequency range between 5 Mλ and 60 Mλ.Our maximum angular resolution of λ/2B is of about 1.7 millisecond of arc (mas) for a longest baseline B of 130 m, which corresponds to about 0.24 au at a distance of 140 pc.Except for 5 objects that are partially resolved (HD 85567, HD 114981, HD 158643, HD 259431, and PDS 27) with minimum squared visibilities larger than 0.4, our squared visibility measurements go below 0.4.A few objects are fully resolved and their visibilities are close to 0 at the longest baselines (HD 37806, HD 45677, HD 98922, HD 100546, HD 190073, and R CrA).For each object, the visibility variation with the spatial frequency allows the extent of the environment to be determined through model fitting (Sect.4).
The closure phases are related to the asymmetry of the environment (Haniff 2007): the closure phase is null for a centrosymmetric target; in the case of asymmetrical objects, the closure phase signal increases with spatial resolution up to the point where the spatial scale of the asymmetry is resolved; beyond this point, the closure phase signal decreases.In looking at our sample, the histogram of closure phase variations over the spatial frequency range, i.e. the peak-to-peak difference between the minimum and the maximum closure phases, shows that only one star among 27 has a variation smaller than 5 • .About half of our targets (12 among 27) have closure phase variations smaller than 10 • ; about a third (10 among 27) has a closure phase variation ranging between 11 • and 25 • ; and the five other targets have strong closure variations that exceed 30 • (Fig. 2).As expected, the strong closure phase signals are observed for well-resolved targets (HD 45677, HD 98922, HD 144432, HD 144668, HD 179218, and R CrA).Conversely, several targets are well-resolved with visibilities squared below 0.2 and do not exhibit closure phase variations higher than 15 •  (HD 37806, HD 58647, HD 97048, HD 135344B, HD 142527,  HD 150193, HD 163296, SAO 206462).

Geometric modelling
We use the same analysis tools as in the PIONIER survey and developed similar geometric models as described in L17.In the following section, we only recall the main steps involved in the fitting method.

Model visibility
At near-infrared wavelengths, the disk emission essentially comes from a compact region near the sublimation radius, which justifies the adoption of simple geometric models capable of capturing the general morphological properties of the environment.This would be less true at mid-infrared wavelengths, where the intensity distribution is spatially more extended.Accordingly, our model includes a point-like central star and a circumstellar environment.The complex visibility at the spatial frequencies (u, v) and at the wavelength λ is, therefore, described by a linear combination of the two components as follows: where the star visibility equals 1 since it is assumed to be unresolved (even for the closest star forming regions, the interferometric angular resolution is not high enough to resolve the stellar photosphere).V c stands for the visibility of the circumstellar environment, f s and f c for the fractional flux contributions of the star and the circumstellar environment, respectively ( f s + f c = 1).With regard to targets for which the scattered light amount is significant (mainly for transitional disks; see L17), an additional flux contribution (called halo and characterized by a flux ratio f h ) can be considered.The halo is assumed to be much more extended than the angular resolution of the interferometer.Its visibility is, thus, null.Because the halo contributes to the fractional flux (i.e., f s + f c + f h = 1), its net effect is that for very short baselines, the visibility is smaller than 1.At very long baselines, the environment might be completely resolved and its visibility V c reaches zero.In that case, the visibility curves display a plateau, which corresponds to the stellar flux contribution f s .

Flux ratio determination
To determine the flux ratio between the star and the environment for the Herbig Ae/Be stars, we use the same photometry datasets as the PIONIER LP.For each target, we fit the SED under the same assumptions as L17: we perform a least-square fit of the observed fluxes in the U, B, V, R, I, J, H, K Johnson-Cousins bands; the stellar photosphere is represented by the relative fluxes in each Johnson-Cousins band for the relevant spectral type, while the dust emission is modeled as a single temperature blackbody.The fitted parameters are the absorption-free stellar flux in V band, the absorption-free dust flux in the K-band, the extinction Av, and the blackbody temperature of the dust component.The K-band environment contribution is expected to be higher than in the H-band one.Being able to accurately determine the stellar flux with only our GRAVITY observations would require to reach the plateau when the circumstellar environment is fully resolved, which is only the case for a few stars.
The SED fit provides us with the fractional flux f c of the circumstellar environment at the central wavelength of the K-band; that value, together with its error bar, acts as a constraint in the interferometric fit, but it is kept as a free parameter during the fitting process for visibility (see Table 2).

Spectral dependence
As described in Sect.3, our GRAVITY FT data contain six visibilities and four closure phases for each of the five spectral channels.We use this spectral information to derive the spectral index of the circumstellar environment and model our visibility with the following formula with λ 0 = 2.15 µm the wavelength of the central spectral channel of the FT, and k the spectral index defined by where k s and k c are the spectral indices of the star and the circumstellar environment, respectively.For this interferometric data fit, the star photosphere is approximated as a blackbody at the stellar effective temperature over the K-band.While differences between a blackbody and a Kurucz model are large in the ultraviolet domain, in the nearinfrared range the differences are small for the effective temperature range of the star between 10000 K to 20000 K: the photospheric absorption lines across the K-band are scarce and shallow, the spectral slope of the continuum is almost identical to the blackbody, the absolute flux levels can differ by 15-30% but only the slope is relevant when fitting our visibility data.Regarding the spectral indices, k s is derived from the effective temperature of the star (see Sect. 4.2) using Table 5 of Pecaut & Mamajek (2013).k c is an additional free parameter for our geometric models and it can be converted into an equivalent temperature of the dust T d , assuming a blackbody emission.

Fitting processes
The interval between different observing epochs is shorter than 1 month for most of our targets, and about two months for three targets (HD 144668, HD 150198, and HD 158643).For HD 259431, the two observing sequences are spaced out by close to a year but our fit is entirely dominated by the 12 observations obtained in 2018.For HD 144668, for which we have seven and nine observing files, respectively, we have checked that the fits lead to consistent sizes and orientation when considering each dataset independently.Since there is, so far, no observational evidence of temporal variability at the spatial resolution of our dataset, we have combined the different epochs -when available -for each object in our survey prior to the fitting process.This provides us with a statistically robust dataset against the number of free parameters, assuming that the near-infrared emission and the disk structures do not vary.
For the purposes of consistency, the observational data are fitted independently following two different numerical approaches: one based on Markov chain Monte Carlo (MCMC) and one based on a Levenberg-Marquardt algorithm.Since both methods provide results consistent within 3-σ, in this paper, we only present those of the MCMC approach.
We use the same tools as those developed for the PIONIER Large Program (see L17 for details) based on two steps of the Shuffled Complex Evolution (SCE) algorithm (Duan et al. 1993), and then one final MCMC (Haario et al. 2006).The circumstellar environment can be modeled either as an ellipsoid or as a ring.The radial brightness distribution of both models varies from a Gaussian profile to a Lorentzian one (through the Lor parameter; see Appendix D).With such models, the brightness distribution is centro-symmetric and the resulting closure phase is null.The ring model (with a half-flux radius of a r ) is based on a wireframe image convolved by an ellipsoid kernel of semi-major axis a k and the same axial ratio cos i as the wireframe.This permits the description of several kinds of rings, from a thin ring (for very small values of a k with respect to a r ) to a very wide ring tending to ellipsoids (for large values of a k with respect to a r ).In addition, the ring model can include an angular first-order modulation involving sine terms to mimic an inner rim of a disk that is viewed as inclined and that appears as an ellipse with a brighter far side.Such skewed ring models were first introduced to attempt to account for non-zero closure phases in the first IOTA observations (Monnier et al. 2006).Indeed, for such a model, the brightness distribution is not centro-symmetric and non-zero closure phases can be computed.
Since the error bars on our observables can be underestimated and/or correlated, we impose floor values to the error estimates (from the reduction pipeline) of respectively 2% for the squared visibilities and 1 • for the closure phases.Even so, the reduced χ 2 r is generally larger than unity; assuming this results from under-estimated error estimates for the observed quantities, we rescale the χ 2 supplied to the SCE and MCMC algorithms such that χ 2 r ≈ 1, correspondingly increasing the error estimates for the derived model parameters.The reader can refer to Sect.3.3 of L17 for details.While this procedure addresses the issue of optimistic error bars on the observables, there are cases where the issue is clearly not the value of the model parameters, but the model itself.
The free parameters of our models are gathered in Table 2.The ellipsoid model has 7 free parameters; the ring model without azimuthal modulation has 8 and the ring model with a first mode azimuthal modulation has 10.Weighting for radial brightness distribution of convolution kernel (a) log a Log of half-flux semi-major axis (with a in mas) log(a k /a r ) Log of ratio between angular radius of kernel and of ring c 1 , s 1 Cosine and sine angular modulation of 1st order (a) Lor ranges between 0 and 1; Lor = 0 for a Gaussian radial distribution and Lor = 1 for a Lorentzian one.See Appendix D.

Results
As the closure phase signals are generally small, we first focused on centro-symmetric models to derive characteristic sizes for the environment.For all targets, we tested three different models: (a) an ellipsoid leaving Lor free, (b) an ellipsoid with Lor fixed to zero for a Gaussian distribution, (c) a ring leaving Lor free.
Finally we tested models without halo by fixing f h to zero.We compared the models by comparing the χ 2 r values.

General findings
The best agreements are obtained for models with free Lor.Gaussian models with halo generally lead to slightly worse χ 2 r while Gaussian models without halo exhibit much worse χ 2 r than the other models and, thus, they are not considered further.For comparison, Table 3 summarizes the parameters for the best ring models with free Lor and Table C.1 summarizes the parameters for the best Gaussian ellipsoid models.For 21 targets (78%), the best-fit angular sizes are similar among all models.For the other 6 targets (HD 37806, HD 97048, HD 98922, HD 114981, HD 139614, R CrA), since the models with the free intensity distributions are better than the Gaussian ones, we discuss in the following the best-fit parameters of Table 3.For all models, the inclinations i and the position angles PA are well-constrained and do not depend on the models within an accuracy of ± 5 • .
Our models do not fit the data properly for three targets (HD 95881, HD 114981, and HD 142666) since the best models have χ 2 r larger than 15.We retain them in our sample despite the poor fits.

Halo contribution
When fitting the environment with a Gaussian intensity distribution, 70% of the models require a halo contribution (Table C.1.): the halo contribution reaches 10% or more for 8 objects (30%), varies between 5% and 9% for 5 objects (19%), and is below 5% for 6 objects (22%).In contrast, when fitting with a free Lor, the best fits generally converge towards models with null halo contribution, except for HD 85567, HD 97048, HD 142666, HD 169142, PDS 27, and V1818 Ori (Table 3).
The halo contribution determination requires very short baselines, as reported in Setterholm et al. (2018).These are generally not covered in our (u, v) planes but the halo contribution in the H-band as constrained by the PIONIER observations is null or below 6% for all of our targets except HD 100546 (11%) and HD 169142 (8%) (see L17).This halo contribution is expected to be no larger in the K-band than in the H-band when interpreted as scattered light (Pinte et al. 2008).

Inner disk morphology
Using our FWHM values and Gaia distances (Table 1), we derived the half-flux radii (Fig. 3-right) of the environments; these range from 0.1 to about 6 astronomical units (au) with a median of 0.6 au.
Comparing the results obtained for the ellipsoid and the ring models shows that rings generally have a lower χ 2 r .The best fits for the ring are, most often, wide rings, that is, with a ratio of width to half-light radius w = a k /a of the order of unity.Considering all the targets, the median of the relative width w equals 0.83 (Fig. 3-left) and the error bars on ring width are generally large.Thus, we cannot derive an accurate width of the ring and observe an inner cavity, meaning one without any dust inside a Table 3. Best-fit parameters for ring models with 1-σ error bars as defined in Table 2. Closure phases are not included in the χ 2 value.Column 1 corresponds the numbering of Table 1.Column 8 provides the half-flux diameter (FWHM = 2 a) and column 10 provides the dust temperature at 2.2 µm derived from the spectral index k c (see Sect. 5.4).
Object sharp inner rim.This result is in agreement with previous PI-ONIER findings, where only a few objects show inner cavities, even when higher spatial frequencies are probed.
Regarding the disk inclination, the fitted values of cos i range from about 0.3 and 0.95 (i.e., i between ∼ 20 • and ∼ 70 • ), confirming that pole-on assumption (i = 0 • ) for characteristic size determination is not realistic for Herbig Ae/Be stars, despite the low optical extinction usually invoked.

Dust temperature at 2.2 µm
Using the spectral dependence of our interferometric data and assuming the spectral index of the stellar photosphere from the star spectral type, we fitted the spectral index of the environment and converted it into dust temperature at 2.2 µm under the following assumptions: (a) the circumstellar contribution is dominated by the dust radiation while the circumstellar environment produces an emission like a black (gray) body, and (b) there is no non-photospheric unresolved emission.A median dust tempera- ture of ∼1470 K with a standard deviation of ∼200 K is measured (Fig. 4).This is comparable with the sublimation temperature of silicates (1500 K) and of carbon (1800 K) compositions.It could be compared with the PIONIER findings that favor sublimation temperatures which are definitely larger than 1500 K for the dust responsible for the H-band emission (see Sect. 6.2).

Which asymmetry do our non-null closure phase signals probe?
Most objects of our sample exhibit non-zero closure phase differences, indicative of an asymmetry of the environment.As seen in Fig. 2, while about 70% of our targets have closure phase variations over the spatial frequency range smaller than 15 • , a few stars exhibit strong closure phase variations (i.e., larger than 25 • ).This is the case for HD 45677, HD 98922, HD 144432, HD 144668, HD 179218, and R CrA which have an inclination i as large as 60 • .We explored various scenarios that could produce such closure phase signals.

Inclination effects
We first tried to quantify the closure phase signals induced by a projection effect which is only due to an inclination that is different from face-on.For this purpose, we used radiative transfer models (Klarmann 2018) to simulate closure phase signals with different grain populations.We fixed the inclination to 45 • and had the grain size vary between 0.1 µm and 1 mm.At maximum, closure phases of ± 10 • are obtained.Looking at our histogram of closure phase variations (Fig. 2), we can conclude that inclination effects with some grain size distributions (see Sect. 6.2) can explain the closure phase variations below 25 • (81% of our sample).For the highest closure phase signals, however, other scenarios should be invoked.

Azimuthal modulation
In order to mimic the environment's macro-structures as hot spots or arcs, or a binary companion (for 5 among the 6 targets exhibiting a high closure phase variations, Vioque et al. ( 2018) mark them with a Bin flag), we used the same MCMC approach as previously, and a ring model with an azimuthal modulation described in L17.We took into account the closure phase signals for the fit and the χ 2 r computation.As expected, the higher the closure phase signals, the stronger the χ 2 r improvement when considering an azimuthal modulation: this improvement reached a factor 14 for HD 45677, 21 for HD 98922, 6.5 for HD 144668 and 3 for HD 144432.The best-fit models exhibit a brighter, more inclined, and flattened rim as displayed in Appendix D of L17.Two targets (HD 45677 and HD 144668) are best modeled with a skew orientation as expected from a pure inclination effect; these targets have a high inclination of 50-60 • .Two others (HD 98922 and HD 144432) whose inclinations are smaller (30-40 • ) are best modeled with a skew angle different from the PA.For HD 179218, the fit does not converge properly, which is not surprising looking at the "flat" visibility curve (Fig. B.5) and the findings of Kluska et al. (2018): the latter authors combine optical, near-and mid-infrared high resolution data and find that the near-infrared emission radius is about 30 times larger than the theoretical sublimation radius for a dust temperature of about 1800 K.The authors invoke quantum heated dust particles (e.g.polycyclic aromatic hydrocarbons grains) at large radii to explain this unusual extended emission with such a high temperature.Azimuthal modulation could explain the highest closure phase signals.

Size-luminosity relation
The size-luminosity diagram explores the correlation between the location of the K-band emission and the strength of the  2018).The numbering refers to Table 1.
radiation field of the central star.The current picture suggests a strong R ∝ L 1/2 correlation at near-infrared wavelengths, which can be explained by the fact that the emission in this spectral range is dominated by the dust located close to the dust sublimation radius, whose position itself is primarily constrained by the luminosity of the central star (Monnier et al. 2005).This correlation invokes the model of an "optically thin" inner gaseous cavity between the star and the inner radius of the disk.Additionally, to explain the deviation of undersized disks from the correlation law, an alternative model of an "optically thick" inner gaseous cavity is proposed, in which the optically thick gas may shield the dust from direct stellar irradiation, thus bringing the sublimation radius closer to the star.
Bringing our GRAVITY observations together with the Gaia DR2 distances (Table 1), we increase, by a factor of two with respect to the latter authors, the homogeneous sample of K-band disk sizes and significantly improve the coverage of the ∼10 3 L region (Fig. 5).Despite some scatter, we confirm the R ∝ L 1/2 correlation around 10 1 L .In the 10 2 -10 4 L region (early A to late B spectral type), we observe a larger scatter of the data points around the L 1/2 slope, which does not clearly favor either the "optically-thick" or the "optically-thin" inner cavity scenario.Such a scatter is in clear agreement with the PIONIER findings (L17).
In the simple model of a passively irradiated disk with an optically thin inner cavity, the dust radius follows (Dullemond et al. 2001) where T g is the temperature of the dust grain with a cooling efficiency , and R( ) is the stellocentric radius of the K-band emission (which can be different from the sublimation radius).
Here, we assumed a backwarming coefficient equal to one (Kama et al. 2009).In Fig. 5, our measurements can be compared to typical trend lines plotted for boundaries of 1300 K and 1700 K in temperature as well as 0.1 and 1 in cooling efficiency.Using the same assumption for the disk model, we estimate the dust grain cooling efficiency for each object of our sample by taking into account the results of Table 3 (see inset of Fig. 5).The majority of our sources exhibit a dust grain cooling efficiency between roughly 0.1 and 1.Although cooling efficiencies larger than one are also possible, in particular for carbon-poor silicate dust and especially forsterite, this is found more at the mid-infrared rather than at near-infrared wavelengths.Hence, in our context, a cooling efficiency larger than one may simply indicate that the assumed disk model is incomplete and wrong.
A few sources in our sample merit an individual commentary: -R CrA (#26) belongs to the low-quality sample of Vioque et al. (2018).The spread in luminosity reported in the literature for this source is larger than the current error bar, which means that this object could move by some extent towards the right part of the diagram.
-HD 85567 (#5), HD 259431 (#24), and PDS 27 (#25) appear significantly below the = 1 trend line with a rather compact emission size in comparison with their bolometric luminosity.On one side, if we consider their young age of ∼0.04-0.4Ma (see their position in Fig. 1) and their strong accretion rate of ∼10 −6 -10 −4 M /yr −1 (Fairlamb et al. 2015;Reiter et al. 2018), these sources are likely still in the active disk phase where the accretion luminosity strongly contributes to the energy output close to the star.The inner cavity might be filled with optically thick gas, in part detected with GRAVITY, shielding the dust.With these considerations, the assumption that the characteristic size of the emission is to first order set by the luminosity of the central star might not be valid.The viscous accretion produces an emission that is much closer to the star than in the case of the passive disk with an "optically-thin cavity", which results into the measurement of a more compact circumstellar component.
-On the other hand, this must be interpreted in the context of the aforementioned large scatter of the ring radius for high luminosity: indeed HD 95881 (#6), HD 98922 (#8), and HD 190073 (#23) have similar luminosities (log L ∼ 3) and masses, and share roughly the same evolutionary track as HD 259431 (#24) and HD 85567 (#5), but their position in the size-luminosity diagram is quite different, and even opposite to the later ones for HD 190073 (#23) which is above the = 0.1 trend line.Although there are only upper limits for the accretion rate found in the literature for HD 95881 (#6) and HD 98922 (#8) (Fairlamb et al. 2015), HD 190073 (#23) is a strong accretor with ∼10 −5 M /yr −1 yet it is nonetheless found with a large K-band characteristic size.

Can we constrain the grain size distribution of the inner disk?
Even if there is a strong degeneracy between the inner rim position and the grain populations in terms of composition and size, and if the inner rim position defined by the near-infrared emission can be different from the physical inner rim of the dust disk for a given grain population (Klarmann et al., subm.), we can use our GRAVITY measurements to look for trends: -The ring half-flux radii determined by GRAVITY and PIO-NIER (given in Table E.1 in Appendix E) are comparable at a 3-σ level for most of our targets (Fig. 6).The maximum of the K-band emission appears farther away than the Hband one, which is consistent with a K-band emission cooler than the H-band one (see Sect. 5.4).The outliers observed at the shortest radii generally correspond to gapped sources for which the halo contributions are expected to be larger than for ungapped sources, and these can be different in the H-band than in the K-band.As mentioned in Sect.5.2, the halo contribution can only be probed accurately with the very short baselines that are generally not covered in our (u, v) planes.
-The best models in the H and K-bands are wide and smooth rings, and the closure phase distribution (Fig. 2) displays a median value larger than 11 • , which favors wedge-shaped inner rims and disfavors mono-size grain distributions since the latter produce very thin, sharp rings as a result of all grains being sublimated at the same stellar distance.As reported by Klarmann et al. (submitted), wedge-shaped inner rims can be reproduced with grain distributions of two different kinds or more, that is, with grains with a cooling efficiency close to 1 that survive near the inner rim and grains with a smaller cooling efficiency and a temperature close to the sublimation temperature that survive further away.Such distributions, meaning those with sub-micron or micron sized grains, produce maximal closure phases of ± 10 • and a halo contribution at a level of 10-20% due to a high level of scattered light.Adding larger grains or considering a MRN grain size distribution (Mathis et al. 1977) leads to a decrease in the closure phases.-The cooling efficiency depends on the chemical composition of the dust grain (Kama et al. 2009), as well as on the grain size and radiative equilibrium temperature.Grains significantly smaller than the observing wavelength become poor emitters1 , thus leading to a reduced cooling efficiency.According to the location of our sources in the size-luminosity diagram (Fig. 5), the inner disk K-band emission is dominated by grains with low cooling efficiency, although the origin of this low value -the intrinsic opacity of the grain or sub-micron size -cannot be definitely disentangled here.This scenario does not exclude, of course, the presence of a large reservoir of grains with cooling efficiency close to unity -due, for instance, to micron-sized (or larger) grains -but this population would not be sensed by GRAVITY if they already had become sedimented towards the optically thick mid-plane.

Can we probe the inner disk morphology?
In the Meeus classification (Meeus et al. 2001), disks of Herbig Ae/Be stars are classified as group I for flared disks and group II for flat disks.Maaskant et al. (2013) used spatially-resolved midinfrared imaging to suggest that group I disks are preferentially gapped structures, with the resulting large puffed-up wall at the inner edge of the outer disk making the comparative size of the mid-infrared emission larger in relation to the flat, preferentially continuous, disks of group II.Using mid-infrared interferometry, Menu et al. (2015) find that the presence of a gap is not a unique feature of group I disks but may also be present in group II disks, albeit with a smaller width.Looking at Fig. 5 from the perspective of a group-I/group-II dichotomy may graphically suggest that, for a given luminosity class, the group-I (blue symbols) characteristic sizes generally appear larger than those in group-II (red symbols), particularly in the 10-100 L range.However, this interesting hypothesis, which would have important implications in terms of disk structuring, is difficult to confirm beyond any reasonable doubt.Future work on the characteristic of disk size and infrared classification of these objects may shed more light on this area.
From the 27 objects in our GRAVITY sample, 19 of them have a mid-infrared disk size counterpart measured with MIDI.Since the mid-infrared sizes reported by Menu et al. (2015) are the result of modelling face-on disks and neglect inclination effects, we compute an a posteriori correction factor that is to be applied: first, we check their reported (u, v) spatial frequency planes to verify if a specific position angle has been preferentially probed by the MIDI baseline.When the (u, v) coverage is qualitatively uniform, we propose to correct the published midinfrared size by the multiplicative factor f = 2 -cos i, with cos i measured with GRAVITY.For two particular sources, HD 95881 and HD 142666, the limited MIDI (u, v) coverages indicate a preferential position angle along the ellipse minor axis, we apply a multiplicative factor of f =1/cos i (Table 4).
In the size-size diagram of Fig. 7, we compare the K-band and N-band sizes for sources classified as group I or group II according to Meeus et al. (2001);van Boekel et al. (2005); Juhász et al. (2010).In the first order, the characteristic sizes of the  K-band and N-band increase proportionally.At the small-scale level (i.e., <2 au in K-band and <8 au in N-band) there is no strong segregation between the groups I and II, and an overlap is observed, as found by Maaskant et al. (2013).However, at a larger scale, the group-I sources appear more clearly on the upper part of the diagram with respect to the group-II sources, as demonstrated by the plotted linear regression laws for the two populations.Fig. 8. N-to-K ratio as function of age and mass.The triangle symbol shows group-II sources, diamond symbol shows group-I sources and the square symbol shows objects with no classification (i.e., HD 45677).In the left plot, the different colors correspond to different luminosity classes as shown in the inset.In the right plot, blue diamonds identify group-I sources and red triangles group-II sources.Identified gapped sources by high angular resolution observations are marked with (g).The object R CrA (#13 in Table 4) is not represented here as no age estimate is reported in Vioque et al. (2018).

Evolution of the inner disk structure
An interesting question underlying to the group I/group II disk classification is the possible evolutionary link between these two different populations.Maaskant et al. (2013) advance the hypothesis that group I (presumably gapped) and II (presumably continuous) disks may not be connected via a common evolutionary path, but represent, rather, two spatial configurations for independent disks present in different objects at a similar age.Menu et al. (2015) also propose an alternative evolutionary scenario for group I and II disks: either each group might follow a distinct evolutionary path from ungapped to gapped disk, or flat gapped disks could later evolve into flared disks with larger gaps.Recent ALMA results provide clear evidence that gaps in disks -whatever the origin of their formation -are not only a feature specific to the inner disk region but one that can also be found in a concentric arrangement out to large distances (Zhang et al. 2018).Since the one-gap scenario might need to be revisited, the same may occur with the original group I/II classification.
We use the absolute ages and masses provided by Vioque et al. (2018) (Table 1) to study the relationship between the age of the system and the ratio of the measured FWHM values in the N and K-bands (N-to-K ratio).
-From Fig. 8-left, we see that all the sources with a N-to-K ratio larger than ∼10 are classified as group-I (in the sense of Meeus 2001), except for HD 142666.However, this last source is found by Schegerer et al. (2009) to have a gap of ∼0.5 au, which makes its position in our diagram relatively coherent.The N-to-K ratio could be used in part as a proxy to assess the group-I classification.-From the plots of Fig. 8, a marked distinction in the N-to-K ratio appears when looking at the luminosity class, and, correspondingly, that of mass.The brightest, and most massive objects (orange and red colors) show a ratio less than ∼10, whereas objects with lower luminosities (blue to cyan) clearly have a stronger scatter in the N-to-K ratio with values reaching above ∼100.
-With no further consideration of the spectral type, we observe for our sample an increasing scatter in the N-to-K ratio for sources older than ∼1 Ma.
The distribution of group I and group II sources does not clearly point toward a segregation between the two populations based on their "absolute" age or mass, with both types of objects found among ranges in question.However, there is a mass bias in the age determination of pre-main sequence stars through stellar interior modeling.As a consequence, in order to correctly interpret Fig. 8, we should focus more on specific iso-mass paths (see Fig. 1) along which the "relative" ages are more relevant.Looking at the targets at about 2 M (see the dotted box in Figs. 1  and 8-right), we see that the more-evolved sources (in this case group I) closer to the Zero-Age Main-Sequence (ZAMS; Fig. 1) exhibit a larger N-to-K ratio (Fig. 8-right), whereas the lessevolved sources (in this case group II) show a smaller N-to-K ratio.Interestingly, for the objects of about 1.5 M in Fig. 8right, we also find that (with the exception of HD 139614), the more-evolved sources toward the ZAMS (in this case group II) display a larger N-to-K ratio than the less-evolved objects (in this case group I).Note that for 1.5 M all our sources are reported as gapped disks.The limited size of our sample makes it difficult to propose a global evolutionary picture for now, but such a plausible and important trend which would imply an increase of the N-to-K ratio with age should be confirmed through the undertaking of additional measurements of K-and N-band sizes aimed at populating the plots of Fig. 8 with, for instance, all the sources of the iso-mass 5 M in Fig. 1.

Gap formation scenarios
Photo-evaporation (Armitage 2011) is one of the key processes of disk dispersal via the formation of gaps and disk inner holes.Heating by the central star through radiations in the Extreme Ultra-Violet (EUV), Far Ultra-Violet (FUV), and in the X-rays induces the splitting of the disk and the formation of a gap with a peak efficiency at the critical radius r c ∼ 0.1-0.2r g (Gorti et al. 2009), where r g = GM * /c 2 s is the radius where the heated gas in the disk becomes gravitationally unbound with regard to the central star.In the EUV/FUV/X-ray photo-evaporation scenario, the formation of gap takes ∼4 Ma (see Fig. 8 in Gorti et al. (2009)) for solar-type masses, and is likely shorter for masses larger than 3 M (see Fig. 12 of Gorti et al. (2009); Ercolano & Pascucci (2017)).Once the gap is formed, the inner disk is rapidly depleted in a very short viscous timescale of ∼10 5 years (Gorti et al. 2009).Over the course of this interval, the result is the formation of the transitional disk with a void inner cavity.
Using our GRAVITY observations, we can directly compare the characteristic size of the K-band emission to the theoretical critical radius (Gorti et al. 2009) where T is the temperature of the EUV/FUV/X-rays heated gas.
With the exception of HD 45677, all our sources lie below the critical radius for a standard T ∼10 4 K temperature (Fig 9).When looking at our gapped sources in Fig. 8 with an age of ∼5-10 Ma, and assuming the aforementioned photo-evaporation scenario, we could have expected these objects to already have been in the transitional stage with the inner rim of the disk beyond the critical radius, but this is not the case according to Fig. 9.According to this model, EUV/FUV/X-ray photo-evaporation would have had a negligible or no effect for these disks.Alternatively, one may surmise that the timescales for photo-evaporation are significantly longer as found for EUV photo-evaporation, or that the process of gap formation in these discs is indeed dominated by the dynamical clearing of young planets.For the objects with an age of ∼1 Ma and younger, the smaller sizes compared to r c are compatible with a scenario where photo-evaporation had just started.Finally, given the approximations in Eq. 5, the position of HD 45677 with respect to the r c slope remains uncertain.However, when combining the different graphic results of Sect.6, this source appears consistent with the evolutionary stage of a young, flat and ungapped disk.This could be corroborated by photometry and SED analysis following Meeus et al. (2001) and van Boekel et al. (2005).

Inner and outer disk (mis)alignment
Connecting the inner and outer parts of protoplanetary disks are essential for understanding all the dynamic effects involved, as well as the implications for planet formation in particular.While the recent images of scattered light produced by SPHERE and images from ALMA in the millimetric range have allowed us to determine orientations for the outer disks accurately, nearinfrared interferometry is a relevant means to measure these inclinations and these position angles of the inner disks.
We aim to compare our inner disk orientations with the outer disk orientations derived from imaging using either SPHERE or ALMA.Since the ALMA and SPHERE images have a variety of spatial resolutions, in order to simplify analysis and avoid the over-interpretation of small variations between inner and outer inclinations and position angles, we determine a uniform uncertainty on these values that is to be applied to all objects.We generate a family of synthetic observations of thin disks with known inclination and position angles for the nearly edge-on (∼80 • ) and nearly face-on (∼10 • ) mimicking the lowest spatial resolution data in our target sample [i.e., ALMA data of HD 97048 with ∼0.5 resolution, (van der Plas et al. 2017)].For both the nearly face-on case (5 • to 15 • ) and the edge-on one (75 • to 85 • ), we vary the inclination in steps of 1 • .After convolving the model thin disks with the ∼0.5 beam and producing the synthetic observations, we fit Gaussians to the resulting images to recover the major and minor axes and calculate the inclination and position angles.The inclination angles of both the edge-on and face-on cases can be reliably recovered to a level within 4 • .While the position angles of the edge-on case can be recovered to within 2 • , the face-on case carried more uncertainty, with only a 5 • accuracy.Since most of our sources are closer to face-on than edgeon, we adopt a uniform and conservative 5 • uncertainty for the inclination and position angles of all sources.Since the SPHERE and ALMA orientations are consistent with each other amid our conservative error bars of 5 • , we consider the average between the two determinations.
We can compare these orientations for five targets (the orientations for HD 169142 using GRAVITY are too poorly constrained).We note differences for HD 135344B and HD 142527 (Fig. 10).Interestingly, recent images in polarized light of the latter star have revealed features like spirals and dark arcs (Avenhaus et al. 2014).A scenario that has been invoked to explain such features is attributed to shadows cast by an inner disk that is misaligned with respect to the outer disk (Marino et al. 2015).The authors cited use radiative transfer models to reproduce the NACO observations and they derive a relative inclination between the inner and outer disks of 70 ± 5 • and a position angle of the inner disk of -8 ± 5 • .While we determine a small position angle (∼ 10 • ), our near-infrared interferometric observations are not consistent with such a high difference in inclination.It is out of the scope of this paper to address the open questions on the dust structures of this transitional disk (Avenhaus et al. 2018).Increasing the sample for comparing the orientations of the inner and outer disks, especially with regard to targets for which arcs, warps, and shadows have been observed and or for which planets have been detected, could be a key to a better understanding of the interplay between the different components of these complex and dynamical environments, and of the origins of the variety of features seen by disk imaging as well.

Summary and perspectives
Near-infrared long baseline interferometry is a powerful means to constrain the continuum emission due to the effects of reprocessed light by the dust grains in the circumstellar environments of young stellar objects.We extend the PIONIER study of Herbig Ae/Be stars to the K-band through a resolution with the GRAVITY instrument of the inner disks around 27 targets of a sample that spans a large range of luminosity, mass, and age.The environments of these stars in the near-infrared continuum appear as wide and smooth rings, exhibiting an asymmetry since the measured closure phases are most often non-zero.To probe the disk morphology, a comparison of the K-band FWHM with N-band ones could be used as a proxy to disentangle flat and flared disks since a transition around the N-to-K ratio of about 10 is observed.The brightest, most massive stars show a ratio less than 10, whereas objects with lower luminosities clearly have a stronger dispersion in the N-to-K ratio with values up to 100.As far as disk evolution is concerned, no universal mechanism is clearly in evidence since the different types of disks (flat or flared, with or without gaps) are found ranging from 0.2 Ma to 15 Ma.In looking at the inner positions derived from our GRAV-ITY measurements, it appears that gap formation processes in the disks are likely related to the dynamical clearing of young planets rather than FUV/EUV/X-ray photo-evaporation.Near-infrared long baseline interferometry can also be used to derive the inner disk orientations that are to be compared with those of the outer disk so as to better understand the origins of the variety of features seen through disk imaging in these complex and dynamical environments.
As perspectives, investigating macrostructures that have proven more complex physically using radiative transfer codes or the presence of undetected close companions to explain the largest closure phases would be a natural follow-up for the current study.Combining GRAVITY and MATISSE observations will be of utmost interest in probing dust distribution and disk mineralogy and understanding the mechanism for disk evolution and the gap formation processes by populating the N-to-K size diagrams with stars of a few solar masses.
. 1.Most of them are part of the PI-

Fig. 2 .
Fig. 2. Histogram of GRAVITY closure phase peak-to-peak variations over spatial frequency range for our sample.

Fig. 3 .
Fig. 3. Left.Histogram of normalized width w for ring model of our GRAVITY observations.Right Half-flux radius a derived from our GRAVITY measurements and Gaia DR2 distances.Insert is a zoom on the shortest radii.

Fig. 4 .
Fig. 4. Bar plot of temperatures of environment dust at 2.2 µm, assuming a blackbody emission.Dashed line shows the median dust temperature.

Fig. 5 .
Fig.5.Ring half-flux radius vs. luminosity for group I (blue), group II (red) and unclassified (black) targets.The dashed and dotted lines show the relationship for the labeled temperature and cooling efficiency of the dust g of an optically thin inner cavity model.The inset shows the distribution of values of the dust cooling efficiency derived from our ring half-flux radii and dust temperatures.The empty symbol (source #26) corresponds to R CrA, a low-quality source inVioque et al. (2018).The numbering refers to Table1.

Fig. 6 .
Fig. 6.Comparison of half-flux radii of ring models in H and K-bands as determined by PIONIER and GRAVITY observations (see TableE.1).Blue diamonds denote group-I sources, while the red triangles denote group-II sources.The gapped sources are marked with (g).The dash-dotted line is the resulting linear regression.The inset is a zoom on the smallest ring radii.

Fig. 7 .
Fig. 7. Size-size diagram of FWHM in N-and K-bands for group I (blue), group II (red), and unclassified (black) sources.The numbering refers to Table 4. Blue and red dashed lines correspond to an error-barweighted linear regression on the group I and group II sources, respectively: y = ax + b with a I =5.87±0.14, b I =-0.50±0.08, a II =3.23±0.04,b II =1.54±0.09.The dotted lines delimit the uncertainty on the slope parameters.The inset shows a zoom of the overlap region (see text for details).

Fig. 9 .
Fig. 9. Half-flux ring radius as function of mass of central star.Dashed line corresponds to the critical radius r c given by Eq. 5 as T = 10 4 K.

Fig
Fig. B.1.GRAVITY observations: visibilities squared (left), closure phases (middle), (u, v) plane (right).The colors code the spectral channel of the FT.The symbols without contours at the bottom of the visibility and closure phase curves display the residues of the best-fit ring models without azimuthal modulation.

Table 2 .
Free parameters of our geometric models.

Table 4 .
Menu et al. (2015)characteristic sizes.Near-IR sizes correspond to half-flux diameter ring model.Mid-infrared data are the halfflux diameter adapted fromMenu et al. (2015).Multiplicative factor f c is applied to the measured N band characteristic diameter.

Table A .
1. Observation log of VLTI/GRAVITY observations.N denotes the number of 5-minute long files that have been recorded on the target.