© ESO, 2017

1. Introduction

Herbig Ae Be stars are intermediate-mass (2−8 M_⊙) pre-main sequence stars. They are surrounded by accretion disks that are believed to be the birth places of planetary systems. The remarkable diversity of discovered exo-planetary systems has renewed the strong interest in understanding the structure of such disks at the astronomical unit (au) scale. The sublimation front at the inner rim contributes to define the physical conditions in the dusty planet-forming region.

Hillenbrand et al. (1992), Lada & Adams (1992) have recognized that the spectral energy distributions (SED) of Herbig Ae Be stars display near-infrared excesses in the 2–7 μm region that are too strong to be reproduced with the flaring hydrostatic disk models that are successfully applied to T Tauri stars (see e.g., Kenyon & Hartmann 1987; Chiang & Goldreich 1997). The initial proposition by Hillenbrand et al. (1992) that accretion heating could provide the near infra-red (NIR) excess flux can be ruled out on the arguments that, firstly, such accretion rates are not measured and secondly, they are incompatible with the presence of an optically thin inner cavity (Hartmann et al. 1993).

The first survey of Herbig Ae/Be stars using near-infrared optical interferometry by Millan-Gabet et al. (2001) concluded that the near-infrared emission size was indeed incompatible with the standard flared disk models. In order to explain both interferometric and photometric measurements Natta et al. (2001) have proposed that the rim vertical structure had to be taken into account in the global stellar energy reprocessing balance. Such a rim is frontally illuminated and therefore naturally hotter. It should emit at the dust sublimation temperature and puff up (Dullemond et al. 2001). Aperture-masking observations of a massive young star by Tuthill et al. (2001) revealed a doughnut shape NIR emission and offered a first confirmation of this idea. Monnier & Millan-Gabet (2002), Monnier et al. (2005) assembled the available interferometric data to confirm that the measured sizes were well correlated with the central stars’ luminosities and consistent with the dust sublimation radius. This trend was observed over more than four decades of stellar luminosity. The advent of spectrally-resolved interferometry has further revealed the presence of an additional hot component emitting from within the sublimation radius in a significant sample of objects (Malbet et al. 2007; Tatulli et al. 2007; Kraus et al. 2008b,a; Isella et al. 2008; Tannirkulam et al. 2008b; Eisner et al. 2009, 2014).

Several teams have proposed alternative explanations for the 2–7 μm excess: an extended (spherical) envelope (Hartmann et al. 1993; Vinković et al. 2006), disk wind (Vinković & Jurkić 2007; Bans & Königl 2012), magnetically lifted grains (Ke et al. 2012), and supported disk atmosphere (Turner et al. 2014).

The availability of spatially resolved observations of the inner rims has prompted research on the processes that might structure them. Isella & Natta (2005) have shown that the dependence of dust sublimation temperature on gas density causes the sublimation front to be curved. Tannirkulam et al. (2007) pointed out that the dependence on grain size of cooling efficiency and dust settling results in a broader NIR emission region. Kama et al. (2009) and Flock et al. (2016), taking into account the pressure dependence of the sublimation temperature, but not grain settling, find rim profiles similar to Tannirkulam et al. (2007), that is, a wedge pointing to the star. Vinković (2014) addressed the issue of the thickness of the rim; considering a comprehensive list of processes such as viscous heating, accretion luminosity, turbulent diffusion, dust settling, gas pressure, he concluded that the resulting inner rim thickness was too small to explain the infrared excess. The reader is referred to Dullemond & Monnier (2010) or Kraus (2015) for reviews of the structure of the inner circumstellar disks.

The goal of this article is to provide a statistical view of the milli-arcsecond morphology of the near-infrared emission around Herbig AeBe stars. For that purpose, we took advantage of the improved performance offered by the visitor instrument PIONIER (sensitivity, precision and efficiency) at the Very Large Telescope Interferometer (VLTI). We established a Large Program (090C-0963) that had the following questions in mind:

• 1.

  Can we constrain the radial and vertical structure of the innerdisk?

• 2.

  Can we constrain the rim dust temperature and composition?

• 3.

  Can we constrain the presence of additional emitting structures (e.g., hot inner emission, envelopes)?

• 4.

  Can we reveal deviations from axisymmetry in the emission that could be linked to radiative transfer effects?

This article is organized as follows. Section 2 introduces the sample and the immediate data reduction, that produces calibrated data. In Sect. 3, we submit the observed data to simple model fits, in a systematic and uniform way, condensing hundreds of visibilities and phase closures into a dozen or fewer parameters for each object. In Sect. 4 we use the processed data to infer some of the physical properties of the sample objects, striving to avoid possible biases. In the conclusion, we wrap up our main findings, place them into perspective with current models, and consider the outlook for follow-up work.

2. Observations

2.1. Instruments

The interferometric observations were made using the four 1.8 m auxiliary telescopes (AT) of the VLTI (see Mérand et al. 2014) at ESO’s Paranal observatory, and the PIONIER instrument (Le Bouquin et al. 2011), a four-telescope recombiner operating in H-band (1.55–1.80 μm). The three AT configurations available at the epoch of the observations provided horizontal baselines in the range 11–140 m. The photometric observations were obtained with two telescopes at the South African Astronomical Observatory (SAAO) Sutherland Observatory: 0.50 m for UBVRI bands, and 0.75 m for JHK bands.

2.2. Sample

Our sample is based on the lists of Malfait et al. (1998), Thé et al. (1994), and Vieira et al. (2003). One should note that a) the definition of HAeBe objects does not rest on quantitative criteria; b) our list of targets does not meet any completeness criterion. Our sample is limited at , where the fainter objects are intended to test the sensitivity limit of PIONIER/VLTI. The program objects are listed in Table 1.

Spectral types.

They were obtained from the SIMBAD database and checked against the Catalog of Stellar Spectral Classifications (Skiff, 2009–2016)¹; in case of discrepancy, we used a recent measurement from the compilation of Skiff, and we quote the original reference. When available, we adopted the determination by Fairlamb et al. (2015), that also includes determinations of bolometric luminosity (the values of T_eff were converted to spectral types using the tables of Pecaut & Mamajek 2013). Statistics of H magnitude and spectral types are shown Fig. 1.

Distances.

Parallactic distances were obtained from the Gaia (Lindegren et al. 2016) and Hipparcos (Perryman et al. 1997) catalogs, and, failing that, from detailed spectroscopic studies; see notes of Table 1.

Binarity.

For T Ori see Shevchenko & Vitrichenko (1994); for HD 53367 see Pogodin et al. (2006); for HD 104237 see Cowley et al. (2013); for AK Sco see Andersen et al. (1989); concerning VV Ser, Eisner et al. (2004) mention that their data do not rule out a binary model; binarity is not supported by our (more extensive) data. HD 36917 is often quoted as SB2, referring to Levato & Abt (1976), who actually mention a composite spectrum, which is not confirmed by Alecian et al. (2013). HD 145718 is listed by Simbad as an Algol binary. However, following Guimarães et al. (2006) and references therein, this remains unconfirmed.

Coordinates of HD 56895.

We provide the following information for the record, and to avoid others a similar mistake. We list in Table 2 the coordinates for HD 56895 or HD 56895B from various sources. It appears that the object listed by Malfait et al. (1998) as HD 56895B has the coordinates (within a couple of arcsec) and the V magnitude of HD 56895, while its NIR magnitudes are those of a nearby (13′′) object with much redder colors. Our interferometric observations were made at the position of HD 56895, not what was intended. As a consequence, this object is not included in the data analysis.

2.3. Interferometric data

The interferometric data were acquired with the PIONIER instrument at the VLTI, during a total of 30 nights in six observing runs between December 2012 and June 2013, under Large Program 190.C-0963. The log of observations is given in Appendix C.

We used all the available configurations offered by the ATs: A1–B2–C1–D0, D0–G1–H0–I1, A1–G1–J3–K0. Most of the data were dispersed over three spectral channels across the H band, providing a spectral resolving power of ℛ ≈ 15, while the faintest objects were observed in broad H-band (ℛ ≈ 5).

Data were reduced and calibrated with the pndrs package described in Le Bouquin et al. (2011). Each observation block (OB) provides five consecutive files within a few minutes. Each file contains six visibilities squared V² and four closure phases φ_cp dispersed over the three spectral channels. Whenever possible, the five files were averaged together to reduce the final amount of data to be analyzed and to increase the signal-to-noise ratio. The statistical uncertainties typically range from 0.5° to 10° for the closure phases and from 2.5% to 20% for the visibilities squared, depending on target brightness and atmospheric conditions.

Each observation sequence of one of our targets was immediately preceded and followed by the observation of a calibration star in order to master the instrumental and atmospheric response. Le Bouquin et al. (2012) have shown that such calibration star(s) should be chosen close to the science object both in terms of position (within a few degrees) and magnitude (within ± 1.5 mag). We used the tool SearchCal² in its FAINT mode (Bonneau et al. 2011) to identify at least one suitable calibration star within a radius of 3 deg of each object of our sample.

Most of our objects are grouped into clusters in the sky. Consequently, the instrumental response could be cross-checked between various calibration stars. This allowed us to discover a few previously unknown binaries among the calibration stars. These have been reported to the bad calibrator list³. We estimated the typical calibration accuracy to be 1.5° for the phase closures and 5% for the visibilities squared. The typical accuracy of the effective wavelength of the PIONIER channels is 2%. Finally, the on-the-sky orientation of PIONIER has been checked several times and is consistent with the definition of Pauls et al. (2005).

2.4. Photometry

We made dedicated photometric measurements of most objects in our study, in the UBVRIJHK bands, at the SAAO Sutherland 0.5 m and and 0.75 m telescopes (now decommissioned); the time parallax relative to the interferometric observations is smaller (at most two years) than for photometry available in the literature. Whenever SAAO measurements could not be made, previous measurements from the literature were used with weighting following their error estimates; see Tables 3 and 4. The NOMAD catalog is not a primary source, and was used only for B and V magnitudes of MWC 297 and VV Ser from the unpublished USNO YB6 catalog. The Wise measurements are used only for plotting on SED graphs.

3. Data processing

In this section, we present steps of intermediate data processing, from reduced and calibrated data (magnitudes, visibilities, phase closures) to meaningful derived parameters. The first subsection below deals with model fits of photometric data, while the remaining subsections deal with model fitting of interferometric data.

3.1. Model fits of the spectral energy distribution

For each object, we performed a least squares fit of the observed fluxes (actually, the logarithms) in the eight Johnson-Cousins (Bessell et al. 1998) bands U, B,..., K, according to a simple model of a stellar photosphere plus dust emission (modeled as a single temperature blackbody). The maximum of λF_λ, for T = 1500,...,1800 K is at λ ≈ 2.4,...,2.0 μm in K band; the other, shorter λ bands used in the SED fit discriminate against cooler dust. (1)where

• k = 1,...,8;

• (m_V−m_k)_s are the intrinsic stellar colors derived from the published (SIMBAD) spectral type and an interpolation in the tables of Pecaut & Mamajek (2013);

• r_k is the extinction ratio A_k/A_V from the extinction law of Cardelli et al. (1989) with R = 3.

The fit parameters are

• F_sV, the absorption-free stellar flux in V band;

• F_dH, the absorption-free dust flux in H band;

• A_V;

• T_dp, the blackbody temperature of the dust component as inferred from photometry.

Three of these (F_sV, F_dH, T_dp) will be used in the Discussion section; also, f_d, the fraction of the observed flux at 1.63 μm contributed by the dust component, will be used in the interpretation of interferometric data, see Sect. 3.4 below.

What value of the extinction ratio R_V should we adopt? Hernández et al. (2004), in a study of 75 HAeBe objects, provide solid evidence for a value R_V = 5 fitting the data better than the standard value R_V = 3.1. This conclusion rests predominantly on objects with a high visual extinction . On the other hand, Montesinos et al. (2009) conclude that R_V = 3.1 results in a better fit to the SED of VV Ser (A_V ≈ 3) and adopt this value for other objects with lower extinction; Fairlamb et al. (2015) take a similar approach; see also Blondel & Tjin A Djie (2006), where among those 60 objects having A_V ≤ 1.0, the preferred value of R_V is 3.1 in 50 cases. Our sample has an A_V distribution with a median value 0.36; accordingly, we decided to perform our analysis assuming R_V = 3.

The numerical method used for fitting the SED is described in Sect. 3.3 below, and the 1σ error bars are derived from the marginal distributions of the Markov chain. The fit results are displayed in Table B.1. A sample fit is shown in Fig. 2.

3.2. Interferometric data: at first glance

Before embarking into numerical model fitting, we should have a qualitative feel for the relationship between measurement results and object properties. We start with a plot of visibilities for the object HD 100453, see Fig. 3. Following common use, the plot is for visibility squared V² – a shorthand for |V|²–, but we will discuss visibilities V. Assuming that the star itself is unresolved at all observed baselines, and that we observe at a single wavelength λ₀ = 1.68 μm, the combined complex visibility of the star and circumstellar matter is: (2)where V_c(u,v) is the visibility that the circumstellar component alone would have, f_c the fraction of the total flux (at λ₀, inside the diffraction-limited field of view of one VLTI telescope) contributed by that component; likewise for f_s and V_s, with generally V_s = 1.

To first order, and ignoring azimuthal dependency, the variation of the combined object visibility V as a function of the radius in the (u,v) plane , reflects the variation of the visibility of the circumstellar component V_c(u,v). The part of the diagram where V² decreases with increasing b contains information on the characteristic size and even the shape of the circumstellar component, see legend sh in Fig. 3. At large enough baselines, the circumstellar component is fully resolved, the only coherent flux is from the star, diluted by the circumstellar flux: , see Fig. 3. Conversely, if one extrapolates the curve of V² to zero baseline, the intercept may be at a value less than one. Actually, the visibility must return to unity at zero baseline, but it does so at a baseline significantly smaller than those sampled by PIONIER/VLTI, i.e. at an angular scale larger than those analyzed by our observations. We will call this fully resolved component the “halo” without otherwise presuming a specific morphology; see also Lachaume et al. (2003, Sect. 3.4).

The visibility curve is used as an a priori criterion to define objects with high quality data “HQ”: a necessary condition being that the squared visibility drop to 80% or less of its zero-baseline intercept. This, together with binarity and other exclusion clauses, appears in the footnotes to Table 1.

Finally, we note that the three wavelength channels appear to be spread diagonally in the b-V² diagram. When observing with a geometrical baseline B, the three spectral channels are sampling different values of . From that cause alone, the various points would populate a unique curve following Eq. (2). To proceed, we introduce in the visibility equation: (a) the dependence on wavelength; (b) the explicit hypothesis that the circumstellar component’s visibility can be factored into a scalar wavelength term and a wavelength-independent spatial part; and (c) we assume the star to be fully unresolved V_s(u,v) = 1: (3)In the H band the stellar flux rises at shorter wavelengths, the dust flux decreases, therefore f_c decreases, and the combined visibility V increases. We refer to this as chromatic dispersion (noted c.d. in Fig. 3). We note that an anisotropic brightness distribution can also cause a spread in the b-V² diagram, which is different from chromatic dispersion because it manifests itself also in monochromatic data.

When the observations reach clearly into the regime where the circumstellar emission is fully resolved, as in the example we have chosen, it is possible to infer from the values of s in each channel the respective values of f_c. Returning to the other observed values of V², using Eq. (2), and under the assumption that V is real positive, it is possible to isolate the visibility of the circumstellar component V_c(u,v). Because these conditions are rarely all met, we cannot generally analyze a program source on the back of an envelope, and we must resort to parametric modeling.

3.3. Parametric model fitting: general methods

We use parametric models to extract morphological and statistical information from our interferometric dataset. Compared with image reconstruction, this has drawbacks, mainly that a functional form is assumed for the source spatial structure. Only a few of our objects have a dataset rich enough to allow for image reconstruction to be attempted; and image reconstruction itself is not totally free from assumptions – regularization, for example. On the other hand, parametric fitting has the advantage that within a family of models, formal error bars can be assigned to the inferred parameter values. Also, in principle, the adequacy of respective models can be gauged from the χ² that they achieve.

Building upon the observations of Sect. 3.2, we construct a parametric model with the following components: (1) the star, assumed unresolved at all observed baselines; (2) circumstellar emission with a model-specific and wavelength-independent spatial structure; (3) an extended halo component, defined by its spatial scale: fully resolved at all observed baselines (ℓ ≳ 20 mas), but inside the PIONIER field of view, defined by its single-mode coupling to the Auxiliary telescopes (ℓ ≲ 200 mas); in the context of this work, this presumes nothing concerning the shape of that component.

All fitting was performed in Fourier space; sky images were derived only for illustration purposes. Within the limits of the H band, analyzed by PIONIER in three spectral bands (1.59, 1.68, and 1.76 μm), the spectral dependence of either the star or the circumstellar component is modeled as a power law, defined by the respective spectral indices k_s and k_c, where k = dlog  F_ν/ dlog  ν. Since at large distances from the star, the H band is in the far Wien region of the grains’ thermal radiation, we assume that the halo component arises from diffusion of starlight (see e.g., Pinte et al. 2008) and has the same spectral index. We verified that making the alternate assumption (same spectral index as the thermal radiation from the near circumstellar component) has a minimal impact on the fit results. Returning to Eq. (2), but now considering a range of wavelengths and including the halo component as discussed above, the visibility can be expressed as: (4)where λ₀ is a reference wavelength (conveniently that of the central PIONIER channel, 1.68 μm) at which the flux fractions of the three components are defined: f_s, f_c, f_h for star, circumstellar, and halo respectively (f_s + f_c + f_h = 1), and, in contrast with Eq. (3), the wavelength dependence is explicitly parameterized. The visibility of the circumstellar component V_c(u,v), is specific to each model, as concerns both the functional form (chosen) and the numerical values of the parameters (the result of the fitting procedure). Several functional forms are introduced below in Sects. 3.5–3.7.

The value of k_s is found assuming that the object’s central star radiates in H band as a blackbody at T_eff. T_eff itself is derived from the spectral type, using Table 5 of Pecaut & Mamajek (2013).

The 1σ error estimates for the observed quantities, visibility squared and closure phase, are computed by the data reduction software, and are derived from the internal scatter during one exposure (typically on 1 min timescale).

The minimization of χ² proceeds in three steps. First the shuffled complex evolution (SCE) algorithm (Duan et al. 1993) performs a global search for a minimum of χ². Next, the breakdown of the value of χ² is examined, and the contributions from V² and φ_cp are allocated weighting factors, such that (a) they are balanced; (b) the reduced χ² is close to unity. The SCE search is then restarted. Finally the Monte-Carlo Markov chain (MCMC) algorithm (Haario et al. 2006), still with rescaled χ², is used to (a) derive 1σ error bars for each fit parameter from the marginal distributions of the Markov chain; (b) improve slightly the best-fit χ². As in Barmby et al. (2009), the rescaling of χ² is motivated by some error bars on observables being likely underestimated, with a risk of underestimating the error bars on the fitted parameters. Moreover, some of the errors on the observables might be correlated (per-night calibration errors). An over-estimated χ² carries the risk of the minimum search algorithm being caught in a local minimum. The unmodified χ² value is listed in the results as an indication of how the model captures the complexity of the data.

3.4. Size-flux degeneracy

Our initial fits of simple models to the interferometric observations (we skip for the moment some details of the model, to be found in Sect. 3.5 below) revealed, as one might expect, that some parameter pairs are correlated. Apart from the obvious correlation arising from: f_s + f_c + f_h = 1, the clearest correlation was generally found to be between the characteristic size and the fractional flux of the circumstellar component. This is illustrated in Fig. 4. The top panel shows the scatter plot (in the MCMC chain) for the variables a, the half-flux radius of the circumstellar component, and f_c, its fractional flux. Points A and B at 10% and 90% of the respective marginal distributions are also shown; note that the respective values differ by a factor of more than 1.5. In the bottom panel, one can see that the visibilities of two Gaussian models for parameter values at A and B are quite similar within the range of observed B/λ, and could be discriminated only by observations at longer baselines.

We see that a degeneracy between the size and the flux fraction of the circumstellar component plagues the model fitting of a partially resolved object, as already noted by Lachaume (2003). This begs the question: are the dust component (identified in the spectral domain) and the circumstellar component (identified by its spatial properties), related? The correlation in Fig. 5 shows that this is indeed the case. A similar plot of f_c + f_h versus f_c (not shown) results in a correlation that is neither better nor worse.

Seeing that f_c deviates the most from f_d (reminder: the photometric flux fraction of the dust component at 1.68 μm, defined in Sect. 3.1 above) for objects where the former is poorly constrained by interferometric data, we will assume, as a working hypothesis, that f_d and f_c measure the same (circumstellar) component of the objects. Accordingly, in order to account for this additional constraint, we introduce in the expression of χ² that drives the fitting of interferometric data an extra term ((f_c−f_d) /σ(f_d))², where the value of f_d is from the photometric fit. Figure 6 shows the improvement in the size standard error when the photometric constrain is imposed in the fit of interferometric data. Such a photometric constraint is applied in all subsequent model fits.

3.5. Ellipsoids

The first and simplest model that we use for the circumstellar emission has a brightness distribution continuously decreasing outwards, with elliptic isophotes.

The Gaussian radial brightness distribution is widely used in modeling. But it is peculiar in having a strong asymptotic decay, faster than any power law or plain exponential, which might not be a good description of reality. So we also use a distribution with an opposite behavior, having an r^-3 asymptotic decrease (the smallest integer exponent with converging flux integral). That is the pseudo-Lorentzian (Lor) distribution as shown in Table 5, which has a simple analytic Hankel (radial Fourier) transform.

Starting from a face-on model with circular symmetry, an inclined model with arbitrary axial ratio cosi and position angle θ can be derived by simple transforms in the (u,v) coordinates.

In exploratory tests, it was found that some objects “favor” (achieving a lower χ²) the Gauss profile, while other objects favor the Lor distribution. Following that, we introduced f_Lor as one of the fit parameters, allowing intermediate radial distributions, such that the circumstellar component visibility in Eq. (4) is (5)We proceed with the fits of the program objects. Table 6 gives the list of fit parameters; the values are listed in Table B.2.

3.6. Ring

This type of model is intended to model the emission from a bright rim, for example, the inner rim of a sublimation-bounded dust disk. We start with a “wireframe” distribution (6)normalized to unit flux. Its complex visibility is (7)Next we introduce azimuthal modulation: (8)where φ is the polar angle, m is the order of azimuthal modulation, and the sum is absent if m = 0. The underlying motivation for such a modulation comes from models such as those of Dullemond et al. (2001), Isella & Natta (2005), where the inner rim of a disk is viewed inclined: the rim contour appears, to first order, as an ellipse, with the far side being more luminous from the vantage point of the observer. Leaving aside the radial structure, the angular modulation in such a model involves only sine terms of odd order. Because we want our analysis to remain model-agnostic as much as possible, we do not include such a restriction in Eq. 8.

Transforming the modulation amplitude to polar representation: c_j + i s_j = ρ_jexp(i θ_j) and likewise for the coordinates in the (u,v) plane: q_u + i q_v = qexp(i ψ), the complex visibility for the modulated skeleton is: (9)Just as for ellipsoids (Sect. 3.5), the wireframe distribution is squeezed by a factor cosi along the direction of the minor axis, and rotated by the position angle. The above-mentioned azimuthal modulation, Eq. (8), is defined with the angle origin along the major axis.

Because real ring-like images (if they exist) have a finite width, the wireframe image is convolved by an ellipsoidal kernel, with a semi-major axis a_k, the same axial ratio cosi as the wireframe, and a radial distribution that is a hybrid between a Gaussian and a pseudo-Lorentzian, as described above in Eq. (5) and Table 5.

Irrespective of the particular parameterization, this model can describe both thin rings (a_k ≪ a_r) and, asymptotically (a_k ≫ a_r) ellipsoids, since the convolution kernel is identical to the ellipsoid model.

The geometrical (see Table 7) parameter pair (a_r, a_k) becomes degenerate in the limit of ellipsoid-like cases where only a_k is relevant. The physical pair is more meaningful, especially when the object is marginally resolved, and the parameter w spans a finite range 0 ≤ w ≤ 1. It is not, however fully suitable for the MCMC search because (by default) the method assumes a uniform prior over the allowed search interval: most of the 0 ≤ w ≤ 1 interval corresponds to a distinct ring structure, leading to a bias towards detection of such a structure. The “fit” parameters in Table 7 cover a wide range, with a scale-free prior. The fit parameters number one more (excluding azimuthal modulation) than for the ellipsoid model; they are listed in Table 8.

Figure 7 shows examples of fits for some program objects, and shows how a given object – observed at a given resolution – requires a certain degree of complexity, but not more.

3.7. Power law with inner boundary

In some, if not all, current models for inner disk rims, one expects the brightness distribution to have a sharper cutoff inwards than outwards. Accordingly, we tested models with a radial distribution: (10)and similar models with an exponent different from −3, or an exponential decay, with or without convolution by a kernel (see Sect. 3.5 above). None of these showed a significant advantage, as measured by the χ², over the already listed model types.

4. Analysis and discussion

In this section we use the results of parametric fits performed in the previous section (Data processing) to address the properties of the objects themselves. Because the observations provide (at best) a 2D (interferometry) or 1D (photometry) projection of the objects, and because the parametric models are arbitrary, we attempt to reduce these shortcomings by comparisons with physical models.

We issue up front two disclaimers: firstly, our goal is not to fine-tune a physical model and declare that it represents reality; at best we can hope to support statements on some large-scale properties of the real objects. Secondly, we discuss our sample globally, and derive generic properties, which makes the best use of our limited statistics, but does not imply that real objects come from a single mold. Depending on the issue being discussed, the sample may be either the 44 non-binary objects “noBin” or the 27 so-called high-quality objects (HQ): in short, those non-binary objects that are suitably resolved, with V² dropping below 0.8 at the long baselines.

4.1. Dust temperatures

4.1.1. Temperatures from photometry

Figure 8 shows the histogram of fitted dust temperatures T_dp, subject to f_d> 0.1 to exclude objects where a fit of dust properties is not meaningful. The histogram of T_dp shows a few possible outliers above 2200 K or below 1200 K; to avoid any excessive influence of such outliers, we make use of the median statistic. The median value of T_dp is ≃1650 K (log ₁₀(T_dp) = 3.22). The median of the error estimates on log ₁₀(T_dp) is: 0.026, consistent with an estimate of the sample dispersion as: median[| log ₁₀(T_dp)−median(log ₁₀(T_dp)) |] = 0.030, and with the dispersion being mostly from measurement uncertainties.

The median value for T_dp appears to be in the range of dust sublimation temperatures, in agreement with the current ideas regarding disks around HAeBe stars; see Dullemond & Monnier (2010) and references therein, or the disk-wind model of Bans & Königl (2012). However, in any such model, one expects dust to exist over a range of temperatures below the sublimation temperature; furthermore, our fitting procedure might introduce some bias. To address these two concerns, we use the MCFOST radiative transfer code (Pinte et al. 2006) to generate the SED from a few simple models of HAeBe objects combining two spectral types (for the central star) with two types of grains. The inner rim radius of each model is adjusted so that the maximum grain temperature would be (approximately) the sublimation temperature T_sub for the respective grain types: 1500 K for silicates, and 1800 K for carbon. These models are used only to discuss dust temperatures, and are distinct from those defined in Sect. 4.2. The values of νF_ν in the eight Johnson bands (U−K) are used as input to the same fitting procedure that is used for deriving T_dp from the observations. A sample fit is shown in Fig. 9. The results are summarized in Table 9; the four values of T_d,fit are plotted in Fig. 8; the models with carbon grains appear to agree better with the observed values of T_dp; in fact, the difference in log ₁₀ (T_dp) between carbon and silicate models, 0.0765, is 2.7 times larger than the (median-based) sample dispersion. We conclude that: (1) the best-fit blackbody temperature for the NIR emission is close to 1650 K; (2) based on simple models, that value (1650 K) is consistent with T_sub = 1800K (carbon?); (3) the data rule out that the grains responsible for the NIR bump have T_sub = 1500K (silicates?).

Our finding that the sublimation-bounded inner rim of Herbig object disk is populated predominantly by grains with a high sublimation temperature can be explained by the survival of the fittest in a hostile environment. A similar conclusion has been reached by Carmona et al. (2014) in their study of HD 135344B.

4.1.2. Temperatures from interferometric data

When spectrally dispersed interferometric data has been acquired, the model fit to the interferometric data provides an estimate of the slope, across the H band, of the radiation intensity from the (partly) resolved component. More precisely, the fit provides the differential slope between the extended and the unresolved (stellar) components; the latter being known from the star’s spectral type.

Define for any component the spectral index (across the H band) as: (11)that can be converted easily to and from blackbody temperatures. The spectral slope of the circumstellar component, determined from interferometric data, k_c, can be converted to an equivalent temperature T_ci.

In Fig. 10 we present a scatter diagram of T_ci versus T_dp. The data are from the Ellipsoid fits to the 44 non-binary program objects, with extra restrictions: (i) dust flux fraction (H-band photometry) f_d ≥ 0.1; (ii) circumstellar flux fraction (H-band interferometry) f_c ≥ 0.1; (iii) object HD 56895 excluded (see Sect. 2.2 above); (iv) object must be observed in dispersed mode. A total of 35 objects meet these conditions.

First, the ranges of T_ci and T_dp are similar and not far from expected grain sublimation temperatures. Second, there is no clear correlation between these two variables. Third, the values of T_ci are predominantly lower than T_dp for the same object. We remind the reader that T_dp is derived by decomposing the light from a Herbig object in the spectral domain, while T_ci relies on a decomposition in the spatial domain, where the circumstellar component is assumed to be dominated by radiation from dust. The identity between the respective stellar and non-stellar parts in each of these decompositions is, however, an assumption.

Following our caveat above on the respective definitions of T_ci and T_dp, one possible cause for the observed discrepancy would be a non-photospheric and unresolved component of the object. For instance the inner gaseous disk that should be present inside the dust sublimation radius, at least if the star is still accreting. The determination of T_ci relies on the differential slope in the H band of the (partly) resolved component relative to the unresolved component. If a non-photospheric emission is present with a flux rising (towards short wavelengths) faster than the photospheric flux, then assigning the photospheric value to the slope of the unresolved component is an under-estimate, which propagates to an under-estimate of the dust temperature. The existing models for such an inner gaseous disk, for example, Muzerolle et al. (2004; their Fig. 11), or Dullemond & Monnier (2010; their Fig. 13) based on the opacities of Helling & Lucas (2009) show that at least potentially, its SED can upset the slope of the photospheric emission.

4.2. Overcoming some limitations of parametric fits via substitution weighing

In the previous section, parametric fits of the observed data were used to obtain quantitative information in preparation for the discussion and statistical analysis of the present section. The values of fitted parameters can be subject to biases: either identifiable like the difference between the radius of a Ring model (mean radius) and the rim radius (innermost), or, worse, unforeseen and undetected.

To overcome as much as possible these issues, we use an approach analogous to substitution weighing (that allows accurate measurements using a balance with unbalanced or even unequal arms). We generate a small number of physical disk models, using the MCFOST code.

Each model is then “observed” at various inclinations and distances. Here, to observe means to generate visibilities squared and closure phases for a set of wavelengths and baselines typical of a program object (HD 144668). A realistic noise is added, based on the estimates derived from extensive analysis of PIONIER data; for details see Appendix A. These data are fed to the same fitting procedure as the real data. The comparison between the fit parameters from, respectively, real and model objects, helps to interpret the fit results in a way that eliminates possible systematic biases from the fitting procedure.

Our models are loosely inspired by two published models: Isella & Natta (2005, hereafter IN05), and Tannirkulam et al. (2007, hereafter THM07); we take full responsibility for the shortcomings of our imitations; our goal is mainly to sample two different morphologies. Each of IN05 and THM07 models derives a consistent solution for radiative transfer, thermal equilibrium, hydrostatic equilibrium, and dust sublimation. Furthermore, THM07 considers two grain sizes and size-dependent settling, as considered in detail by Dullemond & Dominik (2004). The smaller grains, having a lower cooling efficiency, reach the sublimation temperature at larger distances than the larger grains.

Our IN05-like model, IN hereinafter, has a single grain radius: 1.2 μm while our THM07-like model, THM hereinafter, has two grain radii: 0.35 μm and 1.2 μm, resulting in a broader radial distribution of H-band emission. In contrast with the IN05 and THM07 originals, our models solve only for thermal equilibrium of dust and radiative transfer, given an arbitrary dust distribution with azimuthal symmetry. See Fig. 11.

Variants of each of the base models were produced for three types of central stars: A8V, A0V, and B2V. Spatial scales ℓ for each model were adjusted ( to first order) to bring the maximum grain temperature close to 1800 K, while densities were adjusted ρ ∝ ℓ^-1 to preserve opacities at analogous points. Anticipating the observed value of the fraction f_rep of stellar flux intercepted by the dusty disk (Sect. 4.7.2), we have tuned the relative thickness z/r of our models to match f_rep ≈ 0.2 for the A8 and A0 variants, while we targeted f_rep ≈ 0.05 for the B2 variants. More details are provided in Appendix A.

Additionally, to test the prevalence of false positives (detection of ring structure) in our fitting procedure, we generated the ELL series of images, that are ellipsoids, having the same values for the stellar flux, the circumstellar flux, and the half-light radius as the THM models; in other words, the “same”, but no ring structure. These images are generated from a prescription for their light profile, not from a radiative transfer calculation, and should not be confused with the Ellipsoid fit model introduced in 3.5.

4.3. Is there evidence for ring structure?

One possible way to answer that question is to compare the values for an Ellipsoid fit and a Ring fit in its simplest form, with no azimuthal modulation. Figure 12 shows the results, over the HQ sub-sample, sorted by size. One can see that: firstly, six objects express a preference for the ring structure: VV Ser, HD 190073, HD 98922, MWC 297, HD 100453, HD 45677; secondly, a large angular size favors, but does not guarantee that the ring structure results in a lower . A natural (but not unique) interpretation of these results is that (i) a significant fraction of our program objects have a ring structure; (ii) the detection of such a structure is difficult or impossible in marginally resolved objects.

The same question – whether our data show evidence for ring structures in the program objects – can be addressed in a different way. Remember that the Ring model fits comprise the Ellipsoid fits as a subset, since as log ₁₀(a_k/a_r) → ∞, w → 1, and the ring structure degenerates into an ellipsoid; see Tables 7 and 8. So, instead of examining the of a ring fit, we examine the value of w, and, since we suspect that the ability to detect ring structures is hampered by the limited angular resolution, we plot w versus a; see Fig. 13 (left panel). The following results emerge: (1) a negative correlation between w and a, i.e. the detection of a ring structure is less likely for small radii, as expected; (2) for a significant fraction of the objects, a ring structure is detected by the fit, but... (3) the rings seem to be wider than expected models with a rounded rim, as proposed by Isella & Natta (2005) and pictured in a number of publications. At this point, we must consider the following possibilities: (a) our parametric fit is biased, the fitted ring is wider than reality; (b) actual rings are wide (w of order 0.5); (c) there are no rings in the H-band images of HAeBe objects circumstellar matter, and the fit algorithm was biased towards finding them.

Figure 13 (right panel) shows, alongside the HQ objects, the fit results for the IN models at various inclinations and distances. Although we cannot compare one-to-one the fitted ring width with the real width (if only because of the lack of a unique definition for the latter), this plot makes it clear that the putative rings of the LP objects are wider than those of the IN model.

Also in Fig. 13, we have plotted the results of fitting THM models. This shows that at least one simple model of dust distribution and radiative transfer can result in a (a, w) diagram, if not identical, at least much closer to the observed one than the IN model.

We remind the reader that the IN model, inspired by Isella & Natta (2005), has a rim with a rounded shape, with a (relatively) large radius of curvature, and a single grain size, while the THM model, inspired by Tannirkulam et al. (2007), has a shape closer to a knife-edge and two grain sizes; its sloping surface under oblique illumination leads to a more radially extended emission.

Finally, performing fits of plain ellipsoid images ELL (defined in Sect. 4.2) we find only a small rate of false ring detection: 90% of the w values are >0.925, and 80% are >0.975; while w < 0.7 occurs only for small angular radii a < 1 mas.

We conclude that, using the substitution weighing approach, we have eliminated possibilities (1) and (3) as listed above, leaving only (2): our interferometric data support the presence, in a significant fraction (at least 50%) of the HQ program objects, of a ring-like structure. Moreover, these ring structures are significantly wider than what results from a simple rounded rim as our IN model, and a fortiori for a flat vertical wall.

The THM model, however, fits the observed (a, w) diagram rather well. Some doubt concerning a sharp inner rim had already been voiced by Tannirkulam et al. (2008b); these authors interpreted the lack of a strong bounce in the visibility curve as supporting the presence of emission inside the sublimation rim. It is noteworthy that, besides THM07, Kama et al. (2009) and Flock et al. (2016) also perform a self-consistent modeling of the rim including grain sublimation, and that most of their models result in a wedge-shaped (as opposed to rounded) rim. The photospheric contours in their Figs. 6 and 7 (bottom rows) bear a striking similarity with those of our THM model (Fig. A.2, bottom row, middle panel), considering that our THM model (not physically self-consistent) was adjusted specifically to match the large value of the observed width of the H-band emission region. Dullemond & Monnier (2010, their Fig. 10 and Eq. (13)) also explicitly derive a physical model with a large radial extension of the sublimation front.

4.4. Do we see azimuthal modulation?

To answer this question, we have performed, over the HQ sub-sample, three types of model fits, beginning with the Ring model with no azimuthal modulation; see Table 10.

Because we keep a record of the partial χ² contributions, we are able a posteriori to choose to include the closure phase (φ_cp) contribution in the final value of χ² for the m = 0 model, allowing a fair comparison with the other two models.

The values of the three models listed in Table 10 are compared in Fig. 14, with the un-modulated model serving as a reference. As in Fig. 12, the objects are ordered by increasing angular size (value of a). One can see that most, but not all, of the better-resolved objects show an improvement of χ² when modulation is introduced; m = 2 modulation seems to provide a smaller incremental improvement than m = 1.

A deviation from centrosymmetry (odd m modulation) is expected when a rim is viewed at non-zero inclination. But a spatial offset between the circumstellar emission and the star has the same effect. So, we propose a slightly more stringent test in the next subsection.

4.5. Is the azimuthal modulation aligned as we expect?

In any model featuring an inner disk rim heated by the central star, one expects, at intermediate inclinations, the far side of the rim to be more luminous than the near side. Such an asymmetry involves, in Eqs. (8) and (9), only sine terms of odd order; other terms are included because the model must be unbiased.

We test the presence of such a far-side dominance using a plot based on the modulation coefficients c₁ and s₁. We remind the reader that the modulation is defined in the rotated frame that is aligned with the apparent shape of the rim, and note that the position angle of that ellipse is (nearly) indeterminate for (near) zero inclination (pole-on). Accordingly, we represent each object at the coordinates (x = | c₁ | sini,y = | s₁ | sini). See Fig. 15. The asymmetry between the near and far sides of a rim image is expected to result in points aligned along the vertical axis. As a reference, the right panel shows the result of fits to in_A0 and thm_A0 models. The scatter diagram on the left panel of Fig. 15 show a majority of points closer to the vertical than to the horizontal axis, similar to, but less marked than the model diagram on the right panel.

One notable exception is MWC 158 (at graph coordinates (0.85,0)); this object has been observed with PIONIER before and during the Herbig Large Program, over a time span of ≈2 years, and has been shown (Kluska et al. 2016) to have a variable structure on ≈2mas scale, possibly corresponding to the orbital motion of a bright condensation around the central star.

Following this special examination of MWC 158, which does not fit into the expected pattern for modulation coefficients, we must examine with equal attention the objects that do conform to the expected pattern. Twoe objects dominate the alignment along the vertical axis in Fig. 15: HD 45677 and HD 144668. The values of w for these two objects are:

Figure 13 (test for false positives with ELL images) shows that a w value below 0.95 corresponds to a ring structure with a high degree of confidence. Therefore, the visual alignment pattern in Fig. 15 does arise from modulated rings, not modulated ellipsoids (in which case it would be essentially an offset).

4.6. Closure phase statistics

Closure phases provide a model-independent way to measure the degree of asymmetry of the objects. We generate plots of closure phase φ_cp versus resolution ρ′, as follows: (1) we measure the resolution by the largest of the three projected baselines of the triangle that defines each closure phase measurement; this is an arbitrary way to subsume a triangle by a scalar; (2) the resulting scatter plot is converted to an ordinary plot of closure phase versus resolution by computing percentiles in a sliding resolution window (width 0.2 in log ₁₀ representation); (3) the effective resolution must take into account the angular radius of the object, and, to first order, its anisotropy. Accordingly, for an object having half-light semi-major angular diameter a at position angle θ, with an axis ratio cosi, and observed over a projected baseline with coordinates (u,v) = (b_x,b_y) /λ, we define the dimensionless resolution σ′ along that baseline: (12)and for a closure phase measured over a triangle, the value of σ′ for the largest side: (13) Figure 16 shows (black lines) plots of the sliding third quartile (see above) of φ_cp versus ρ′ either for the totality of the 44 non-binary objects, or excluding HD 45677. Monnier et al. (2006) used similar plots, but based, at each resolution value, on the maximum closure phase, which may have a low probablity of being measured in actual observations. We show two variants, with the measure of dimensionless resolution based either on the largest or the smallest side of each closure triangle.

The plot that includes HD 45677 shows higher values of φ_cp above ρ′ ≈ 1. HD 45677 is the object with the second largest angular radius (≈10mas) in our sample, almost four times larger than the next largest object. Image reconstruction (Kluska et al. 2014) and parametric models up to azimuthal order m = 5 agree to show substructure in the sublimation rim of HD 45677; we believe that the closure phases contributed by HD 45677 data above ρ′ ≈ 1 belong to such substructure rather than to the ring asymmetry. Finally, the object with the largest angular radius, HD 179218 (≈25mas) does now show a clear closure signal (see Fig. D.2).

Closure phase peaks around ρ′ ≈ 0.5, as already noted by Monnier et al. (2006), and the peak value for our data set is 13° or 7.5° respectively for the two versions of the graph. Aiming to relate these values with properties of the sublimation rim, we also show in Fig. 16 similar plots, for the MCFOST models, IN and THM, already used to calibrate our parametric fits. Borrowing from the results of Sect. 4.7.1, we have applied to these simulated observations the same cosi> 0.45 cutoff as found for actual data. One can see that none of the six model variants matches the observed φ_cp versus ρ′ plot. Possibly in_A0 comes closest, but the IN model has been discarded, on the basis of the ring width, in favor of the THM model. On the other hand, the THM model fails to reproduce the observed closure phases. One can also note that the value of ρ′ where φ_cp peaks is not uniform, and can vary between ≈0.3 and ≈0.7. Clearly the closure-resolution diagram is a promising sieve to reject or validate models. But, in keeping with our disclaimer in 4.2, we do not attempt, in this paper, to tune a model to fit all our observables.

4.7. Disk thickness

We attempt to gain information on the disk thickness (z/r) via two approaches: the statistics of inclinations and of reprocessed flux fraction.

4.7.1. Inclination statistics

Figure 17 shows the histogram of cosi fit values for either all non-binary objects, or the HQ subsample. A flat histogram over the full (0,1) range would be expected for a population of thin rings with random orientations. Apart from the small-number fluctuations between the two fit models, we note that: a) there is no significant pileup around cosi = 1, as would result from a contamination by a population of intrinsically spherical objects; b) there is a lower cutoff around cosi ≈ 0.45. We interpret this cutoff as due to self-obscuration, i.e. objects do exist at higher inclination, of course, but are absent from the sample because either the central star or the H-band dust emission is obscured. Let η = cosi_cutoff, then the relative thickness at which obscuration sets in is given by: ζ = z/r = η/ and ζ(η = 0.45) = 0.50. We find no significant correlation of the fitted cosi with the spectral type.

Fitting the THM models reveals a small (<0.02) positive bias in the values of cosi, which does not affect our interpretation.

The observed cutoff may be a consequence of flaring and obscuration, by the outer disk of either the NIR radiation of the inner disk or the central star itself. We do not attempt to mimic this situation with our MCFOST models because our models extend radially only far enough (typically 10 × rim radius) to model the NIR emission up to K band. If indeed flaring is responsible for the observed cutoff, we can simply state (without the need for numerical modeling) that the radial direction having an extinction A_V ≈ a few has a slope (with respect to the equatorial plane) of order z/r ≈ 0.50.

4.7.2. Reprocessing of stellar flux

Another approach to the disk relative thickness is to consider the fraction of the stellar luminosity reprocessed to near-infrared. More specifically, we use the fits of the SED as discussed in Sect. 3.1, to derive: a) the de-reddened stellar bolometric flux F_{bol ⋆}; b) the flux of the blackbody component F_NIR. Remember that the data used in the fit do not extend longward of K band, so only hot dust is considered. We plot the histogram of the reprocessed energy fraction f_rep = F_NIR/F_{bol ⋆}, for the set of non-binary objects; the result is shown in Fig. 18. One can see that the histogram has a cutoff around 0.5. One object, CO Ori, is outside the limits of the histogram, at f_rep = 1.64; that object has a significant extinction A_V ≈ 1.2 and a spectral type G5, atypical for the HAeBe class.

The value of f_rep would have a straightforward geometrical interpretation in the case of a disk rim fully opaque to incident stellar radiation, but optically thin for re-emitted dust radiation; in that case, its value is simply the solid angle, as seen from the star, intercepted by the rim. Because of (a) the self-obscuration of the disk in the NIR; and (b) some dependence of the emergent radiation on the inclination relative to the local normal direction, one can expect the reprocessed fraction f_rep to depend on the inclination.

We can, however, compare the median value of the observed sample with the similar metric for a model disk to derive a typical value. The median values of f_rep are: for the 44 LP non-binary objects: 0.18; for the in_A0 MCFOST models: 0.19; for the thm_A0 models: 0.18. The median for MCFOST models is restricted to inclinations for which the extinction of the central star A_V < 3. Remember (Sect. 4.2) that the thickness of the IN and THM models was adjusted to achieve such an agreement. The scatter of f_rep seen in Fig. 18 probably reflects the diversity of the objects included in our sample.

There is, in principle, a flaw in this procedure. In the previous subsection we argued that the flaring of the outer disk restricts the observed HAeBe disks to the inclination range 1 ≥ cosi ≥ 0.45, while we are now calibrating the observed statistics of f_rep against models without any outer flaring, that remain observable over a wider range of inclinations, typically 1 ≥ cosi ≥ 0.2. However, the model with the strongest dependence of f_rep versus inclination (a vertical-wall rim) is known, even before the present work, to be incompatible with the visibility curve. The favored THM model has only a weak dependence: Δ f_rep/ Δ cosi ≈ 0.1. This means that selection effects on cosi do not affect significantly our determination of the thickness z/r of the inner disk.

Another question is: why should the geometrical interpretation of the reprocessed fraction be correct? After all, the albedo of the grains in our models is close to 0.5 over the range of stellar photon wavelengths. But the diffusion is strongly forward-peaked; in MCFOST, with the Mie treatment, for a = 1.2 μm grains, and at 450 nm, the phase function drops to half the peak value within ≈5 ° of the forward direction; so, diffusion will not hinder absorption. The situation is of course different for 1.65 μm stellar photons diffused off grains (typically ≈0.1μm) in the outer parts of the disk’s surface.

Summarizing the results of the previous and present section, we find typical values z/r ≈ 0.2 for the inner disk, and z/r ≈ 0.45 for the outer disk are needed to match the statistics of f_rep and cosi. The value z/r ≈ 0.2 is known from previous work (see e.g., Vinković et al. 2006, Mulders & Dominik 2012 and references therein) to be larger by a factor of approximately two than predicted by hydrostatic models. See discussion in Vinković (2014). One possibility is that magnetic support increases the disk thickness over the hydrostatic value, as proposed by Turner et al. (2014).

4.8. Luminosity-radius relation

The disks of Herbig AeBe stars are thought to have an inner radius bounded by dust sublimation. This can be tested quantitatively via the luminosity-radius relation (Monnier & Millan-Gabet 2002; Monnier et al. 2005). We present in Fig. 19 the radius-luminosity diagram for the 27 objects of the HQ sub-sample. It is prepared as follows: (a) our SED fit provides a de-reddened m_V; (b) using Simbad spectral types (Table 1) and bolometric corrections from Pecaut & Mamajek (2013), we derive bolometric apparent magnitudes m_b; (c) using parallaxes from Hipparcos (Perryman et al. 1997) we convert m_b to absolute luminosity L_b and angular to linear radius. Some of our HQ central stars have benefited from detailed spectral studies: 15 of them by Fairlamb et al. (2015), and one by Montesinos et al. (2009); see Table 1. For these objects the value of L_b derived as outlined above and the distance were superseded by the (presumably more reliable) values from these two references. Note that an error on the distance of an object moves the corresponding point in the diagram parallel to the theoretical log L−log a line; so, accurate distances are not critical to establish a possible correlation.

On the same plot, we have drawn the (or rather a) theoretical luminosity-radius relation. Such a relation is defined by Tuthill et al. (2001); Monnier & Millan-Gabet (2002) mention and neglect backwarming, i.e. treat the equilibrium temperature of an isolated dust grain, while Kama et al. (2009) explicitly consider backwarming in the range from C_bw = 1 to C_bw = 4, the latter value being implicitly used by Dullemond et al. (2001). The rim radius can be written as: (14)where ϵ = Q_abs(T_sub) /Q_abs(T_∗), the cooling efficiency, is the ratio of Planck-averaged absorption cross-section for the considered grain species.

The line drawn in Fig. 19 expresses Eq. (14) with T_sub = 1800K, ϵ = 1 and C_bw = 1. The sublimation temperature is fixed (see Sect. 4.1.1 above). The values assigned to the other two are, however, arbitrary, and the relation as represented is best called a reference, rather than theoretical, relation.

Several factors can cause the actual objects to deviate from that reference: (1) ϵ is expected to be ≤ 1, reaching unity for large (a_grain>λ_NIR/ (2π)) grains. Tuthill et al. (2001) argue that the largest grains (largest ϵ) are the last survivors, but one needs to assume that grain growth indeed proceeds to ≈0.5μm size; (2) the value adopted for the reference line, C_bw = 1 is the smallest possible; (3) the observed radius, with limited spatial resolution, is not the inner rim radius, but an effective radius for the H-band emitting region, that will in all cases be larger than the inner rim radius.

All three of the listed factors tend to make the measured radius larger than the reference radius (at a given luminosity). To investigate this quantitatively, we display in Fig. 20 the histogram of the ratio r_norm = a_fit/a_sub,1800 of these two radii. The median value is median(r_norm) = 2.0. Note that the sublimation temperature of the grain material is not used as an adjustment variable.

Still using our procedure of substitution weighting, we process, with the same fitting procedure, images from three types of MCFOST models: (1) the IN model achieves values of r_norm in the range 1–2, mostly due to backwarming; (2) the THM model achieves larger values of r_norm than the IN model, due to the larger radial extent of the H-band emitting region; it does not, however, seem to match the observed values; (3) a third model, THMC, is similar to THM, but with smaller grain sizes and correspondingly smaller ϵ values: 0.147 μm (ϵ = 0.14) and 0.235 μm (ϵ = 0.28); the range of r_norm values is correspondingly larger.

Taken together, the three models IN, THM, and THMC can account for a spread of r_norm over the range 1–4. HD 100453 (r_norm = 4.5) exceeds that range by a small amount, while the five objects with r_norm < 1 that can be identified in Fig. 19 are more problematic. The objects with the four lowest values of r_norm, ranging between 0.18 and 0.59, all have dust temperatures (see Sect. 3.1) below 1760 K, so that anomalous refractory grains cannot be invoked. Anticipating the results displayed in Fig. 23, we also see that an H-band emission component located inside the sublimation rim cannot plausibly explain a fit radius smaller than the reference radius.

4.9. Extended component

Our data show the presence of a fully resolved flux in a significant fraction of our sample. This is quantified by the f_h fitted parameter, which is a fraction of the total H-band flux in the instrument’s field of view. This can be attributed to a spatially extended emission⁴; see for example, Fukagawa et al. (2010, their Table 5).

Based on ISO spectra Meeus et al. (2001) defined two types of Herbig Ae Be stars: group I are stars for which the mid-infrared continuum spectrum can be modeled with one power law and one additional black-body to reproduce the rising mid-IR continuum; group II objects only require the power law. In subsequent mid-IR spectroscopy studies van Boekel et al. (2005) and Juhász et al. (2010) noticed that both groups could be sorted apart in a diagram involving the near-infrared (J, H, K, L and M) to far infrared flux ratio (IRAS 12, 25 and 60 μm) and the IRAS color m₁₂–m₆₀⁵. These parameters have been considered as indicators of respectively the ratio between the flux intercepted by the inner rim and the outer disk and the degree of flaring.

We explored whether we could find a correlation between such indicators and our own estimator of extended flux (see Fig. 21). The latter was constructed as the ratio between the extended flux fraction f_h and the global excess fraction 1−f_s. We do find a clear indication that objects with a higher ratio of near infrared to infrared flux display small fully resolved flux while higher mid infrared color indices are correlated with higher extended flux. This fits well in the scenario where indeed group I and group II are objects with significantly different large scale near to far infrared emission.

We remind the reader of the discussions (Maaskant et al. 2013; Menu et al. 2015) concerning the exact nature of the group I/II classification and the attempts to introduce notions of disk evolution (e.g., dust settling) and planet formation (e.g., gap clearing). It has even been proposed that the near-infrared fully resolved flux might be caused partly by scattering off the inner rim of an outer disk in a transitional disk containing a large gap (Benisty et al. 2010b; Tatulli et al. 2011). We note that the f_h parameter provides a limited spatial information and do not wish to speculate more on its exact origin. In particular we recall that single-mode observations have a very restricted field of view, for PIONIER a Gaussian of roughly 200 mas field of view, which introduces a bias in the gap detection capability. But, nevertheless, it would certainly be interesting to correlate it with other planet formation evidence.

4.10. Inner gaseous disks

The possibility of an inner gaseous disk located inside the dust sublimation rim of Herbig objects has drawn some attention in recent years, from several different viewpoints. Firstly, as logically unavoidable if the dusty disk must feed accretion at or near the stellar surface (Muzerolle et al. 2004). Second, a number of authors have directly detected gas inside the dust sublimation radius of various Herbig objects via spectrally resolved interferometry in the Brγ line, for example: Kraus et al. (2008a), Eisner et al. (2009), Caratti o Garatti et al. (2015), Kurosawa et al. (2016), Mendigutía et al. (2015), Ilee et al. (2014). Third, continuum (or low spectral resolution) observations, for example, Eisner et al. (2007), Isella et al. (2008), Kraus et al. (2008b), Tannirkulam et al. (2008b,a), Benisty et al. (2010a), infer the presence of an additional component inside the dust sublimation radius, that can be uniform, ring-like, or compact; note that the authors of the last paper argue in favor of an inner refractory dust disk. Fourth, some authors have investigated the opacities and/or the emitted spectrum from a gaseous disk as might be found inside the dust sublimation radius; see Malygin et al. (2014), Vieira et al. (2015), with the former authors concluding that the thermal equilibrium problem may be multi-valued.

In view of the diversity of spatial distributions that have been proposed for the inner disk continuum emission, and that cannot at present be tied to a physical model, we could not submit the hypothesis of dust rim plus inner gas disk to the same systematic tests as we did in the previous subsections for the dust-rim-only hypothesis. The issue is not the presence of gas inside the dust sublimation radius, for which the Brγ observations provide direct evidence. But, line opacities being typically several orders of magnitude larger than continuum ones, this does not prove that an inner disk should be a significant contributor to the continuum image and visibilities.

Nevertheless, we did a simple test of a model comprising:

• a ring similar to the Ring model used above, except that an elasticconstraint enforces l_kr = −0.3 ± 0.15, corresponding to a relative width , intervals corresponding to Δ(χ²) = 1. Such rims are narrower than any of the fit results in Fig. 13;

• a disk with a structure like the Ellipsoid model used above, and again with an elastic constraint log ₁₀(a_disk/a_ring) = 0 ± 0.1, both a values being half-light radii, and the interval given, again, for Δ(χ²) = 1.

We show just two results of this experiment. Figure 22 shows the improvement in χ² from an Ellipsoid fit to a ring+disk model. Comparing with Fig. 12, one can see that the two are quite similar. Figure 23 compares the values obtained for the ring radius in the present versus the Ring model fit; again, the values are quite similar: the ratio has a mean of 0.975 and a standard deviation of 0.107.

What can we gather from this limited test of the ring+disk model? The results are quite close to those of the Ring model. With our data, we have typically just enough resolution and precision to distinguish between a sharp rim and any other structure with a more damped visibility versus resolution curve, whether that is a wide rim, or a narrow one combined with an inner disk, central clump, or smaller ring.

5. Summary and conclusions

We have presented interferometric H-band observations of 51 Herbig AeBe objects, of which 44 have a non-binary star, and 27 have been classified high quality (HQ) on the basis of the resolution achieved; a dedicated and homogeneous set of photometric measurements was also acquired.

The spatial resolution achieved for a significant part of our sample allowed us to go beyond just size determinations, but, for most objects, was not sufficient to produce images. In order to make the fullest use of our data set, we examined the statistics of structural parameters obtained by fitting simple geometrical templates: ellipsoids and rings. We separately generated a small number of physical models of dust disks irradiated by the central star; they were observed with a uv coverage similar to actual sample object, and these synthetic observations were in turn fitted exactly like the real observations. This allowed us to (i) eliminate possible biases of out fitting process; (ii) make statements on physical models, and not just on structural parameters.

We have established, or confirmed, a number of properties of the dusty environment of HAeBe objects, that we summarize below under headings corresponding to the four goals spelt out in the Introduction.

5.1. Can we constrain the radial and vertical structure of the inner disk?

Radial structure.

We have confirmed the presence of a ring-shaped structure (that tends to be taken for granted) in the Herbig objects of our sample, in several ways. First, for a substantial number of objects, fitting a ring structure results in a lower reduced χ² than fitting an ellipsoid, this being more marked for the better resolved objects. Second, performing a parametric fit that allows both an ellipsoid or a ring shape, with a continuous transition, for a number of objects, a ring shape is preferred. Also, when a ring shape is preferred, it is a wide ring, with w > 0.4, where w is the ring width normalized by the global half-light radius. This excludes “fat” rim geometries with a substantial curvature radius and favors wedge-shaped rims, concurring with several recent physical models.

Luminosity-radius relation.

We confirm the known luminosity-radius relation. The observed spread can be explained in part by variations in the rim geometry and in the maximum grain size, but a few objects have radii that are smaller than expected. The observed spread in sublimation temperatures can only be a minor contributor to the luminosity-radius scatter.

Vertical structure.

The statistics of the reprocessed fraction (from the star’s bolometric flux to the dust’s NIR flux) allow us to infer a thickness z/r ≈ 0.2 for the inner (meaning close and hot enough to radiate in the NIR) region of the dusty disk. The statistics of inferred inclination show a histogram populated only at cosi > 0.45, that we interpret as resulting from self-obscuration by the outer parts of a flared (z/r ≈ 0.5) disk.

5.2. Can we constrain the rim dust temperature?

We find that the dust responsible for the NIR emission of our sample objects has a sublimation temperature definitely larger than 1500 K, the favored value being 1800 K. This result is obtained by using radiative transfer models to calibrate out both the spread of temperatures in any real object, and possible biases in our fit procedure. Relating a sublimation temperature to dust properties is a difficult task, considering the range of the physical conditions and possible time-dependent effects. Accordingly, we abstain from discussing dust composition.

5.3. Can we constrain the presence of additional emitting structures?

Extended component.

We find in a significant fraction of our sample a spatial component that is fully resolved at the shortest spacings (B/λ ≈ 5 × 10⁶), i.e. an angular scale larger than 40 mas. We also find that the presence and relative flux of that component are correlated with the Herbig group, either the original categorical form or based on infrared colors. Diffusion of starlight off a flared disk might be the cause of such extended components, but there may be other viable mechanisms.

Inner disk.

Our conclusion (above) that rims are wide results essentially from the (near-)absence of bounce in the visibility curve beyond the first minimum. A number of authors have interpreted similar observations as resulting from an extra component (hot gas, very refractory grains?) inside the dust sublimation radius. Such an alternate explanation, while viable, is loosely constrained, lacking a predictive physical model.

5.4. Can we reveal deviation from axisymmetry?

Azimuthal modulation.

Fit models that include azimuthal modulation result in lower values of ; such azimuthal modulation is expected from self-shadowing of an inclined rim. To probe further that agreement between observations and models, we also tested the alignment of the modulation; we find a suggestive, but not overwhelming, agreement with the expected alignment.

Phase closures.

This item is closely related to the previous one, except (a) model-independent, and (b) only odd azimuthal orders contribute. Phase closures versus resolution show, as expected, a maximum near dimensionless resolution ρ′ ≈ 0.5. The value of the maximum is of order φ_cp ≈ 10°; one of our physical models (IN) comes close to this value, but the other one (THM) peaks much lower, at φ_cp ≈ 2°; we leave this point open for future work.

5.5. Outlook

Our understanding of the structure and physics of Herbig AeBe objects would benefit, in our opinion, from: (a) an increase in angular resolution of NIR interferometric observations by a modest factor (between two and three), that would provide a decisive improvement for all the objects that are marginally resolved with current data; (b) the continuation and extension of theoretical or numerical models of the inner disk regions, notably with the following issues in mind: the shape of the rim; its vertical structure and hydrostatic equilibrium; a comprehensive treatment of the dusty rim and the inner gaseous disk; the deviations from axisymmetry.

---

Acknowledgments

Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 190.C-0963. Generous time allocations by SAAO are gratefully acknowledged. We are grateful to Francois Van Wyk for diligently carrying out the infrared part of the photometry at Sutherland. J.K. acknowledges support from a Marie Sklodowska-Curie CIG grant (Ref. 618910, PI: Stefan Kraus). J.D.M. and F.B. acknowledge support from NSF-AST 1210972. Part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This research has made use of: the Jean-Marie Mariotti Center Aspro2 and SearchCal services⁶; the SIMBAD database, operated at CDS, Strasbourg, France; and data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia).

References

Appendix A: MCFOST models used for double-weighing the parametric fits

Here we describe in more details the MCFOST models used in conjunction with the parametric fits. We populate the inner rim region with grains having the highest sublimation temperature, meaning carbon grains (see Sect. 4.1.1), and the lowest equilibrium temperature, meaning large grains; altogether, the grains at the inner rim are the fittest for survival. We use in the MCFOST program amorphous carbon grains, with the optical constants from Rouleau & Martin (1991). Figure A.1 shows, versus grain radius a, the cooling efficiency ϵ, the ratio of the Planck-weighted averages of the absorption cross-section at, respectively, the dust and star temperatures.

Our IN05-inspired model, labeled IN, has a curved inner rim and a single grain radius: 1.2 μm. Our THM07-inspired model, labeled THM, has two grain radii: 1.2 μm and 0.35 μm; the smaller grains, because of their smaller emissivity in the NIR, cannot survive quite as close to the star as the larger grains, but, where they survive, their equilibrium temperature is higher than for larger grains, which helps to achieve a radially extended H-band emission. As explained in Tannirkulam et al. (2007), the respective spatial distributions of the two grain species result naturally from differential sedimentation and differential sublimation. In the present work, however, the spatial distributions are ad hoc, subject to the condition that the maximum temperature for each species is equal to the assumed sublimation temperature T_sub = 1800 K. Both the geometry and the range of grain sizes contribute to the THM model having a broader radial distribution of H-band emission than the IN model. The information provided in Fig. A.2 complements that already given in Fig. 11.

We simulate observations of each MCFOST model as follows: (1) we vary the inclination with values of cos(i) at the centers of 10 uniform bins; (2) we vary the distance in steps of 10^0.1 subject (for each inclination) to the magnitude limits of our HQ sub-sample: 4 ≤ H ≤ 8; note that we claim neither that our sample is complete within these limits, nor that the distribution of distances in our simulated sample matches the observed sample, but we believe that this serves sufficiently well our purpose: to reveal trends while avoiding possible systematic biases caused by parametric fitting; (3) for each combination of inclination and distance, the uv plane is sampled following the coverage of a typical Large Program object (HD 144668); noise is added; values of visibility squared V² and phase closure φ_cp are derived; (4) inclinations for which the V band extinction of the central star is in excess of three magnitudes are ignored; while the limit is somewhat arbitrary, heavily obscured central stars will not appear in lists of HAeBe objects.

Simulated noise is added to the model’s complex visibilities; the resulting noise for the (real) observables V² and φ_cp has median values slightly larger than the median values for the error estimates derived by the data reduction software for the real LP objects. We feel that matching the median error between real and simulated data is adequate for our purposes, and that a discussion of noise models for interferometric data is outside the scope of this work. We do not attempt to simulate the calibration errors or other systematic errors.

Finally, the simulated observations of the MCFOST models are subject to the same parametric model fitting as the observational data.

Appendix B: Results of photometric and interferometric fits

Appendix C: Log of observations

Appendix D: Summary plots for program objects

All Tables

All Figures