© ESO 2020

1 Introduction

The study of protoplanetary disks around young stars is fundamental to understanding planet formation: the initial conditions for planet formation are indeed determined by the disk properties, its dust and gas densities, its composition, structure, and dynamics. Once formed, the planets in turn influence the disk structure, potentially causing gaps, warps, and a variety of features that were revealed by millimeter-interferometry (e.g., Pérez et al. 2014; Andrews et al. 2018) and scattered-light imaging (e.g., Benisty et al. 2015, 2018; Avenhaus et al. 2018). However, these observing techniques cannot access the morphologies of the inner disk regions (located within 5 au around the central star) as they cannot reach a sufficiently high angular resolving power (~ milli-arcsecond) even thoughALMA starts to reach sub 10 au resolution for some targets. These inner regions are nonetheless of fundamental importance because they are the place where most of the planets are believed to form in or migrate to (Mordasini et al. 2009a,b; Alibert et al. 2011). Recent numerical studies have shown that the inner astronomical units of disks might harbor a strong viscosity transition between a so-called dead-zone, which is nearly impermeable to a magnetic field, and the hotter ionized disk front, showing active accretion powered by magneto-rotational instability (see, e.g., Kretke et al. 2009; Flock et al. 2016). This inner edge of the dead-zone might be the preferential location for dust pile-up and planetesimal formation through different instability mechanisms (Kretke et al. 2009; Flock et al. 2017; Ueda et al. 2019) and an effective planet filtering frontier (Faure & Nelson 2016).

Further evidence for the presence of perturbed inner disks comes from their known photometricand spectral variability associated with the presence of a circumstellar disk, in particular, in the infrared. Many authors have shown that 60–100% of stars with disks are variable (e.g., Espaillat et al. 2011; McGinnis et al. 2015; Poppenhaeger et al. 2015; Flaherty et al. 2016) and that the variability patterns arediverse (scale height variation, variable accretion, or ejection). Observations of the outer disk can also hint at variability in the close-in regions, as shown recently in the case of a rotating shadow in the disk of TW Hya (Debes et al. 2017) or HD 135344B (Stolker et al. 2017). In many cases, the period of the variability is compatible with Keplerian periods for radii up to the first inner astronomical unit around the star. Directly constrainingthe structure of the inner disks is an important key to better understand how planets form. This requires spatially resolving the inner astronomical units to possibly detect a companion or a disk perturbation (Flock et al. 2017)

The morphology of disks in the innermost regions can be constrained with multitelescope interferometry in the infrared-wavelength range. With a baseline of 140 m, such as the longest one currently available at the Very Large Telescope Interferometer (VLTI), and by observing in the near-infrared (1.65 μm), a resolution of ~0.3 au can be reached for young stars belonging to the closest star forming regions (like Taurus at 140 pc). The first interferometric near-infrared size measurements showed that the emission sizes are proportional to the square root of the luminosity of the central star (e.g., Monnier & Millan-Gabet 2002), suggesting that the near-infrared observations probe the dust sublimation front of the disk. Interestingly, some objects, which were studied more extensively, display a radially extended emission inside the theoretical dust sublimation radius that cannot be explained by an inner rim alone. The presence of a gaseous disk (Kraus et al. 2008) or refractory dust (Tannirkulam et al. 2008; Benisty et al. 2011) inside the conventional silicate sublimation radius was suggested to account for both an extended emission and a large flux excess at near-infrared wavelengths. Recent interferometric studies of the Br_γ hydrogen line have pointed towards the presence of a disk wind in several objects (e.g., Malbet et al. 2007; Garcia Lopez et al. 2015; Kurosawa et al. 2016). Such a wind can lift dust off from the disk and create an additional excess emission in the near-infrared (Bans & Königl 2012). In order to study the general properties of these regions, a comprehensive survey was needed.

Recently, Lazareff et al. (2017, hereafter L17) observed a large sample of intermediate-mass stars in the H-band continuum in a systematic and homogeneous fashion, with the Precision Integrated-Optics Near-infrared Imaging ExpeRiment (PIONIER) instrument, which was the first one to combine four telescopes at the VLTI. The observations of 51 Herbig AeBe stars were analyzed using parametric models to constrain the main characteristics of the dust emission at the au-scale. The main findings from this study are as follows: (1) The sublimation temperature of the dust emitting in the H-band is higher (~ 1800 K) than the classical silicate sublimation temperature (~1500 K); (2) the inner rims are smooth and radially extended, and (3) their brightness asymmetry is suggestive of inclination effects that make the rim at the far side brighter that the side closer to the observer. Those findings are confirmed by a survey of Herbig stars in the K-band performed with the GRAVITY instrument at the VLTI (GRAVITY Collaboration 2019).

There is ample observational and numerical evidence that inner disks might show a much more complex morphology than the one used in L17 to analyze interferometric observations. The putative presence of instabilities, vortices (Fung & Artymowicz 2014; Faure et al. 2015; Flock et al. 2017), or ring-like structures related to dead-zone dust pile-up self shadowing (Ueda et al. 2019) justifies the need for a different approach than parametric modeling. In this paper, we analyze a subsample of 15 objects observed by L17, with good enough (u, v)-coverage to allow for image reconstruction in order to obtain model-independent information on the disk morphology. This approach aims at investigating any structure that cannot be accounted for by simple geometrical model fitting and is a novel way to analyze the dataset without being dependent on a geometrical model. We therefore compared the results of the image reconstructions with those from parametric model fitting. We also searched for additional structures at the inner disk rim or outside of it that could be linked with dynamical instabilities, perturbations by a companion, detection of an inner cavity, or disk self-shadowing. For that purpose, we employed several image analysis tools that could be used in future optical and infrared interferometric imaging surveys. To perform image reconstruction, we used the Semi Parametric Approach for image Reconstruction of Chromatic Objects (SPARCO; Kluska et al. 2014) that enables the use of the spectral channels altogether to separate the central star from the environment.

The paper is organized as follows. We describe the observations and data selection in Sect. 2 and the image reconstruction process that is applied to each individual object in Sect. 3. We also present the reconstructed images of the extended structure in Sect. 4 with associated analysis graphs. In Sect. 5, we classify the objects and we discuss the origins of the detected asymmetries.

2 Observations and target selection

The dataset was extracted from the VLTI Large Program on Herbig Ae/Be stars (190.C-0963, PI: Berger). This program made use of PIONIER, which is a four-beam interferometric combiner working in the H-band at 1.65 μm (Le Bouquin et al. 2011). The presented dataset was observed between 2012 Dec. 19 and 2013 June 06 and comprises 15 targets (see Table 1). The associated datasets are shown in Figs. F.5–F.8. We refer to L17 for details regarding data reduction. This large program has been divided in two subprograms: survey observations on a large sample to assess the statistics (L17) and a detailed study on a small object sample to perform image reconstruction on the most resolved targets (this work). An interferometer measures partial information on the visibility, a complex number containing the fringe contrast, and phase that relates to the image via a Fourier transform. The dataset considered in this paper shows a clear decrease of visibility amplitudes with baseline, meaning that the targets are spatially resolved. The three spectral channels across the H-band display a similar trend where the blue channels have higher squared visibility values because there is a higher stellar contribution (e.g., Kluska et al. 2014).

The quality of an image reconstruction depends on the ratio between the smallest angular resolution element and the size of the emission and on the completeness of the (u, v)-plane. During the observations, targets that appeared to be well resolved were more extensively observed to optimize the (u, v)-coverage and perform image reconstructions. The (u, v)-coverages of the selected 15 targets are displayed in Fig. 1.

3 Image reconstruction approach

The image reconstruction is performed in a Bayesian framework where the best image minimizes a mathematical distance (), which is a combination of two terms: (1)

where is the distance to the data (here the χ²), is the regularization function (here the quadratic smoothness function defined as where S is a smoothing operator implemented via finite differences¹), and μ is the regularization weight that determines the balance between the two distances. The quadratic smoothness regularization is one of the most reliable regularizations (Renard et al. 2011). Also, the fact that this regularization is quadratic limits the effects of local minima. The value of the regularization weight (μ) depends on the (u, v)-coverage and on the morphology of the target (Renard et al. 2011).

The image reconstruction process is the one used in Kluska et al. (2016) and Hillen et al. (2016), coupling the SPARCO method (Kluska et al. 2014) with the multiaperture image reconstruction algorithm (MiRA, Thiébaut 2008). In this paper, we used a simple model for the star in the form of a point source with a visibility of unity across all baselines. The final complex visibility (V_tot) is therefore computed by this relation: (2)

where λ₀ = 1.65 μm is the reference wavelength, λ is the wavelength of observation (depending on the observed spectral channel), is the stellar-to-total flux ratio (between 0 and 1), d_env is the spectral index² of the circumstellar environment (d_env = −4 is the Rayleigh–Jeans regime; d_env = 1 corresponds to the spectral index of a black-body with a temperature of 1400 K), and V_img is the visibility of the image computed via a Fourier Transform.

It is important to note that and d_env are called the chromatic parameters and, together with the regularization weight μ, they have to be determined for each target separately. In order to do so, we briefly recall the main steps involved. First, to determine the stellar-to-total flux ratio and the spectral index of the circumstellar environment d_env, we performed a grid of reconstructions on those parameters. We then selected the parameters with the lowest (Appendix A). Second, to determine the best regularization weight μ of the quadratic smoothness regularization, we used the L-curve method (Appendix B; Hansen 2001). This method allows one to find the regularization weight that is optimal between a regime dominated by the data, but which is very noisy, and a regime dominated by the regularization, but for which the image does not fit the data. The two regimes form two asymptotes and the value in the knee of the curve represents the optimal value of the regularization. We investigated the effect of different regularization weights for two targets in our dataset with different angular sizes (HD 45677 and HD 100453, see Appendix B and Fig. B.2). The effect of the regularization weight does not change the overall image features. Last, a bootstrap process was then applied (Appendix C) to assess the significance of the image features.

For each object, Table 2 lists the pixel size, the chromatic parameters, the regularization weight log μ, and the corresponding reduced of the image³.

Once the images were computed, we used the following three tools to analyze them. First, asymmetry maps: those maps contain the part of the emission that is not point-symmetric with respect to the central star. Second, radial profiles: these profiles are made by azimuthally averaging the brightness of the deprojected image. Third, asymmetry factors (f_asym and ϕ_asym): they describe the strength and the orientation of the asymmetry with respect to the central star.

To justify the pertinence of these tools for our analysis, we tested them on a dataset generated from a radiative transfer model of a star with a disk (see Appendix E).

Fig. 1 (u, v) coverage of the targets. The colored circles represent the (u, v)-coverage; the colors correspond to the wavelength (blue for 1.55 μm, orange for 1.65 μm, and red for 1.75 μm).

4 Results

4.1 Reconstructed images

Figure 2 displays the reconstructed images for all the Herbig AeBe stars of our selection. Most of the targets have a that is around unity. Four of them, however, have relatively larger (HD 45677, MWC158, MWC297, and R CrA). This could be due to the break of a hypothesis used in the image reconstruction, for example that the flux in the image has the same spectral index or that the object is static during the observations. The break of the latter hypothesis is discussed in Sect. 5.7. Three targets, HD 45677, HD 98922, and R CrA, show the presence of a cavity (ring-like morphology) with a brightness deficit close to the central star. The other objects show an elongated morphology.

Five images exhibit flux at large distances from the star (HD 45677, MWC158, HD 98922, HD 100453, and HD 145718), but their significance per pixel is less than 3σ. This flux can either result from image reconstruction artifacts due to the lack of observations at small spatial frequencies (HD 45677) or from a drop in the squared visibility curve at very short baselines (MWC158, HD 98922, HD 100453, HD 145718) caused by an over-resolved component (Monnier et al. 2006; Pinte et al. 2008; Anthonioz et al. 2015; Kluska et al. 2016; Klarmann et al. 2017, L17). The morphology of the over-resolved component cannot be determined by the current dataset because of the lack of very short baselines. Without any additional spatial information (from, e.g., single telescope aperture masking observations), the reconstruction of the over-resolved flux distribution is highly influenced by the (u, v)-plane and is not treated further in this work.

Fig. 2 Reconstructed images of the Herbig AeBe targets separated in three groups. Top: centrally depressedobjects; middle: centrally peaked and asymmetric objects; bottom: centrally peaked and asymmetric objects. The cyan star marks the position of the central star, which has been subtracted, the black contours indicate the significance at 1 (dashed line), 3, and 5σ (solid lines), and the blue ellipse shows the size of the beam. The green dashed ellipses represent the theoretical dust sublimation radius (see Sect. 5.5).

4.2 Radial profiles

To probe the radial structure of our targets, we computed radial profiles for each image. To correct for the object orientations, we used the inclinations and position angles we obtained by fitting the images in the Fourier space (see Appendix D) We segregated the 15 targets into two categories based on the following profiles.

4.2.1 Centrally depressed

This category comprises the targets with an inner cavity seen in their radial profile (see Fig. 3). HD 45677, HD 98922, and R CrA show centrally depressed profiles, that is, increasing with radius towards a maximum and then decreasing. We note that the profile of HD 45677 displays a decrease before an increase of the brightness close to the star, indicating the presence of an additional component inside the central depression that is marginally resolved.

4.2.2 Centrally peaked

The emission is resolved but no inner cavity is detected at the angular resolution of our observations (the radial profile are monotonically decreasing). We note that 80% of the objects, that is, the 12 targets (HD 37806, MWC158, HD 100453, HD 100546, HD 142527, HD 144432, HD 144668, HD 145718, HD 150193, HD 163296, and MWC 297) display a monotonically decreasing profile.

4.3 Asymmetry analysis

To investigate any departures from axisymmetry in our targets, we computed the asymmetry maps for each image (Fig. 4). We also computed asymmetry factors f_asym and ϕ_asym (Table 3), which quantify the level of asymmetry and its orientation. These factors are described in Appendix E. We note that f_asym describes the amplitude of the asymmetry. A ϕ_asym close to 0° or ± 180° means a preferential orientation along the minor axis, whereas a ϕ_asym close to ±90° indicates an asymmetry along the major axis.

All of the images deviate, to some extent, from point symmetry. In order to study the significant asymmetries only, we focus on the objects that have a significant nonzero closure phase signal. It is important to note that 60% of our sample (9/15; HD 45677, MWC158, HD 98922, HD 100546, HD 144668, HD 150193, HD 163296, MWC297, and R CrA) present at least three nonzero closure phase measurements beyond 5σ, which indicate a departure from point-symmetry. In the rest of the asymmetry analysis, we focus on those objects only.

The objects that are centrally depressed present a clear asymmetry: they are brighter on one side of the major-axis than on the other. The radial profile of HD 45677 (see Fig. 3) has two radial components: a half-ring and an emission very close to the star. In the asymmetry map, the inner centrally peaked emission asymmetry is reversed with respect to the centrally depressed one. The half-ring brightness maximum is towards the northwest and the inner part maximum is in the opposite direction (i.e., towards the southeast; see Fig. 4 top-left).

Two centrally peaked objects (HD 144668 and MWC297) clearly show the same kind of asymmetry map as the centrally depressed objects even though their inner rims are not resolved. These asymmetry maps are similar to the model one (shown in Fig. F.4). Their morphology is therefore likely to be rim-like as well. MWC158, HD 100546, HD 150193, and HD 163296 have an irregular asymmetry map, which is probably due to a more complex inner rim morphology (as demonstrated for MWC158, Kluska et al. 2016). The asymmetry factors for these objects can be found in Table 3 and leads us to distinguish the following two trends.

Fig. 3 Radial profiles of the reconstructed images of Fig. 2. Black solid lines: radial profiles from images with 5σ error bars in color shades. The green line shows the azimuthally averaged size of the interferometric beam.

4.3.1 Asymmetry direction perpendicular to the emission position angle

The clearest example of such an asymmetry is the ring of HD 45677. This ring is very asymmetric when describing an arc in the image and in the asymmetry map. This is reflected in the asymmetry factors of f_asym = 0.45 ± 0.01 and ϕ_asym = −7° ± 1°, which have the highest amplitude among our targets. The other asymmetric objects that have an asymmetry angle (ϕ_asym) less than 30° from 0° are HD 98922, HD 144668, HD 163296, MWC297, and R CrA. R CrA has the second largest asymmetry amplitude factor (f_asym = 0.26 ± 0.02) and has a similar asymmetry map to HD 45677. The asymmetry map of HD 144668 displays flux on one side of the major-axis and the asymmetry amplitude is relatively high (f_asym = 0.24 ± 0.01). HD 98922 has the highest asymmetry angle of this subsample (ϕ_asym = 27° ± 6°). This is probably due to the presence of flux on the southeast part as displayed by the asymmetry map. MWC297 shows an asymmetry map and asymmetry factors that are compatible with a brightness enhancement on one side of the major axis. Finally, we note that HD 163296 has a low asymmetry angle, but it displays a relatively complex asymmetry map with several maxima and minima. It is interesting to mention that HD 142527 and HD 144432, despite having only two and one closure phase measurements at 5σ over zero, have a f_asym of 0.10 ± 0.01 and 0.11 ± 0.01 and a ϕ_asym of 171° ± 9° and 1° ± 11°, respectively. This, together with their asymmetry maps, would make them compatible with an asymmetry direction that is perpendicular to the emission position angle.

4.3.2 Asymmetry direction not oriented perpendicular to the position angle

We find that three out of the fifteen objects showing nonaxisymmetric emission display a skewness orientation that cannot be related to the disk position angle in a clear way. In this group, MWC158 is the source displaying the most extreme asymmetry factors (f_asym = 0.24 ± 0.01 and ϕ_asym = 72° ± 2°), but it was already identified as a temporally variable source in previous interferometric observations (Kluska et al. 2016). The emission of HD 100546 and HD 150193 show skewness that cannot be related to any of the two major or minor axes.

Fig. 4 Asymmetry maps of the targets classified the same way as in Fig. 2. The blue star represents the central star, the solid black contour represents the 3σ significance profile, and the dot-dashed line represents the orientation of the major-axis from parametric fitting. The targets with three closure phase points or more with a significance of 5σ are indicated as “asym”.

5 Discussion

To discuss and interpret the images, we classify the images (Sect. 5.1) before discussing individual characteristics of the circumstellar emission, such as the detection of an inner cavity (Sect. 5.2) or extended structures (Sect. 5.3). We also compare the most resolved targets with the parametric models from L17 (Sect. 5.4). Finally, we discuss some features of the reconstructed images with respect to the dust sublimation radius (Sect. 5.5), compare the asymmetry factors with those of axi-symmetric models (Sect. 5.6), discuss the origin of the noncentrosymmetry (Sect. 5.7) as well as inner and outer disk misalignment (Sect. 5.8).

5.1 Classification

To synthetize all the morphological features that we retrieved, we classify the targets based on the observational properties (see Table 4), that is, on the radial profiles and azimuthal analysis (Sects. 4.2 and 4.3). The centrally peaked objects that are marginally resolved appear to be centrosymmetric. The lack of significant nonzero closure phase signals can be interpreted as a strong flux contribution from the central star, which would smooth out the asymmetry from the disk. Another interpretation would be that a low disk inclination reduces the contrast between the far and the near sides of the rim or that we lack precision to detect the asymmetry. In support of the high stellar-to-total flux ratio hypothesis, we find that four out of the six symmetric objects have a strong contribution of the stellar flux ( > 50%), while all the noncentrosymmetric objects have a weaker contribution from the star ( < 50%). The objects showing a ring-like emission, which we associate to the dust sublimation rim, show azimuthal features that are compatible with inclination effects. Three centrally peaked objects (HD 144668, HD 163296, and MWC297) also show an inclination-like signature. They are likely to have an inner disk cavity as well but it is not resolved by our observations because of their angular size, radial extensions, and/or disk inclination. The (non)detections of inner cavities are discussed in more detail in Appendix F.2.

The objects showing irregular asymmetry do not show an inner cavity. This can be due to a lack of angular resolution (see Appendix F.2), but it might also indicate that a ring-like morphology is not possible with such irregular brightness distributions. Having a clear inner rim would lead to inclination effects in the case of a strong or moderate inclination angle. Two of these objects (HD 100546 and HD 150193) still need to be confirmed as irregular since their ϕ_asym has relatively large error bars (see Sect. 5.7). However, HD 100546 is surrounded by a perturbed circumstellar environment with gaps (Benisty et al. 2010; Tatulli et al. 2011) and possible bars and a warp in the first tens of astronomical units (Mendigutía et al. 2017; Walsh et al. 2017), which possibly explain the irregular structure we detect. MWC158 is already known to have a strongly variable environment (MWC158; Kluska et al. 2016).

Fig. 5 Inner hole detection in image reconstructions from (a) (continuous lines) synthetic observations of radiative transfer models and (b) (dots) our targets. Models with two radial distributions and two inclinations are included, see legend and text. Some dots are offset to avoid overlap.

5.2 Detection of inner cavities

The detection of a central cavity in the brightness distribution can be impacted by various factors, one of them being the ratio between the object’s size and the reached angular resolution. As a measure of the angular resolution, we used one-half of the smallest fringe angular wavelength achieved in the observations: θ_res = λ_min∕2B_max, and we examined the detectability of the cavity as a function of the ratio θ_1∕2 ∕θ_res where θ_1∕2 is the isophotal half-flux radius. The isophotal half-flux radius (θ_1∕2) is defined asthe size of an ellipse, which has the same orientation (inc and PA) as the imaged target and that contains half of the image flux.

We do not detect an inner cavity in twelve out of the fifteen targets despite them having a relatively extended circumstellar emission compared with the smaller angular scale resolved by the (u, v)-coverage (θ_res). In order to understand this, we made image reconstructions on radiative transfer models (using MCFOST; Pinte et al. 2006) of an A0 star surrounded by a disk having a sharp (Isella & Natta 2005) and a radially extended (or smooth; Tannirkulam et al. 2007)inner rim (more details in Appendix F) for two inclinations: one almost pole-on (inc= 18°) and a relatively high inclination (inc = 49°) corresponding to the higher inclinations in our sample (see Table 7). We show the detection of the inner cavity as a function of the θ_1∕∕2∕θ_res ratio on Fig. 5. We did not plot the detection against the inner cavity size because we want to compare it with our actual observations. In the actual dataset, in case of an undetected inner cavity, its size is unknown.

We can see that for pole-on orientations, the ratios of θ_1∕2 ∕θ_res are comparable though higher for smooth rim models as expected. For higher inclinations, the inner cavity of the models were resolved for ratios of θ_1∕2∕θ_res between twoand three with sharp rim models, which were resolved at lower values. The inclination of a disk, therefore, has a strong influence on the detection of the inner disk cavity.

We also plotted our targets on this diagram so we can discuss the inner cavity (non)detections. The inner cavity of two (HD 45677 and R CrA) of the most resolved targets were detected because of their size. HD 98922 is at the theoretical limit of detection, but its pole-on orientation and likely sharp inner rim contribute to the inner cavity detection. Several of the nondetections are probably due to a small size compared with the achieved angular resolution (HD 37806, HD 142527, HD 144432, HD 163296, and VV Ser, with θ_1∕2 ∕θ_res ≤ 1.5). Several targets are more resolved (θ_1∕2∕θ_res > 1.5), but most of them have high inclination (i > 45°). There are two targets (HD100453 and MWC297) that are relatively well resolved (θ_1∕2∕θ_res = 2.1) and with inclinations of 44.2° and 38.3°, respectively, for which no inner cavity is detected. For these targets, it is more likely that the radially extended emission from the inner rim causes the nondetection, which is compatible with conclusions from L17. For the other targets, the reasons for the nondetection are unclear (e.g., lack of angular resolution, absence of inner cavity, or high inclination) and longer baselines are needed to investigate the presence of inner cavities.

5.3 Detection of extended structures

In the reconstructed images (Fig. 2) in our sample, it appears that there are no significant features (above 3σ) outside the inner emission. From the azimuthally averaged radial profiles, there is also no significant emission outside the core emission. Here, we discuss the reason for the lack of extended structures in the images. To place an upper limit on the detection of such an emission, we added a ring to the restored images of our targets, generated a new dataset with the same (u,v)-plane, reconstructed an image, and plotted the radial profile. We did that with two different ring radii (5 and 10 mas) and different fluxes (from 1 to 15% of the total flux at 1.65 μm, the ring has the same spectral index as the image).

The radial profiles on reconstructed images using the dataset of one representative target, HD 144668, are displayed in Fig. 6 for ring-to-total flux ratios of 1 and 15%. The profiles with an added ring at 1% are within ±1 standard deviation of the original profile. This is true for most of the reconstructions up to 10% flux. For the 15% flux ring, the detection in the radial profiles are at ~5 and ~ 2σ for the 5 and 10 mas radius rings, respectively. No bump detected is detected in the radial profiles on our targets and the average detection limit between 5 and 15 mas is around 10% (excluding the most resolved targets HD 45677 and R CrA).We can therefore rule out any structure that is larger than the core emission that contributes to a ~ 10% flux level or less. We note that Ueda et al. (2019) predict that, under certain conditions, dust pile-up at the dead-zone inner boundary might cast a shadow on the inner astronomical units, which would lead to a second apparent ring-like structure of a few astronomical units in diameter. However, the predicted flux levels are on the order of the percent at best. We clearly lack the dynamical range to detect such a faint structure. The near-infrared emission in our images is therefore dominated by three sources: the star, the core emission close to the star, and, for some targets, an over-resolved emission for which we can only evaluate the integrated flux through model fitting (L17). Several possible origins can explain the latter, such as scattered light from the external parts of the disk (see, e.g., Benisty et al. 2010) or emission from quantum heated particles as pointed out by Klarmann et al. (2017).

Fig. 6 Radial profiles of reconstructed images on synthetic data where a ring has been added. The black curves show the radial profiles in the images restored from the original data. The blue and orange curves show the radial profiles in the restoredimages with an additional ring of radius 5 and 10 mas, respectively. The shaded areas show the limits at ±1 standard deviation of the profiles (of the same color). The standard deviation of the profile from theoriginal data is not shown. Top: the ring has a flux of 1% of the total flux at 1.65 μm. Bottom: the ring has a flux of 15% of the total flux at 1.65 μm.

5.4 Comparison with model fitting from L17

Whereas the least resolved targets do not show specific morphology features, some of the most resolved ones displayfeatures that are not included in the parametric models of L17. To make a fair comparison, we performed image reconstruction on datasets that were created from the best fits of L17 of ring models, including one so-called m1 azimuthal mode that corresponds to sinusoidal azimuthal modulation with a 2π-period. The model dataset has the same (u, v)-coverage as the real data on a given target and a Gaussian noise based on the error bars of the real dataset.

In Fig. 7, we compare the image reconstructions on the real dataset with the image reconstructions on the synthetic dataset on best fit m1-models from L17 for the nine best resolved targets (HD 45677, R CrA, HD 100453, HD 100546, HD 163296, MWC297, MWC158, HD 144668, and HD 98922). The m1-models can reproduce inclined inner rims of protoplanetary disks and especially any inclination effects. The orientations of the targets are similar to the ones of the images from the models. HD 98922 is the only target for which the orientation is not correct. This is due to the complex azimuthal structure of the extended emission where two maxima are present. A model with second order modulation retrieves an orientation of the disk that is compatible with what we find in this work (see Fig. D.1).

We can see that for HD 100453, HD 163296, and MWC297 the images are similar. For HD 45677, the overall structure is similar. However, because the disk rim is well resolved, we can discuss its details. The emission from the rim shows a dip in the north direction that was not recovered in the image from the model that has a smooth azimuthal brightness distribution. The dip is detected at 3σ significance. Also, the sides of the rim are sharper than in the model. This is not surprising as the azimuthal modulation in the m1-model isdone by a sinusoidal azimuthal modulation and may not reproduce azimuthal distributions of actual inner rims. Moreover, the flux close to the central star is not as strong as in the real image. The inner rim is therefore perturbed (for example by density or scale-height perturbations) and has a specific geometric profile that makes the closest part of the disk sharply self-shadowed. A hypothesis to test would be the presence of a close companion, which contributes to the marginally resolved flux close to the star, that carves and perturbs the disk rim, which is located three times farther than the theoretical sublimation radius (see Sect. 5.5 and Table 5). The actual origin of those features are worth investigating in future work dedicated to this target.

The R CrA model image does not show an inner cavity, whereas it is clearly present in the image on the real dataset. Recently, a binary inside the disk inner rim of R CrA was deduced from photometric times-series (Sissa et al. 2019). However, we do not see signs of disk truncation as the inner rim is compatible with dust sublimation. Also, the asymmetry is compatible with inclination effects which do not display perturbations due to the inner binary. For HD 100546, the image shows a smaller emission for the real dataset. There are also more structures in the image on the real dataset. HD 100546 is a transition disk and shows perturbation in the outer parts of the disk. It is further discussed in Sect. 5.7. MWC158 has an environment which is asymmetric and the position of the maximum is along the major axis. This asymmetry is reproduced by the model. However, in the real images there are structures that are not present in the smoothly changing azimuthal profile of the ring model. Finally, for HD 144668, we see that the sinusoidal azimuthal variation of the ring model produces emission coming from the disk that is too extended.

From this qualitative analysis, we can see that the m1 ring models reproduce the global target morphology. However, the sinusoidal azimuthal modulation and the presence of a single ring are not sufficient in reproducing the morphology of the several resolved targets. The structure of the inner rim is more complex and can be more easily reproduced by the image reconstruction technique. It is therefore likely that an additional physical process (e.g., disk instability or the presence of a companion; Flock et al. 2017) produces this kind of asymmetry. A first improvement to the study would be to ensure that observations of a given object on all telescope configurations are all carried out well within the expected Keplerian orbital time at the inner rim. This would allow a possible smoothing out of the asymmetric emission to be lifted. A second improvement would be to push for a higher angular resolution at VLTI using the longest baseline (≈ 200 m). This implementation is currently under study at the European Southern Observatory (ESO). Finally we note that some objects are sufficiently to the north so that they could be observed with VLTI and the Center for High Angular Resolution Astronomy (CHARA), which would provide a remarkable (u,v) coverage and angular resolution.

Fig. 7 Comparison between the image reconstructions on each target from the dataset (left panels) and the synthetic data based on the best-fit m1 model from L17 (middle panels) and the difference (right panels) in the 3σ significance area of the image. In the top-right corner of the model image reconstruction, there is the high resolution image of theparametric model. The green star represents the position of the central star.

5.5 Comparing the images to the theoretical dust sublimation radius

We can measure the sizes of circumstellar emission from the images by taking the half-light radius corrected for the disk orientation (see Appendix D). We assessed the accuracy of the size measurement against the true sizes from radiative transfer models in Appendix F by reconstructing images on the synthetic data generated from those models. The recovered sizes match the true sizes of the models with a sharp inner rim as long as the inner rim is resolved enough (i.e., θ_1∕2 ∕θ_res ≥ 1.2). However, we find a bias for sizes from extended inner rim models (see Appendix F). We tend to recover sizes that are smaller than the true sizes for very resolved targets (i.e., θ_1∕2 ∕θ_res ≥ 2.5).

Most of the targets (12/15) are relatively well resolved (i.e., θ_1∕2 ∕θ_res ≥ 1.3). In Fig. 8, we plotted the sizes of the targets against their luminosity for both this work and L17. For the imaged object, we see good agreement between both size determinations, that is, using parametric models and image reconstruction. For some targets (HD 145718 and R CrA), the difference is significant (ratio larger than 1.5 between the two sizes) and this is likely due to model bias (as the parameteric model for R CrA was an ellipsoid instead of a ring and the stellar-to-total flux ratios are different for HD 145718). We note that the size difference is in the opposite direction of the size bias from image reconstruction that we identified in Appendix F.

In order to have a reference to interpret the sizes of the circumstellar emission from the reconstructed images, we computed the theoretical dust sublimation radii using: (3)

(see, e.g., Monnier & Millan-Gabet 2002; Kama et al. 2009, L17), where C_bw is the back-warming coefficient that we fixed to four, ϵ is the cooling efficiency dust grain that we set to one, L_* is the stellar luminosity that we took from L17, and T_sub is the sublimation temperature that we set to 1500 K. Since the physical size and the square root of the stellar luminosity are proportional to the distance, we can estimate the angular radius of sublimation without being sensitive to the distance error for a given target. The angular sublimation radius for each target is defined in Table 5 and is shown on the reconstructed images (Fig. 2). However, this theoretical sublimation radius can only be used as a reference since the different parameters involved (C_bw, ϵ, or T_sub) are not absolutely constrained and they may vary from object to object.

While most of the targets are reasonably aligned with the reference size-luminosity relation defined by Eq. (3), two targets are outliers: HD 45677 and MWC297. HD 45677 has a resolved cavity and has a size that is more than three times larger than the theoretical dust sublimation radius. This could be due to an inner binary truncating the disk, for example, even though no evidence for such a companion has been found in literature. For the other two targets that have a resolved inner cavity (HD 98922 and R CrA), the dust sublimation radius matches the radius of the centrally depressed emission well. MWC297 shows that its emission is more compact than the typical dust sublimation radius. It is possible that for this target, which is a B1.5 star (Fairlamb et al. 2015), the emission is not coming from the dust sublimation but from another source, for example, free-free emission as it is the case for Be stars. For all the other targets, the size-luminosity diagram confirms the findings of L17, where the targets near-infrared circumstellar sizes are proportional to the square root of stellar luminosity.

Fig. 8 Size-luminosity relation. The green circles are the sizes measured from images, red crosses are the sizes from parametric models as determined in L17, and cyan triangles are the physical spatial resolution reached by our observations. The two size estimations and the spatial resolution are linked by a vertical dotted line. The dashed gray line represents the theoretical sublimation radius.

Fig. 9 fasym versus | arcsin(sin(ϕasym))| for observed (green) and model (blue and orange) objects.

5.6 Asymmetry orientation compared with axisymmetric models

We have investigated the asymmetry that we could find by analyzing images that were reconstructed on synthetic datasets. We used the synthetic datasets generated from MCFOST radiative transfer models used in L17. Those models are described in Appendix F.

We selected data that have a significant closure phase signal. We then applied the same image reconstruction approach to obtain images for each of the model inclinations and distances. In Fig. 9, we present the asymmetry factors found from the models. Instead of plotting ϕ_asym, we plotted |arcsin(sin(ϕ_asym))|, which is the absolute angle between the minor-axis and the asymmetry, independent of its absolute orientation. We can see that | arcsin(sin(ϕ_asym))| is below 20° for most of the models. The sharp models are more represented as they have a higher closure phase signal. They also have a larger f_asym.

We also placed the observed asymmetric objects on this diagram for comparison purposes. First, we can see that the objects that have an inclination-like asymmetry are in the left part of the diagram and are in the area that is populatedwith the axisymmetric models. While MWC158 cannot be reproduced by the models, HD 100546 and HD 150193 are located outside of the parameter space reproduced by the models and they are therefore classified as asymmetric in our analysis. However, they have a large error bar. HD 98922 has an asymmetry angle (ϕ_asym) of 27°±6° and is located outside the range span by the axisymmetric models. This is due to the asymmetry located at the southeast of the star. This is a hint for disk asymmetry that needs to be followed up by further observations. Further investigation of these targets would help to clearly identify if those stars are surrounded by nonaxisymmetric disks, especially since we note that six out of ten objects that are the most resolved (θ_1∕2∕θ_res ≥ 1.5) have an f_asym > 0.1, meaning that potentially about half of the unresolved targets could have an asymmetric morphology. This is in contrast to the five less-resolved targets that all have f_asym ≤ 0.11.

5.7 Origin of the noncentrosymmetry

An obvious cause for the noncentrosymmetry of the stellar environment is the disk inclination. This seems to be the case for five targets (HD 45677, HD 98922, HD 163296, MWC297, and HD 144668). However, in this section, we further investigate the origin of environment noncentrosymmetry that could not be due to inclination effects.

The disk dust sublimation rim can be perturbed by instabilities due to irradiation from the central star (Fung & Artymowicz 2014) or by Rossby wave instability occurring at the inner edge of the dead zone (Flock et al. 2017). Any perturbation of the inner disk affecting its vertical structure should generate a detectable departure from central-symmetry that would add to the natural one arising due to the inclination effect. It’s orbital motion around the central star should leave a time-dependent signature. Such a bright spot on the inner disk rim produces an asymmetric image but also orbits the central star.

The three objects with irregular departure from axial symmetry discussed in Sect. 5.1 can therefore trace dynamical phenomena in their disk, which induce temporal variability. A clear and strong change in the disk morphology of MWC158 was previously noted (Kluska et al. 2016).

HD 100546 is a well known transition disk, with evidence for asymmetric features, such as spiral arms in the outer disk (seen in scattered light), and for a possible protoplanet candidate (e.g., Quanz et al. 2013, 2015; Currie et al. 2015; Garufi et al. 2016; Follette et al. 2017; Rameau et al. 2017) that could propagate to the inner regions. HD 163296 is known to host a jet that launches a knot every 16 yr (Sitko et al. 2008; Ellerbroek et al. 2014). The origin to this regular periodicity can be a perturbation by a companion (Terquem et al. 1999; Whelan et al. 2010; Estalella et al. 2012) or a warped inner disk (Lai 2003). More recently, a slight misalignment of the inner disk was suggested by modeling the direct imaging observations (Muro-Arena et al. 2018). We also detect a disk misalignment for HD 150193 (see Sect. 5.8), suggesting a link with the perturbations we detect in the inner disk.

We discuss how an asymmetry located at the detected circumstellar structure would move during the observing time span. Because several different configurations on the VLTI are required to cover the (u, v)-plane, interferometric measurements can be collected during a period of a few weeks. We therefore compare the temporal span of our observation with the Keplerian period at the position of the detected asymmetric circumstellar structure (Table 6). To convert the angular distances to the physical ones, we used distance from Vioque et al. (2018) who used parallaxes from the second data release of the Gaia mission (Gaia Collaboration 2016, 2018). Most of the objects completed at least half of their orbit. Therefore, our observations sample a non-negligible part of the orbital time-scale of these objects. The exceptions are HD 45677 and MWC297 for which the observation samples less than 15% of the orbital period. The irregular features we observe can be orbiting the star during the time span of the observations (as it was the case for MWC158). To explore this hypothesis, well-sampled interferometric observations both in time and spatial frequencies are necessary to follow the evolution of the disk morphology.

5.8 Innerand outer disk misalignment

Several disks around young stars display shadows that are attributed to misalignment between the inner and outer parts of the disk (e.g., Marino et al. 2015; Benisty et al. 2017; Casassus et al. 2018). These disks are mainly transition disks, whichhave a gapped disk structure. In these gapped disks, the observed dark spots on the outer disk are interpreted as shadows, which are being cast by the inner disk on the inner edge of the outer disk. The origin of such a misalignment is still under debate, but the main invoked cause is the presence of a companion in the gap (see e.g., Nealon et al. 2019).

A direct way of probing the disk misalignment is to compare the orientation of the inner disk from this work with values of the outer disk coming from near-infrared scattered light imaging or millimeter interferometric observations. We computed the angle difference (α) between the normal vector of the inner disk we determined and the one from the outer disk that we found in literature (see Table 7; Marino et al. 2015). There is an angle degeneracy with the inclination as one can only measure the length ratio between the major and minor axis and not the inclination that can be negative or positive. We therefore computed both angles (α₁ and α₂). In order to conclude on a misalignment between the inner and the outer disks, we propagated the errors on inclination and position angles to have an error on the angle difference between the two disks axes (σ_α). When the errors on the inclination and position angles were not known, we assumed an error of 5° (Min et al. 2017).

We find orientations of the outer disk for five targets. Four out of five of these targets show significant misalignment between the outer and the inner disk. Disk misalignment was already observed for HD 100453 (Benisty et al. 2017; Long et al. 2017; Min et al. 2017) and HD 142527 (Marino et al. 2015). In the following, we discuss the misalignments per target.

5.8.1 HD 100453

For HD 100453, a misalignment angle (α) of ~ 72° was deduced (Benisty et al. 2017). This angle is compatible with one of our solutions (α₂ = 73.2 ± 9.6).

5.8.2 HD 142527

From the shadow modeling of HD 142527, a misalignment angle of 70° ± 5° was suggested (Marino et al. 2015). Our measurement does not confirm this value. A modeling of the whole system, including orientations of the inner disk from this work, is needed to assess if the shadows can be due to the inner disk misalignment that we detect (see, e.g., Nealon et al. 2019).

5.8.3 HD 100546

HD 100546 is a well studied transition disk with evidence of a complex environment in the inner disk. From polarimetric direct imaging observations in the visible with SPHERE, Mendigutía et al. (2017) observed bars across the gap. ALMA observations of CO display a strong residual signal once a Keplerian model was fit (Walsh et al. 2017), suggesting out-of plane kinematic gas motion occurring in the inner parts of the disk. A possible warping of the disk was previously suggested (Pineda et al. 2014). We find misalignment angles of 14.6° and 83.0° ± 2.6°. For the latter value, we would expect to detect sharp shadows but no evidence for them has been found so far (Sissa et al. 2018). Nevertheless, the orientations of the inner disk that we derived are not compatible with the inner feature seen with ALMA or the bar morphology seen with SPHERE.

5.8.4 HD 150193

HD 150193 is also such an object with suspected nonaxisymmetric inner disk morphology. The misalignment is about 28.9° or 70.0° ± 10.3°. Unlike the other targets showing misaligned disks, HD 150193 is not a known transition disk (e.g., Banzatti et al. 2018). However, a stellar companion was detected outside the disk (Fukagawa et al. 2003; Carmona et al. 2007; Garufi et al. 2014; Monnier et al. 2017), which is able to dynamically perturb it as in the case of HD 100453.

5.8.5 HD 163296

Finally, we do not detect a significant misalignment of the disk of HD 163296. A misalignment of ~ 3° was suggested by a recent modeling work to reproduce both direct imaging SPHERE and millimetric ALMA observations (Muro-Arena et al. 2018). Our observations are not precise enough (σ_α = 9.4) to confirm this interpretation. Moreover, this object can have a nonaxisymmetric morphology at the inner rim (see Sect. 5.7).

6 Summary and conclusions

We summarize here our most important findings from this work:

• 1.

  We probed the inner astronomical units of the circumstellar environments of Herbig Ae/Be stars from the VLTI/PIONIER survey by applying the SPARCO image reconstruction method on 15 objects. It is the first imaging survey of the inner parts of protoplanetary disks through near-infrared interferometry. This imaging method is motivated by a need to investigate the morphologies of the disk inner regions and it is complementary to model fitting (L17) but it is less dependent on a strong prior hypothesis.

• 2.

  We note that 20% (3/15) of our targets show a ring-like morphology with a detected inner cavity. The rest of the objects (80%) are consistent with a centrally peaked morphology at the observed angular resolution.

• 3.

  We classified the objects into different categories based on their radial and azimuthal morphology. We notice that, as expected, the less resolved objects and the ones where near-infrared flux is dominated by emission from the star do not show strong closure phase signals. The objects showing ring-like morphologies are dominated by inclination-like asymmetric emission.

• 4.

  Our comparison of the most resolved targets with their geometric models from L17 indicates a more complex morphology for seven targets.

• 5.

  We ruledout structures outside the very inner disk emission. We worked out an average upper limit of 10% of the flux at 1.65 μm that could be located between 5 and 15 mas from the star.

• 6.

  We thenanalyzed the retrieved images by producing asymmetry maps and radial profiles. We note that 60% (9/15) of our targets show noncentrosymmetric emission, and 40% of them (6/15) are consistent with inclination-like noncentrosymmetry. The other ~20% (3/15) display irregular features that are linked with nonaxisymmetric and possibly time-variable morphologies. Those morphologies could be due to a disk warp or instabilities at the disk inner rim, such as Rossby wave instability, creating a vortex or the presence of a companion.

• 7.

  We probe misalignments between the inner and the outer disk regions. We confirm the misalignments of HD 100453 and HD 142527, even though, for the latter, the orientation of the inner disk is different from what was suggested from the shadows cast on the outer disk. We confirm a potential disk warp in HD 100546. We also suggest that there is a disk misalignment in HD 150193.

• 8.

  Finally, we raise the question of the non-negligible dynamical evolution of the objects during interferometric observations that can disrupt the process of an image reconstruction if the time span of observation is comparable to the dynamical time-scale.

The simultaneous availability of PIONIER, GRAVITY, and MATISSE at VLTI offers the opportunity of multiwavelength studies covering the H, K, L, M, and N bands. This allows the observer to probe different regions in the disk, both vertically and radially. One of the first outcomes of such a study would be to determine the degree of flaring of the inner rim as predicted by Flock et al. 2017, for example, in relation with various processes at play (sublimation, wind, magnetic field). Moreover, the combination of visibility and SED fitting from 1.5 to 13 μm is a powerful tool to determine dust composition, size, and mineralogy. It would be particularly interesting to confirm L17’s hint that a significant fraction of dust responsible for near-infrared emission might be of a refractory nature. Exploiting the 200 m available at VLTI will allow one to locate the disk inner rim and compare it to the dust composition. Additionally, GRAVITY’s polarimetric capability should be used to try to map the scattered light and further constrain the nature of the dust grains. It might even be possible to use polarization to enhance the contrast of the environment with respect to the star and better reveal the perturbed disk. We also note that combining time-resolved observations in the continuum, in the Brackett Gamma line together with spectro-polarimetry should be a powerful tool to test the magnetospheric accretion scenario. In the case of Herbig AeBe stars, as for T Tauri stars, one might expect that the inner dust is lifted into the magnetospheric funnels, thereby offering the possibility to relate continuum asymmetry with the exquisite line spectro-astrometry that could be related to possible accretion impacts andmagnetic field topology. Finally, a dedicated search for putative inner companions responsible for disk perturbations should be carried out by pushing the dynamic range of the interferometric observations (better precision, many observing points) in combination with the high resolution infrared spectrometers that will be available soon.

Acknowledgements

We want to acknowledge the referee for his/her comments. This work is supported by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR-10-BLAN-0511). JK acknowledges support from an Marie Curie Career Integration Grant (SH-06192, PI: Stefan Kraus) and a Philip Leverhulme Prize (PLP-2013-110). J.K. acknowledges support from the Research Council of the KU Leuven under grant number C14/17/082. F.B. acknowledges funding from National Science Foundation awards number 1210972 and 1616483. S.K. acknowledges support from an ERC Starting Grant (Grant Agreement No. 639889). F.M. acknowledges funding from the EU FP7-2011 under Grant Agreement No 284405. F.M. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013. C.P. acknowledges funding from the Australian Research Council via FT170100040, and DP180104235. PIONIER is funded by the Université Joseph Fourier (UJF, Grenoble) through its Poles TUNES and SMING, the Institut de Planétologie et d’Astrophysique de Grenoble, the “Agence Nationale pour la Recherche” with the program ANR EXOZODI, and the Institut National des Science de l’Univers (INSU) via the “Programme National de Physique Stellaire” and “Programme National de Planétologie”. The authors want to warmly thank all the people involved in the VLTI project. This work is based on observations made with the ESO telescopes. It made use of the Smithsonian NASA Astrophysics Data System (ADS) and of the Centre de Données astronomiques de Strasbourg (CDS). All calculations and graphics were performed with the open source software Yorick. This research has made use of the Jean-Marie Mariotti Center ASPRO2 and SearchCal services co-developed by CRAL, IPAG and FIZEAU.

Appendix A: Determination of the chromatic parameters

The determination of the chromatic parameters ( and d_env) of each object was obtained by a selection a posteriori of the reconstructed images: a 100 × 100 grid was built using image reconstructions with different chromatic parameters. For each couple of chromatic parameters, the value of the total function f was minimized by the image reconstruction algorithm. The probability density function was then normalized to one because the densities are negligible at the edges. The grids are displayed in Fig. A.1.

Fig. A.1 Significance maps for the chromatic parameters of each object. Contours indicate 5σ likelihood of the parameters for the actual dataset (black). Blue and red triangles represent the chromatic parameters from L17 for ring models with first order and second order modulations, respectively.

One bias was identified with the SPARCO method: when the environment is marginally resolved, the parameter that corresponds to the unresolved contribution can carry the flux of the star and a part of the flux from the unresolved environment. Therefore, it yields an overestimation of the stellar-to-total flux ratio () and of the environment spectral index (d_env) since the assumption that all the flux in the parametric model is in the Rayleigh-Jeans regime is no longer correct. It leads to a whole range of acceptable values (see Fig. A.1; Kluska et al. 2014, 2016). An interesting case is HD 100453, where the parametric models from L17 predict a ring-like structure, but the reconstructed images from these models do not show such a feature. It is important to note that (1) if the data do not resolve the inner cavity, it is not seen in the image, as illustrated in the imagereconstruction of the model shown in Appendix F and that (2) the reconstruction is not based on a simultaneous photometric measurement. However, trying an image reconstruction using the chromatic parameters from L17 did not result in an image with an inner cavity. This object is very similar to the case of a model treated in Appendix F that has a transition radial profile between clearly resolved inner cavities and unresolved ones.

We also noticed a change in the allowed values of the chromatic parameters with the strength of the regularization. An increase in μ favors a smooth image and quenching noise, but it potentially fills a central cavity in the circumstellar brightness distribution. A side effect is that the extra flux assigned to the center of the circumstellar distribution is taken from the stellar flux, that is, a decrease of . Furthermore, because of the degeneracy in the (, d_env) parameter pair, which is visible in the oblique likelihood contours of Fig. A.1, the decrease of is coupled with an increase of d_env. We show, however, below and in Fig. B.2, that the actual impact of these effects is negligible.

However most of the chromatic parameters are well constrained and they are in agreement with a photometric fit of the parameters. Unfortunately, the photometric observations were not usually simultaneously carried out with our observations; therefore, it is difficult to state whether the determination of the parameters from the interferometric data is wrong or if the object has varied. In the paper, we keep the chromatic parameters determined from the interferometric data alone.

Appendix B: Regularization weight determination

The optimal regularization weight depends on the (u, v)-coverage and the complexity of the shape of the object. It can be determined using the L-curve method (Renard et al. 2011; Kluska et al. 2016). Once the chromatic parameters were fixed (see previous section), image reconstructions were performed to build a grid of the regularization weight μ. In the plot representingf_rgl in function of f_data, two regimes appear to create asymptotes that give an “L” shape to the plot. The first regime is dominated by the likelihood and the second by the regularization. The optimal regularization weight corresponds to the one located at the bend of the L-curve. The L-curves for each object are displayed in Fig. B.1. If the found regularization weight differs significantly from the one used to determine the chromatic parameters, then we have iterated by redoing the chromatic parameters determination with the new regularization weight and the L-curve was checked again. All the parameters used for the image reconstruction including the chromatic parameters and the regularization weights used for each object are summarized in Table 2.

Fig. B.1 L-curves for all objets. The bend of the curve corresponds to the optimal regularization weight.

Fig. B.2 Illustration of the effect of a slight variation of the regularization weight (μ) around the chosen value for HD 45677 (top) and HD 100453 (bottom).

We also performed image reconstructions with slightly different parameters (see Fig. B.2). There are some slight changes in the morphology of the images as the less regularized images appear less smooth that the most regularized once, which is as expected. The sizes of the image do not change for the most resolved targets (HD 45677 and HD100453). This is the degeneracy between the size, the stellar-to-total flux ratio, and the spectral index that was also seen for geometrical models (L17). The details we see for HD 45677 (dip of flux in the northern part of the ring and sharp turns at the extrema of the circumstellar emission) stay in the image even if they appear dimmer in the most regularized case as the image becomes over-regularized.

Appendix C: Image pixel bootstrap

We estimated the significance of each pixel in the image using a bootstrap method: 500 datasets were generated from the original one by drawing the baselines and the triangles one by one. Each baseline (or triangle for the closure phases) can be redrawn several times. We ended this process when we had drawn as many baselines and triangle as are in the original data. Image reconstruction was then performed for all of these datasets.

The average flux value of the pixel was divided by its standard deviation in order to give an estimation of its significance. The average images are less sensitive to the (u, v) coverage and the rest of this work is based on them. The difference between the of the averaged images from the bootstrap and the reconstructed images directly from the dataset is negligible.

The flux of the extended structures, even if it is not significant pixel by pixel, is required to correctly fit the data because it contributes to the short baselines. The shapes of these extended structures vary depending on the (u, v) coverage and it is therefore not reliable.

Appendix D: Determining the inclination and position angles from the reconstructed images

To determine the inclinations and position angle directly from the images, we fit a two-dimensionnal Gaussian to the amplitude of the Fourier transform of the image. This method is less model-dependent than fitting a model to the dataset itself since, depending on the object, different models can be fitted and account for the object morphology (e.g., a Gaussian, a ring, ring modulations, shift between the star, and the environment). As the central part of the image is similar if it is a ring or a Gaussian, the obtained orientations (inclination and position angle) are not dependent on the exact type of morphology. Whereas the inclination of the Gaussian is directly the inclination of the image, the position angle has to be shifted by 90°. We comparethe values of the inclination and position angle with the those from ring models with first order modulation from L17. The inclinations obtained from the image are similar to those obtained by parametric fitting in L17 (see Fig. D.1).

Fig. D.1 Comparison between inclinations (left) and position angles (right) obtained from the image and parameteric fitting in L17. The dashed line represents the one-to-one line.

Appendix E: Generic tools to analyze the images: Asymmetry maps, azimuthal and radial profiles, and asymmetry factor

To analyze the reconstructed images, we developed several tools that allow us to quantify the morphology displayed by the image. These tools are described below.

E.1 Radial profile

The radial profile is the azimuthally averaged profile of the image, which takes the disk orientation on the plane of sky into account. In order to do so, we defined concentric ellipses, which were oriented using the inclination and position angle of the disk model. In the case of the observational data, we fit the inclination and position angles in the Fourier space of the reconstructed images (see Appendix D). More specifically, the radial profiles were obtained by azimuthally averaging the flux of 150 consecutive annuli ranging from the first pixel to the outer edge of the image.

Fig. E.1 Sketch explaining the asymmetry factor computation. We note that Δfmajor (and Δfminor, respectively) is the difference of flux between the two sides of the major-axis (the minor-axis respectively).

E.2 Asymmetry map

Asymmetry maps consist in subtracting the 180 degrees-rotated image to the image itself and taking the positive part (for more information see Kluska et al. 2016). The resulting image shows all the flux that is not point-symmetric and that contributes to a nonzero closure phase signal.

E.3 Asymmetry factors (f_asym and ϕ_asym):

To quantify the asymmetry in both amplitude and orientation, we used the asymmetry factors (f_asym and ϕ_asym) that we define as follows (as inspired by Colavita 1999):

where A, B, C, and D are the different quadrants as indicated in Fig. E.1. We note that ϕ_asym is defined so that ϕ_asym = 0° or ±180° indicates asymmetry along the minor axis and ϕ_asym = 90° or − 90° indicates an asymmetry along the major axis in the anticlockwise and clockwise direction, respectively. If f_asym is close to unity, the amplitude of the asymmetry has a maximum of double the average brightness and a minimum of zero. On the other hand, if f_asym is low, the asymmetry amplitude is not strong or is more complex to be captured by these factors.

Appendix F: Test on synthetic data from a radiative transfer model of a disk

As an example of the tools we have defined, here, we show their application on radiative transfer models. We used the following two radiative transfer models (computed with the radiative transfer code MCFOST; Pinte et al. 2006, 2009): one with a sharp inner rim (Isella & Natta 2005) and one with a radially extended inner rim (smooth rim) as defined by Tannirkulam et al. (2007). Those models are described in L17. Their goal is to produce models of the disk inner rim with a radially sharp or extended near-infrared emission from the inner edge and not necessary to claim that they are representing the reality. Both disk models are axisymmetric and the inner disk rim is set by the radiative transfer model with an inner rim at 0.1 au. The difference between the two models is the scale height (0.037 au for sharp rim and 0.05 au for smooth rim at 1 au) and the dust size distribution with a single size of 1.2 μm for the sharp rim model and sizes between 0.18 and 2.23 μm for the smooth rim model. We can see the differences between the two models in Fig. F.1, which shows the images at 1.65 μm at two different inclinations. The smooth rim models clearly have more extended emission from the inner rim.

Fig. F.1 Images of the radiative transfer models at 1.65 μm at an inclination of 49 degrees. Left: sharp rim model, right: smooth rim model. The position of the star, which was subtracted from the image, is represented by the blue cross.

F.1 Inner cavity detection and radial profile

We performed image reconstructions on data generated from the radiative transfer models using the same methodology as for the observed targets. We used the (u, v)-coverage obtained on HD 144668 and simulated noise on the complex visibility. We reconstructed images of the model taking different distances and inclinations. We then investigated the presence of a central cavity, the radial profile, and sizes. The results of the detection of the central hole are presented in Fig. 5. The radial profiles (Fig. F.2) were used for this inner hole detection. When the radial profile is not monotonically decreasing, then the inner hole is considered to be detected. We can see that the distance, and therefore the angular size of the source, influence the radial profile and the detection of the inner hole. We also see that for higher inclination, it is more difficult to detect the presence of an inner hole since the minor axis of the projection of the inner rim on the plane of the sky is too small to be resolved, whereas the major axis is resolved.

F.2 Size determination

We wanted to test the size determination from the images (described in Appendix D) and test if it is reliable to determine the emission size. We compared the determined sizes from the reconstructed images of the models with different inclinations that are representative of our targets (a low inclination with inc = 18° and moderate with inc = 49°) and distances in Fig. F.3. We did this with the two models. We can see that for the sharp rim models, the size determined from the image matches the model size well for both inclinations as long as the size is resolved by the interferometer. For the smooth rim model, however, the sizes determined from the images match for relatively small sizes (θ_1∕∕2 < 2.5 mas), but they are then biased towards smaller sizes from θ_1∕∕2 > 2.5 mas. The investigation of the origin for such a bias is beyond the scope of this paper, but it is useful to keep it in mind.

Fig. F.2 Radial profiles from image reconstructions on the sharp rim models for high inclination (49°). The solid lines represent the radial profiles for models at different distances, which are color-coded. The green horizontal line above the profiles represents the beam size.

F.3 Asymmetry map

In Fig. F.4, we present an example of an asymmetry map from a reconstructed image (with a of 1.2) of a radiative transfer model that has a sharp rim and an inclination for 49° and at a distance of 100 pc. The asymmetry factors for this model are f_asym = 0.14 ± 0.01 and ϕ_asym = 5°±1°.

Fig. F.3 Comparison of the emission size measured in the image with the theoretical size. The red area indicates where the angular resolution given by the (u, v)-coverage is not sufficient in resolving the emission. The dashed line represents the 1:1 relation between the theoretical size and the size determined from the image.

Fig. F.4 Top: image reconstruction of a sharp model with an inclination of 49° degrees and located at 100 pc. The blue star represents the position of the star. The contours are 1 (dotted), 2, and 3-sigma pixel significance levels. Bottom: corresponding asymmetry map with the contour representing the 2σ significance level.

Fig. F.5 Squared visibilities of the dataset and the images. For each target, the top panel shows the squared visibilities for the data (black circles) and the image (blue triangles).Bottom panel: residuals normalized by the error bars on each data point.

Fig. F.6 Same as Fig. F.5.

Fig. F.7 Closure phases of the dataset and the images. For each target, the top panel shows the closure phases for the data (black circles) and the image (blue triangles). Bottom panel: residuals normalized by the error bars on each data point.

Fig. F.8 Same as Fig. F.7.

References

All Tables

All Figures

	Fig. 1 (u, v) coverage of the targets. The colored circles represent the (u, v)-coverage; the colors correspond to the wavelength (blue for 1.55 μm, orange for 1.65 μm, and red for 1.75 μm).
In the text

Fig. 2 Reconstructed images of the Herbig AeBe targets separated in three groups. Top: centrally depressedobjects; middle: centrally peaked and asymmetric objects; bottom: centrally peaked and asymmetric objects. The cyan star marks the position of the central star, which has been subtracted, the black contours indicate the significance at 1 (dashed line), 3, and 5σ (solid lines), and the blue ellipse shows the size of the beam. The green dashed ellipses represent the theoretical dust sublimation radius (see Sect. 5.5).

In the text

	Fig. 3 Radial profiles of the reconstructed images of Fig. 2. Black solid lines: radial profiles from images with 5σ error bars in color shades. The green line shows the azimuthally averaged size of the interferometric beam.
In the text

	Fig. 4 Asymmetry maps of the targets classified the same way as in Fig. 2. The blue star represents the central star, the solid black contour represents the 3σ significance profile, and the dot-dashed line represents the orientation of the major-axis from parametric fitting. The targets with three closure phase points or more with a significance of 5σ are indicated as “asym”.
In the text

	Fig. 5 Inner hole detection in image reconstructions from (a) (continuous lines) synthetic observations of radiative transfer models and (b) (dots) our targets. Models with two radial distributions and two inclinations are included, see legend and text. Some dots are offset to avoid overlap.
In the text

Fig. 6 Radial profiles of reconstructed images on synthetic data where a ring has been added. The black curves show the radial profiles in the images restored from the original data. The blue and orange curves show the radial profiles in the restoredimages with an additional ring of radius 5 and 10 mas, respectively. The shaded areas show the limits at ±1 standard deviation of the profiles (of the same color). The standard deviation of the profile from theoriginal data is not shown. Top: the ring has a flux of 1% of the total flux at 1.65 μm. Bottom: the ring has a flux of 15% of the total flux at 1.65 μm.

In the text

Fig. 7 Comparison between the image reconstructions on each target from the dataset (left panels) and the synthetic data based on the best-fit m1 model from L17 (middle panels) and the difference (right panels) in the 3σ significance area of the image. In the top-right corner of the model image reconstruction, there is the high resolution image of theparametric model. The green star represents the position of the central star.

In the text

	Fig. 8 Size-luminosity relation. The green circles are the sizes measured from images, red crosses are the sizes from parametric models as determined in L17, and cyan triangles are the physical spatial resolution reached by our observations. The two size estimations and the spatial resolution are linked by a vertical dotted line. The dashed gray line represents the theoretical sublimation radius.
In the text

	Fig. 9 fasym versus | arcsin(sin(ϕasym))| for observed (green) and model (blue and orange) objects.
In the text

	Fig. A.1 Significance maps for the chromatic parameters of each object. Contours indicate 5σ likelihood of the parameters for the actual dataset (black). Blue and red triangles represent the chromatic parameters from L17 for ring models with first order and second order modulations, respectively.
In the text

	Fig. B.1 L-curves for all objets. The bend of the curve corresponds to the optimal regularization weight.
In the text

	Fig. B.2 Illustration of the effect of a slight variation of the regularization weight (μ) around the chosen value for HD 45677 (top) and HD 100453 (bottom).
In the text

	Fig. D.1 Comparison between inclinations (left) and position angles (right) obtained from the image and parameteric fitting in L17. The dashed line represents the one-to-one line.
In the text

	Fig. E.1 Sketch explaining the asymmetry factor computation. We note that Δfmajor (and Δfminor, respectively) is the difference of flux between the two sides of the major-axis (the minor-axis respectively).
In the text

	Fig. F.1 Images of the radiative transfer models at 1.65 μm at an inclination of 49 degrees. Left: sharp rim model, right: smooth rim model. The position of the star, which was subtracted from the image, is represented by the blue cross.
In the text

	Fig. F.2 Radial profiles from image reconstructions on the sharp rim models for high inclination (49°). The solid lines represent the radial profiles for models at different distances, which are color-coded. The green horizontal line above the profiles represents the beam size.
In the text

	Fig. F.3 Comparison of the emission size measured in the image with the theoretical size. The red area indicates where the angular resolution given by the (u, v)-coverage is not sufficient in resolving the emission. The dashed line represents the 1:1 relation between the theoretical size and the size determined from the image.
In the text

	Fig. F.4 Top: image reconstruction of a sharp model with an inclination of 49° degrees and located at 100 pc. The blue star represents the position of the star. The contours are 1 (dotted), 2, and 3-sigma pixel significance levels. Bottom: corresponding asymmetry map with the contour representing the 2σ significance level.
In the text

	Fig. F.5 Squared visibilities of the dataset and the images. For each target, the top panel shows the squared visibilities for the data (black circles) and the image (blue triangles).Bottom panel: residuals normalized by the error bars on each data point.
In the text

	Fig. F.6 Same as Fig. F.5.
In the text

	Fig. F.7 Closure phases of the dataset and the images. For each target, the top panel shows the closure phases for the data (black circles) and the image (blue triangles). Bottom panel: residuals normalized by the error bars on each data point.
In the text

	Fig. F.8 Same as Fig. F.7.
In the text