© ESO 2021

1 Introduction

It is widely accepted that stars of practically all masses are initially surrounded by accretion disks that play a fundamental role during their formation (e.g., Beltrán & de Wit 2016, and references therein). Thus, the understanding of early stellar evolution is directly related to our knowledge of the accretion mechanism transporting material from the disk to the star (see, e.g., the review in Mendigutía 2020). In Classical T Tauri stars (CTTs, 0.1 < M_*∕M_⊙ < 2), the consensus is that accretion occurs through the stellar magnetic field that connects the inner gas disk and the central star, according to the magnetospheric accretion scenario (MA, Uchida & Shibata 1985; Koenigl 1991; Shu et al. 1994). Magnetic fields in massive young stellar objects (MYSOs, M_*∕M_⊙ > 10) are, in principle, much smaller (e.g., Briquet 2015, and references therein), for which their gaseous disks could reach the stellar surface and directly accrete through a hot and dense boundary layer (BL, Lynden-Bell & Pringle 1974). The possible transition between MA and BL probably occurs within the Herbig Ae/Be (HAeBe, mostly with 2 < M_*∕M_⊙ < 10) regime, as suggestedby different lines of evidence (see, e.g., Ababakr et al. 2017; Wichittanakom et al. 2020; Mendigutía 2020, and references therein). In particular, although the flux excess in the near-ultraviolet (nUV) Balmer region of the spectra (the “Balmer excess” ΔD_B) of CTTs and most HAeBes is interpreted from gas accreting magnetospherically, several early-type HBes show such large ΔD_B values that could not have been be reproduced from MA shock modelling, and which may accrete through a BL (Mendigutía et al. 2011; Fairlamb et al. 2015).

Overall, HAeBes are the most massive stars that still show an optically visible pre-main sequence evolution, given that MYSOs remain embedded in their natal envelopes until they reach the main sequence. Moreover, HAeBes are relatively bright and close, which also make them ideal targets for optical and near infrared (nIR) interferometry (see, e.g., the review in Kraus 2015, and references therein). Such an observational technique is unique in that it gathers spatially resolved information at a very close radial distance from the stellar surface (~ 1 au), which is relevant for understanding the star-disk interaction and the accretion process.

Interferometric studies have revealed that for most HAeBes, there is a correlation between the dust inner sizes probed through the spatially resolved nIR continuum emission and the stellar luminosity (Monnier & Millan-Gabet 2002; Monnier et al. 2005). This empirical correlation takes the form of $r_{\text{in}}\propto L_{*}^{1/2}$, with significant scatter, which is consistent with optically thin inner disks and dust destruction radii of relatively large grains (~ 1 μm), with sublimation temperatures in the range 1000–2000 K. Intriguingly, some high-luminosity, early-type HBes were found to fall below their expected positions in the size–luminosity diagram, given that their inner dust sizes are smaller than inferred from their stellar luminosities. As discussed, for instance, in the reviews by Millan-Gabet et al. (2007) and Kraus (2015), the most plausible reason to explain the undersized nIR continuum emission of such HBes is the presence of gas very close to the star. This gas would shield the dust from stellar irradiation, allowing it to survive closer to the star and to emit in the nIR on top of the dust emission, making the sources appear more compact. Indeed, the undersized HBes show inner dust sizes that are more consistent with classical models having optically thick gas that reaches the stellar surface (Millan-Gabet et al. 2007; Kraus 2015, and references therein). This is in line with the BL accretion scenario, rather than that of the MA.

Studies based on spectro-interferometry, which are capable of spatially resolving optical and nIR emission lines such as Hα and Brγ, reveal that the spatial distribution of the hot atomic gas is more compact or more extended than the dust sublimation front (depending on the specific HAeBe), which has been associated to accretion and wind processes (e.g., Kraus et al. 2008; Mendigutía et al. 2017, and references therein). In fact, although there are works devoted to specific stars that suggest that their comparatively smaller inner dust radii could be explained by the presence of inner gas emitting in the nIR (e.g., MWC 147 in Hone et al. 2019), no general trend has been reported.

Alternative explanations for the different location of HAeBes in the size–luminosity diagram, based on the dust properties, have also been proposed. Benisty et al. (2010) tentatively excluded gas as the origin of the nIR continuum emission, at least for HD 163296, proposing instead very refractory grains as the only cause. However, they also noted that such grains should somehow survive temperatures well above 2000 K, inside the conventional dust sublimation zone. More recently, GRAVITY Collaboration (2019) proposed that the SED shape according to the Meeus et al. (2001) classification scheme could be related to the different inner dust sizes from K-band interferometry. According to the Meeus et al. (2001) classification Group I stars show a rising continuum from the IR to the sub-millimeter region that can be fitted by a power-law component plus a cool black body, whereas Group II sources only need the power-law component. Studies based on high-resolution imaging associate Group I sources with the presence of “transitional” disks with large gaps and cavities (Maaskant et al. 2013; Honda et al. 2015; Garufi et al. 2017). This possible link between Group I sources and transitional disks led GRAVITY Collaboration (2019) to propose the hypothesis that the Group I characteristic sizes in the K-band are generally larger than those in Group II for a given luminosity bin, and a possible trend may be present in their data in the 10–100 L_⊙ range. Although no definitive answer could be provided in that work, the confirmation of the hypothesis by GRAVITY Collaboration (2019) may give an important, alternative clue on the divergent location of HAeBes in the size–luminosity diagram. In fact, the classification of SEDs based on the Meeus et al. (2001) scheme is unrelated to the presence of inner gas based on accretion rates (Mendigutía et al. 2012; Guzmán-Díaz et al. 2021, GD21 hereafter).

In this work, we incorporate GRAVITY/VLTI (GRAVITY Collaboration 2017) K-band interferometric data of five HBes in the analysis of the size–luminosity correlation. The particular properties of these sources, described in following sections, may be helpful for understanding the origin of such a correlation and its possible relation with the presence of inner gas, the accretion mechanism, and the SED properties. Section 2 describes the observations and data reduction. Section 3 presents the analysis of this dataset. The re-assessment and interpretation of the size–luminosity correlation is presented in Sect. 4. Our main conclusions are summarized in Sect. 5.

2 Sample, observations, and data reduction

Our sample of stars was chosen from Fairlamb et al. (2015) to represent two extreme cases within the HBe regime: sources with negligible amounts of accreting gas interior to the dust disk in accordance with their null values of ΔD_B (V590 Mon, PDS 281, and HD 94509), and those with large amounts of hot gas very close to the star based on the presence of strong ΔD_B values that cannot be modeled with MA (DG Cir and HD 141926). In addition, the sample includes the two types of SEDs according to the classification in Meeus et al. (2001) and spans a wide-enough range of stellar luminosities (1 < log L_*∕L_⊙ < 5) to essentially cover the whole size–luminosity diagram for HAeBes. Additional specific properties of the stars studied in this work are described in Sect. 3.

Each star was observed twice during one night between the end of 2018 and the beginning of 2019 with GRAVITY in the wavelength region of ~ 2.0–2.4 μm and high spectral resolution (~4000). The four 8.2 m Unit Telescopes (UTs) were used for all stars, except for HD 141926, which is bright enough to use the four 1.8 m Auxiliary Telescopes (ATs). The single-field mode was used, meaning that the interferometric signals of the targets were recorded on the fringe tracker (FT) and the science channel (SC) detectors. Spatially unresolved calibrators with similar K-magnitudes and spectral types were also observed with the same instrumental configuration and typical CAL-SCI sequences. The seeing was between 0.4” and 0.6” in all cases. Table A.1 in Appendix A shows the observing log, and Figs. A.1 to A.5 show the uv coverage of the observed sources.

The data were reduced with the GRAVITY data reduction pipeline (Lapeyrere et al. 2014, version Esoreflex-2.9.1). The same standard procedures were applied both to the targets and calibrators to carry out dark, flat-field, and bad-pixel correction, as well as toextract the interferometric observables. The squared visibilities of the targets were derived by dividing those reduced observables by the unresolved calibrators, whose spectra were also used to remove telluric lines and for the final wavelength calibration. Figures A.6 to A.10 show the interferometric observables around the Brγ line for each target star and baseline: fluxes, squared visibilities, and differential phases.

Finally, Fig. A.11 shows the closure phases measured for each star and telescope triplet. The average closure phase for all sources is ~ 0° regardless of the triplet, with a scatter around that value < 5°. The only possible exception is V590 Mon, whose average closure phase is ~ 20°, although the noise has a comparable size. Additional details inferred from the previous observational results are described in the next section.

3 Interferometric results

Although simple Keplerian disks are capable of reproducing Brγ spectro-interferometric observations in some HAeBes (Kraus et al. 2012a; Ellerbroek et al. 2015; Mendigutía et al. 2015a), hydrogen emission lines are usually interpreted in terms of magnetically driven accretion or winds (e.g., Muzerolle et al. 2001; Kurosawa & Romanova 2013; Tambovtseva et al. 2014, 2016, 2020). However, our sample includes two stars whose strong nUV excess cannot be modelled from MA, and there is a lack of alternative BL models of line emission (Mendigutía 2020, and references therein). Therefore, in the following we will focus on deriving the sizes associated to the nIR continuum adjacent to Brγ, only providing qualitative information on this spectral line. The specifics for each star in the sample are discussed in the next subsections.

3.1 V590 Mon

V590 Mon (aka LkHa 25, Walker 90) is a HBe star with a stellar temperature of ~ 12 500 K and mass of ~ 3 M_⊙ (Fairlamb et al. 2015; Moura et al. 2020). A distance of 689${}_{-49}^{+57}$ pc is inferred from the recent Gaia EDR3 parallax (Lindegren et al. 2020). After applying the zero-correction described in that work, such a distance has been derived simply by inverting the parallax, given that the relative error is < 0.1 (Bailer-Jones et al. 2021). No excess in the nUV, Balmer region of the spectrum (ΔD_B = 0 magnitudes)was reported by Fairlamb et al. (2015), for which V590 Mon can be considered as virtually non-accreting based on this direct accretion probe. Although V590 Mon has been considered an accreting source in the past, this could be associated with the significant variability of this star on timescales of decades (Pérez et al. 2008; Moura et al. 2020, and references therein). Indeed, the more recent work by Moura et al. (2020) shows that the Hα emission currently shown by V590 Mon is associated to outflowing material in a disk wind and that accretion has a negligible contribution to its Hα line profile. A companion at > 5′′ from the central object and fainter by 6.6 magnitudes in the K-band was reported by Thomas et al. (2007). At such a relatively large separation, the binary nature of V590 Mon could not be confirmed by Wheelwright et al. (2010) based solely on optical spectro-astrometry. The SED of V590 Mon can be classified as a Group I according to the Meeus et al. (2001) scheme, based on the rising slope at IR wavelengths (Fig. 2 in Moura et al. 2020).

Our GRAVITY data reveal that the single-peaked Brγ emission and the adjacent continuum are resolved at all baselines (Fig. A.6). However, there is no difference between the visibilities of the line and the continuum, for which the corresponding emitting regions have comparable sizes. The continuum visibilities in the K-band adjacent to the Brγ line observed at different baselines were modeled using the LITpro software tool (Tallon-Bosc et al. 2008). In particular, we tested parametric geometrical models consisting of a central point source plus an elongated dust distribution with a common center and the shape of a ring, a uniform disk, and a 2D Gaussian disk, given that these are the most common structures normally used to interpret interferometric data of HAeBes. The K-band excess estimate in Vioque et al. (2018), based on the same stellar parameters from Fairlamb et al. (2015) that we adopt here, leads to weight contributions of the central star and the dust distributions of ~ 10% and ~ 90%, which were taken as an initial guess for the fitting procedures. The final weights, as well as the position angles, inclinations, and inner dust sizes resulting from the modelling are shown in Table 1. Reduced values of χ² normalized by the number of variables in each model (i.e., ${\chi }_{\text{r}}^2$), are used to quantify the goodness of the fits and are also tabulated.

The Gaussian disk provides the best fit in terms of a smaller ${\chi }_{\text{r}}^2$ value, followed by the ring and the uniform disk. However, all of them result in the same point source contribution within errorbars, larger than the stellar contribution suggested by the K-band excess reported in Vioque et al. (2018). On top of a possible contribution from unresolved circumstellar emission, this difference may be ascribed to nIR variability or to a stellar temperature different from the one reported in that work (see Sects. 3.2, 3.3, 3.4, and 3.5). Both the Gaussian disk and the ring models provide an intermediate inclination angle consistent with that resulting from the recent modelling of the Hα emission of V590 Mon (40°–55°; Moura et al. 2020). The inclination and major axis position angle (PA, measured from North to East) resulting from the uniform disk model have associated errorbars larger than the own values and, as such, the values are unreliable. The PAs inferred from the Gaussian disk and the ring differ significantly, although the dependence on the initial values assumed is negligible inthe Gaussian case, but strong in the ring case. The limited uv coverage in our observations prevents us from establishing firm conclusions concerning the PA (see Fig. A.1 and Sect. 3.4). In therest of this work, the inner dust size inferred from the Gaussian disk model is assumed, leading to 2.99 ± 0.02 mas, or 2.1 ± 0.2 au at the distance to the star. Figure 1 shows graphically that such a model is the one that better reproduces the squared visibility data of the nIR continuum adjacent to the Brγ emission, as expected from its comparatively smaller ${\chi }_{\text{r}}^2$ value.

All previously discussed models are centro-symmetric and, as such, they can only reproduce a closure phase close to 0°. However, V590 Mon shows the largest errorbar associated with the closure phase, which can range from ~ 0° to ~ 20° (Fig. A.11 and Sect. 2). Such a large errorbar also makes difficult to test whether there is a possible variation of the closure phase with the spatial resolution provided by different baselines, a variation expected if there is a real asymmetry of the source. Still, we tested the eventual presence of a binary that could reproduce the possible closure phase signal of ~ 20°. This companion should be different from the above mentioned in Thomas et al. (2007), given that the wide separation reported in that work is orders of magnitude larger than the GRAVITY/VLTI field of view (limited to the Airy disk of each individual aperture, i.e., 60 mas for the UTs in K-band). Although different binary configurations (with and without associated circumstellar structures) are capable of reproducing closure phases close to 20°, they significantly fail to reproduce the distribution of visibilities versus spatial frequencies shown in Fig. 1. Another possibility for breaking the symmetry is to use models with a single star and an associated circumstellar structure not centered at the same position. However, in order to reproduce a closure phase ~ 20° models require that the center of such circumstellar structure falls well outside the GRAVITY field of view, which is physically unrealistic. According to GRAVITY Collaboration (2019), closure phases ~ 25° or lower can be reproduced from a combination of inclination effects and different dust grain distributions that introduce some azimuthal modulation in the circumstellar profile. In the case of a confirmed closure phase close to 20° in future observations of V590 Mon, this may be a reasonable explanation that is consistent with the anomalous circumstellar dust extinction inferred for this star from nUV spectra (Sitko et al. 1984).

Fig. 1 Squared visibilities measured for the continuum adjacent to Brγ in V590 Mon as a function of the spatial frequency (solid dots with errorbars). The best fit is indicated with a solid line and corresponds to a Gaussian disk on top of the central star. Two additional fits corresponding to ring and uniform disk models are plotted with dashed and dot-dashed lines.

3.2 PDS 281

PDS 281 (aka SAO 220669) is a HBe star with a stellar temperature of ~ 16 000 K (Fairlamb et al. 2015; Wichittanakom et al. 2020). Its mass is ~ 9 M_⊙ and it is located at a Gaia EDR3 distance of 868 ${}_{-9}^{+9}$ pc (GD21). Based on the Hα line, Vieira et al. (2003) classified PDS 281 as a relatively evolved HAeBe star. Indeed, the Hα emission resulting after the subtraction of the photospheric contribution is weak, the Brγ line appears fully in absorption, and ΔD_B = 0 magnitudes, indicating negligible disk-to-star accretion (Fairlamb et al. 2015, 2017). Still, the forbidden [O I] line at 6300 Å is in emission (Fairlamb et al. 2017), possibly indicating the presence of winds. Vioque et al. (2018) reported an IR excess where the circumstellar contribution represents around 30% of the observed flux at the K-band, with such a value based on stellar parameters compiled from Fairlamb et al. (2015). The more recent estimates of the stellar parameters based on Gaia distances in Wichittanakom et al. (2020) leads to a null IR excess up to ~ 10 μm according to GD21. Such a “transitional” SED would be more typical of a relatively evolved young star. GD21 also classified the SED of PDS 281 as Group I.

Figure A.7 shows that the Brγ absorption line and the adjacent continuum for PDS 281 are unresolved in our GRAVITY visibilities. In addition, the differential and closure phases are 0° for all baselines (see also Fig. A.11), indicating a common photocentre and a symmetric source for both the line and the continuum. The fact that there is no difference between the visibilities of the absorption line and the continuum suggests that the latter may also be photospheric, in agreement with the null excess at the K-band reported in GD21. An upper limit to the corresponding diameter is given by the best spatial resolution achieved, that is, λ/2B, with λ ~ 2.17 μm and B the largest projected baseline for PDS 281 during the observations (130.16 m). The resulting upper limit radius is 0.86 mas (half the angular resolution), or 0.74 au (159 R_⊙) at the distance to the source. Such an upper limit is well above the stellar radius of PDS 281 estimated by Wichittanakom et al. (2020) and GD21 (~ 10 R_⊙), which cannot be resolved with the VLTI. Assuming that the K-band continuum emission is actually photospheric, the size of the inner dust emitting region is in principle much larger than 0.74 au and could be determined observationally by using interferometry at ~ 10 μm, the shortest wavelength where the IR excess is apparent according to GD21.

3.3 HD 94509

HD 94509 has been commonly catalogued as a HAeBe star (e.g., Vioque et al. 2020), although it has also been considered an evolved B-star with emission lines (e.g., Cowley et al. 2015, and references therein). The stellar mass is ~ 6 M_⊙, and it is located at a Gaia EDR3 distance of 1609${}_{-43}^{+45}$ pc (GD21). The nUV excess in the Balmer region of the spectrum was reported to be null in Fairlamb et al. (2015), indicating negligible amounts of gas accreting onto the central star. Still, Fairlamb et al. (2017) report relatively strong emission in Hα, Brγ, and [O I]6300. The nIR SED of HD 94509 was considered to be in excess with respect to the underlying photosphere(Vioque et al. 2018) based on a stellar characterization compiled from Fairlamb et al. (2015), but again the recent estimates in Wichittanakom et al. (2020) and GD21 lead to a null excess up to ~ 3 μm. The decreasing IR slope in the corresponding SED shown in GD21 suggests that HD 94509 may be a Group II star, although photometry at wavelengths longer than 20 μm is needed to confirm this.

Based on our GRAVITY visibilities for this object (Fig. A.8) the continuum emission remains unresolved for all baselines. The closure phase is ~ 0° (Fig. A.11), for which the source can be considered centro-symmetric. As for PDS 281 above, an upper limit for the radius of the continuum emitting region can be estimated from the longest baseline in our observations (130.16 m) resulting in <0.86 mas, or <1.38 au (298 R_⊙) considering the distance to HD 94509. Once again this upper limit is well above the previous determinations of the stellar radius by Wichittanakom et al. (2020) and GD21 (8–9 R_⊙), which couldnot be resolved with the VLTI if the nIR continuum is actually photospheric.

Concerning the double-peaked Brγ emission, it is also unresolved for most baselines, except for some residual visibility signal in baseline U4U3, perhaps indicating an emitting region of the blue-shifted peak larger than that of the continuum. In addition, there are clear signals in the differential phases measured in baselines U4U2, U4U1, and U3U1, (with some residuals in U2U1 and U3U2 too) indicating that the photocentres of the blue and red peaks of the Brγ emission are displaced with respect to each other and with respect to the continuum. A detailed modeling of the Brγ emission in this star is beyond the scope of this work.

3.4 DG Cir

DG Cir (aka VdBH 65b, HBC 596) is a ~ 3.4 M_⊙ HBe star at a Gaia EDR3 distance of 861${}_{-14}^{+15}$ pc (GD21) embedded in a nebulosity (Tisserand et al. 2013). This star is one of the few HAeBes with a Herbig-Haro source associated, presumably caused by substantial mass loss (Ray & Eisloeffel 1994, and references therein). Indeed, accreting gas is also present, based on such a large nUV excess that cannot be modeled from MA (ΔD_B ~ 0.80 magnitudes; Fairlamb et al. 2015). Excess is also apparent at short, nIR wavelengths, with an SED shape belonging to Group I (GD21).

Both the single-peaked Brγ emission and the adjacent continuum are resolved at all baselines, although there is no difference between the visibilities of the line and the continuum (Fig. A.9), indicating similar sizes of the corresponding emitting regions. Null values of the differential and closure phases indicate common photocentres and centro-symmetric distributions, respectively. The modeling of the continuum was carried out with LITpro in a way similar as for V590 Mon above (Sect. 3.1). The K-band excess reported by Vioque et al. (2018) is based on the same stellar temperature adopted here (Wichittanakom et al. 2020) and indicates that the central star contributes with ~ 3% of the total flux at those wavelengths. Such a weight was assumed as an initial guess for the modeling. On top of the ring and Gaussian/uniform disk distributions that have been tested, the fits significantly improve by adding a background, halo contribution, as happens in other HAeBes (see, e.g., Lazareff et al. 2017; GRAVITY Collaboration 2019). Table 2 shows theresults of the different model fits. Although the Gaussian model best fits the observations in terms of a smaller ${\chi }_{\text{r}}^2$ value, all models provide the same values within the errorbars. The final contribution inferred from the GRAVITY data for the point source is significantly larger than what is expected from the K-band excess mentioned above from Vioque et al. (2018). Apart from unresolved circumstellar emission and potential nIR variability, uncertainties on the stellar characterization assumed may explain the difference. In fact, the photometric stellar characterization recently made by GD21 leads to a stellar contribution in the K-band of ~ 20% (see the corresponding SED in GD21), which is in agreement with our estimates from the interferometry.

On the other hand, Ray & Eisloeffel (1994) reported that DG Cir has associated the Herbig-Haro optical flow HH 140. From a visual inspection of the image in that work, knots A, B, C, and D constituting HH 140 extend along the NW-SE direction at a PA ~ 135°–160°. This PA is displaced < 60° from the one inferred here and, thus, the outflow would not be roughly perpendicular to the disk as one would expect. In contrast, the intrinsic polarization shown by this star has a position angle ~ 40° (Rodrigues et al. 2009). If such a polarization is tracing the disk major axis, then this would be closer to the perpendicular direction from the HH flow, as expected. However, our fits worsen if the PA of the major axis is forced to bein the NE direction and the best fit goes to the largest angle possible within that quadrant (90°). The uv coverage in our observations may bias the PA inferred given that we are probing the NW-SE direction (Fig. A.4), and additional observations covering the NE-SW direction would be useful to properly constrain that value. In turn, the intermediate inclination angle inferred for DG Cir, closer to edge-on, is consistent with the presence of P Cygni profiles in optical spectra (Allen 1978; Vieira et al. 2003) and with that source being classified as an UXor-type variable by the AAVSO¹. The uv coverage does not affect the estimated K-band inner dust size, and 0.5 ± 0.2 mas, or 0.4 ± 0.2 au at the distance to the star, is adopted hereafter as a reliable value in agreement with the models fitted to the observations.

Finally, Fig. 2 overplots the best fit model in Table 2 to the measured squared visibilities versus the spatial frequencies. The rest of the models in that table are not overplotted because there is no visual difference between them.

Fig. 2 Squared visibilities measured for the continuum adjacent to Brγ in DG Ciras a function of the spatial frequency (solid dots with errorbars). The best fit is indicated with a solid line and corresponds to a Gaussian disk on top of the central star, and a background contributing ~ 10%. Two additional fits corresponding to ring and uniform disk models are not plotted as they are practically indistinguishable from the previous.

3.5 HD 141926

HD 141926 (aka PDS 399, Hen 3-1110, 1H 1555-552) was initially catalogued as a Classical Be star (e.g., Jaschek & Egret 1982), but it was later identified as a HAeBe source (Vieira et al. 2003) and considered as such in more recent works (see, e.g., the references in this subsection). HD 141926 is a ~ 17 M_⊙ HBe star at a Gaia EDR3 distance of ~ 1324${}_{-31}^{+32}$ pc (GD21). This source has X-ray emission associated (e.g., Maccarone et al. 2014; Liu et al. 2000, and references therein), which adds to the increasing number of HAeBes with such detections (Stelzer et al. 2009, and references therein). Concerning the nUV, Fairlamb et al. (2015) found significant excess over the photosphere (ΔD_B ~ 0.20 magnitudes)that could not be reproduced from MA but may be associated with a large amount of accreting gas close to the central star (for a compendium of the emission lines detected in this source see also Fairlamb et al. 2017). An IR excess is also present from short, nIR wavelengths, with the inner dust disk contribution ranging between 30 and 40% of the total observed flux in the K-band, depending on the stellar parameters adopted (Vioque et al. 2018; GD21). The overall SED is classified as Group II (GD21).

Our GRAVITY visibilities reveal that the emitting regions of both the double-peaked Brγ line and the adjacent continuum are unresolved (Fig. A.10). The lack of significant signals in the differential and closure phases (Fig. A.11) reveal common photocentres and symmetric distributions for both the emission line and the continuum. An upper limit to the corresponding radius is derived as for PDS 281 and HD 94509 above, but in this case under the consideration that the largest baseline in our observations was 32.0 m (the ATs were used instead of the UTs, Sect. 2 and Table A.1). The corresponding dust radius obtained is <3.49 mas, or <4.62 au at the distance to the star.

4 Analysis of the size–luminosity correlation

In order to put our previous results in context with the size–luminosity relation, Table 3 includes relevant data for our stars and for all HAeBes with previous K-band interferometric sizes from the literature (see the references in the caption of that table). Spatially resolved stars discarded as HAeBes in Vioque et al. (2020) have been excluded from the list. The angular inner dust sizes listed in the table were reported in the corresponding references based on model fits of the visibility vs spatial frequency curves, as in our analysis in the previous section. Such angular sizes have been transformed here into spatial sizes using the distances inferred by GD21 based on recent Gaia EDR3 parallaxes (Lindegren et al. 2020, see also Sect. 3.1 for V590 Mon). Following GD21, a few distances in the table have been tagged when there is some potential degree of spuriousness, based on uncertainties and flags associated to the new Gaia EDR3 parallaxes. However, as discussed in GD21 the criterion adopted in this work is conservative and the probability of having a correct Gaia EDR3 astrometric solution according to Rybizki et al. (2021) is larger than 0.9 for all stars except for Z CMa A. All distances will be used here, but the reader should be aware of these caveats. Stellar temperatures of most stars were obtained in Wichittanakom et al. (2020) from observed spectra. Adopting the relation with the spectral type in Kenyon & Hartmann (1995), 9520 K is the effective temperature dividing the sample in hotter HBes and colder HAes. Four stars show temperatures below 7200 K and thus spectral types later than A9. For simplicity, they have been considered as HAes based on previous catalogues (see, e.g., the discussion in GD21) and their relatively large stellar masses (≥1.5 M_⊙). The rest of the stellar parameters (L_*, R_*, and M_*) were mostly derived by GD21 taking those temperatures and the Gaia EDR3 distances as a departure point. The associated errors are mainly based on SED fitting uncertainties, as explained in that work. It is noted that the errors tabulated in Table 3 are considered symmetric for simplicity, without affecting the analysis in this work. Concerning the distances, almost all sources have small associated errors of the parallaxes (<1–5%), for which the errors are indeed symmetric for them. We encourage the reader to consult GD21 for more specifics on the distance, the rest of the stellar parameters, and their associated errors. The final samplecomprises 48 HAeBes spanning almost 5 orders of magnitude in stellar luminosities (0.5< log L_*∕L_⊙ < 5.2). As such, the following analysis of the K-band size–luminosity relation will be based on the largest sample of HAeBes to date, almost doubling previous samples used for similar analyses (GRAVITY Collaboration 2019).

Additional information for the same stars concerning the nUV Balmer excess, the Hα and accretion luminosities, the mass accretion rate, the presence of atomic gas inside the dust radius inferred from spectro-interferometry, and the SED group according to the Meeus et al. (2001) classification scheme, is listed in Table 4 when possible. This information is further discussed in the following sections.

Figure 3 shows the size–luminosity plot from the data in Table 3. The expected inner dust disk sizes corresponding to two different scenarios are indicated with solid and dashed lines. The solid lines refer to the model where the inner disk is optically thin. The stellar luminosity establishes the inner dust radius, located at the region where the sublimation temperature (T_d) is reached. In particular, it is assumed that r_in (au) = 1.1 × ${(Q_{\text{R}}{L}_{*}/1000L_{\odot })}^{1/2}\times {(T_{\text{d}}/1500K)}^{-2}$ (Tuthill et al. 2001; Monnier & Millan-Gabet 2002; Monnier et al. 2005), where Q_R = Q_abs(T_*)∕Q_abs(T_d) is the ratio of the dust absorption efficiencies Q(T) for radiation at color temperature T of the incident and reemitted field, respectively. The rough boundaries limiting the expected inner dust disk from this model have been overplotted in Fig. 3 with solid lines, obtained assuming values for T_d and Q_R of 1000 K and 10 (upper limit), 2000 K and 1 (lower limit). On the other hand, the classical model assumes an optically thick gas disk reaching the central star and located in-between the stellar surface and the inner dust (Hillenbrand et al. 1992). Such a gas disk partially shields the stellar radiation, allowing the dust to survive closer to the star, hence shrinking the value of r_in. In this casethe disk temperature T(r) is assumed to be proportional to r^−3∕4. The proportionality constant depends on the specific stars and is derived from their corresponding R_* and T_* values, thus for each object $r_{\text{in}}=R_{*}{(T_{*}/T_{\text{d}})}^{4/3}$. Finally, the log r_in values from the previous expression linearly correlate with the corresponding log L_* values for the stars in the sample, and the bestlinear fits for T_d = 1000 K and 2000 K are overplotted in Fig. 3 with dashed lines. These are the boundaries where the dust destruction radius based on the classical model must lie.

It has been remarked that in the classical model, the gas reaches the stellar surface, as in the BL scenario of accretion. Thus, we will use the term “optically thick-BL model” hereafter. In contrast, the above described optically thin model is commonly associated to MA, for which we will refer to it as the “optically thin-MA model”. We note, however, that the link in this case is not straightforward, given that in the MA scenario the inner gas may also be optically thick for sources with strong accretion rates (Muzerolle et al. 2004).

According to Fig. 3 (top left), all HAes fall within the region predicted by the optically thin-MA model in the size–luminosity diagram, although the location of 5 out of 23 of such sources (~ 22%) could also be consistent with the optically thick-BL scenario considering errorbars. In contrast, the distribution of HBes is more heterogeneous. Considering the errorbars, as well as the upper and lower limits, up to 18 and 13 out of 25 of such sources (<72% and <52%) fall within the regions predicted by the optically thin-MA and optically thick-BL models, respectively. This difference between theHAes and many HBes is recovered and expanded here from the initial studies of the size–luminosity relation (Monnier & Millan-Gabet 2002; Monnier et al. 2005), as it has been somewhat lost in more recent surveys (Lazareff et al. 2017; GRAVITY Collaboration 2019). The reason is that these surveys do not include HBes with L_* > 10⁴ L_⊙, explaining why the previously mentioned difference could not be clearly observed.

In the following, we test the hypothesis that the presence or absence of gas in between the dust and the star or, otherwise, the different accretion scenarios might explain the locations of HAeBes in the size–luminosity diagram; in particular, the observed scatter ≥ ± 0.5 dex in r_in for a given L_*. Other possibilities related with the dust distribution based on the SED shape and with variability is also considered. Finally, the potential influence of the distances to the sources on the size–luminosity relation is analyzed.

Fig. 3 Size-luminosity relation based on data compiled in Tables 3 and 4 for HAeBes. Each panel emphasizes a different aspect, as indicated in the legends. Top left: HAes and HBes; the five with new GRAVITY data presented in this paper are highlighted in red. Top right: the HAeBes with Balmer excess measurements above and below the median value (0.07 magnitudes). Bottom left: the HAeBes with spectro-interferometric measurements indicating that the Brγ/Hα emitting region is smaller than the dust emitting region and the ones without such evidence. Bottom right: Group I and Group II sources based on the SED shape. (All panels): triangles indicate upper and lower limits, and the solid and dashed lines indicate the expected boundaries of rin for an optically thin-MA and optically thick-BL inner gaseous disks, respectively (see text for details).

4.1 Inner dust size and the presence of gas

As mentioned before, the main hypothesis for explaining the comparatively smaller r_in values shown by some HAeBes (HBes in particular) with respect to others in the size–luminosity diagram is the presence of large amounts of gas in-between the dusty region and the central star, whose shielding effect on the stellar radiation would allow the dust to survive closer to the central source. We note that although predictions from two different models reflecting optically thin and thick inner gaseous disks have been overplotted in Fig. 3, a smooth transition between the former and the latter scenarios is expected as the mass accretion rate onto the central stars increases (and thus, so does the inner gas density: see Muzerolle et al. 2004, for details). In turn, optically thick gaseous disks reaching the central star as in the models by Hillenbrand et al. (1992) are also associated with the BL mechanism of accretion, instead of MA. Different gas probes are used next to test the different aspects of the previous hypothesis.

4.1.1 Inner gas probed through the Balmer excess

The main reason for using GRAVITY to observe the five HBes analyzed in Sect. 3 is to test the above-mentioned hypothesis based on two types of HBes that it predicts would be located in very different regions of the size–luminosity diagram. On the one hand, HBes with undetectable Balmer excesses indicating the lack of measurable accretion and thus a negligible amount of gas very close to the central star (V590 Mon, PDS 281, and HD 94509) should have comparatively large r_in values. On the other, HBes with large amounts of gas shocking onto the central source as inferred from their large Balmer excesses that cannot be modeled from MA (DG Cir and HD 141926) should have comparatively small r_in values.

The top left panel of Fig. 3 shows the location of the five HBes analyzed in this work in the size–luminosity diagram. Although V590 Mon and HD 94509 show r_in values larger than the rest of the stars in similar stellar luminosity bins, as expected from the negligible amount of inner gas inferred from their null values of ΔD_B, PDS 281 has a comparatively undersized inner dust size, contrary to the expectations. However, it must be remarked that the K-band sizes derived for PDS 281 and HD 94509 may refer to upper limits of the stellar photospheres and not to the inner dust disks, which could be considerably larger than plotted if the nIR excesses actually start at wavelengths longer than ~ 2 μm (see the corresponding discussions in Sects. 3.2 and 3.3). On the other hand, DG Cir and HD 141926 do not show significantly smaller r_in values compared to HAeBes with similar stellar luminosities, despite the fact that they do not accrete magnetospherically, and the presence of inner gas very close to the star should be large based on their ΔD_B values. These exceptions indicate that the presence or absence of accreting gas closer to the star than the dust or the different accretion mechanisms cannot be the only parameters that determine the distribution of HAeBes in the size–luminosity diagram. The previous analysis is extended in the top right panel of Fig. 3 by considering all stars in Table 4 with previous measurements of ΔD_B (col. 2). It must be remarked here that ΔD_B is a distance-independent observational quantity from which accretion rates can be directly inferred through appropriate modelling but, regardless of the type of modeling, for a given stellar luminosity bin, large values of ΔD_B may be immediately associated with large accretion rates and the other way around (see Mendigutía 2020, and references therein). However, a general trend where the stars with large excesses (indicated in blue in that panel) above the median ΔD_B in our sample (0.07 magnitudes) have comparatively smaller r_in values is notobserved. The same absence of a trend is obtained when only the ten most extreme cases with the largest and the smallest values of ΔD_B are considered. In addition, apart from DG Cir and HD 141926, VV Ser is the only HAeBe in Table 4 whose large excess cannot be modelled from MA (Mendigutía et al. 2011). However, VV Ser does not show a comparatively undersized inner dust size either. We note that DG Cir, HD 141926, and VV Ser represent ~ 25% of the HBes in Table 4, with measurements of ΔD_B, a percentage similar to the general fraction of HBes that cannot be fitted from MA (based on the stars identified in Mendigutía et al. 2011 and Fairlamb et al. 2015; see Mendigutía 2020). In contrast, the fraction of “undersized” HBes in the size–luminosity diagram is considerably larger than 25%, as described at the beginning of Sect. 4.

In summary, the hypothesis that the presence of a large amount of inner gas or the accretion scenario determines the location of HAeBes in the size–luminosity diagram cannot be confirmed based on measurements of the nUV Balmer excess.

4.1.2 Inner gas probed through the Hα line and accretion

Alternatively, the presence of large amounts of gas relatively close to the central star can be inferred from the luminosity of the Hα emission line, an observational quantity that is available for almost all stars with K-band continuum interferometry in Col. 3 of Table 4. The Hα luminosities come from Wichittanakom et al. (2020), which have been updated using the Gaia EDR3 distances in Table 3. In order to make a comparison between log L_Hα and log r_in we must keep in mind the fact that the stellar luminosity correlates with the luminosity of all nUV-optical-nIR emission lines including Hα (Mendigutía et al. 2015b), namely, more luminous stars – with larger r_in values – also show stronger Hα luminosities. In particular, based on the data in Tables 3 and 4, the stars in the sample follow a relation L_Hα ~ L* (linear correlation coefficient ρ = 0.9). Thus, an appropriate comparison should first remove the previous dependence between L_* and L_Hα dividing the latter by the former. The left panel of Fig. 4 compares the inner dust sizes with the ratios between the Hα and the stellar luminosities. There is no correlation between the presence of gas directly inferred from the Hα luminosities and the inner dust sizes, once L_Hα is normalized to L_*.

The previous distribution is essentially the same if we use model-dependent accretion luminosities, L_acc, instead of Hα luminosities. The reason is that log L_acc and log L_Hα are related by a simple linear transformation (Mendigutía et al. 2011). In particular, accretion luminosities in Table 4 have been derived from the empirical correlation with the Hα luminosity in Fairlamb et al. (2017)², log(L_acc∕L_⊙) = 2.09(± 0.06) + 1.00(±0.05) × log(L_Hα∕L_⊙). Similarly, given that the mass accretion rate is derived from L_acc through the stellar mass and radius (Ṁ_acc = L_acc R_*/GM_*), when Ṁ_acc values are inferred either from MA or from BL they also tend to be larger for stars with larger r_in values. However, this trend is based on the fact that more massive, that is, more luminous, stars are stronger accretors too (Mendigutía et al. 2015b; Wichittanakom et al. 2020, and references therein). From the data in Tables 3 and 4, Ṁ_acc ~ L_* (ρ = 0.9). Thus, once again, the dependence with L_* has to be removed before comparing accretion luminosities and rates with the dust inner sizes. The middle and right panels of Fig. 4 show that once normalized to L_* there is no correlation (ρ < 0.5) between the accretion luminosities or mass accretion rates and r_in, in agreement with our analysis above based on model-independent Balmer excesses and Hα emission lines.

In short, the previous tests do not support that strong accretors generally have comparatively smaller r_in values than weak accretors considering similar stellar luminosity ranges, despite the fact that the amount of inner gas and their densities arelarger for stronger accretors (Muzerolle et al. 2004).

Fig. 4 Inner dust size versus the Hα luminosity, the accretion luminosity, and the mass accretion rate normalized to the stellar luminosity. Linear correlation coefficients (ρ) in each panel quantify the lack of correlation. Upper and lower limits of rin are indicated with triangles.

4.1.3 Inner gas probed through spectro-interferometry

A direct probe of the presence of gas inside the dust emitting region, regardless of whether it is accreting or not, comes from spectro-interferometry. Column 6 in Table 4 indicates the comparatively small sub-sample of HAeBes for which there are spectro-interferometric observations in Brγ and/or Hα from the literature. The bottom left panel of Fig. 3 indicates the location of this sub-sample in the size–luminosity diagram when that is classified in two types of sources. First, the ones for which there is interferometric evidence indicating that atomic gas is emitted in a circumstellar region more compact than that of the dust (blue symbols). Secondly, the rest of the sources for which there is no such evidence (red symbols). Although for the sources in this second group the region where gas is emitted could be equally or more extended than the dust emitting region, the exact situation is unknown for unresolved sources like PDS 281, HD 94509, and HD 141926. There is no general trend and stars with compact emission have both large and small values of r_in, and the same happens to the sources in the second group. However, the sample of HAeBes with spectro-interferometric measurements is still too scarce to infer general conclusions, with less than 38% of the wholesample in Table 4 having such measurements.

Additional spectro-interferometric measurements of larger samples of HAeBes, for instance, using the GRAVITY/VLTI and VEGA/CHARA capabilities to spatially and spectrally resolve the Brγ and Hα lines, are required. These will be of great help to provide an alternative observational test to the hypothesis of the presence of inner gas as a general cause explaining the different locations of HAeBes in the size–luminosity diagram. We note, however, that the objects for which there is interferometric evidence of a gas emitting region more extended than the inner dust (classified above in the second group) may still have gas inside the dust emitting region too. Thus, a careful modeling of the spectro-interferometric observations is necessary to put constraints to the location of the bulk of the gas emitting region relative to the dust.

Finally, the presence of dense gas can also be probed through molecular transitions. In particular, the CO ro-vibrational emission in the wavelength range ~ 2.25–2.45 μm observed in a few HAeBes is probably coming from the innermost disk regions, based on spectral modelling (see, e.g., Carr 1989; Bik & Thi 2004; Ilee et al. 2014). However, none of the five stars in our sample show signatures of CO emission in the GRAVITY data, which is in agreement with the low detection rate in HAeBes (Ilee et al. 2014, 2018). In fact, spatially resolved interferometric observations of CO are only available for one HAeBe star to date (see the works on 51 Oph by Tatulli et al. 2008; GRAVITY Collaboration 2021), which prevents us from making a comparative study.

4.2 Innerdust size and the SED properties

As mentioned in the introduction, the shape of the SED can be associated with the distribution of the circumstellar dust around the central star and eventually to the presence of gaps and holes. The bottom-right panel of Fig. 3 compares the distribution of Group I versus Group II stars in the size–luminosity diagram. Such a classification is available foralmost all sources in the sample (last Col. in Table 4). Our more complete dataset does not support the hypothesis by GRAVITY Collaboration (2019) linking Group I “transitional” stars with larger K-band characteristic sizes at least for the 10–100 L_⊙ luminosity range. Considering such a range, both Group I and Group II scatter along relatively small and large values of r_in values. Moreover, the situation seems to be the opposite to that proposed by GRAVITY Collaboration (2019) for L_* > 10³ L_⊙ (i.e. mainly HBes), where Group I stars appear undersized with respect to Group II sources. However, this conclusion is based on objects that mostly belong to Group I (9 out of 14). In turn, the recent work by GD21 shows that Group I and Group II stars are roughly equally distributed along the HAe and HBe regimes, for which the HBes analyzed here must be biased towards Group I. Additional Group II HBes are thus necessary to test whether there is a distinct location between Group I and Group II stars in the size–luminosity diagram for this regime.

GD21 also classified the SEDs of most HAeBes known to date according to the shortest wavelength where the IR excess is apparent, which is in principle associated to the radial distance from the star where the dust emission originates. They find that while the sources with IR excesses starting at the J or H photometric bands (“Group JH") dominate in the HAe regime, the HBe sources with L_* > 10³ L_⊙ mostly showIR excess starting at the K photometric band or longer wavelengths (“Group K"). In contrast, the sample of HBes with L_* > 10³ L_⊙ in Tables 3 and 4 only includes two Group K stars – PDS 281 and HD 94509, both belonging to our GRAVITY sample. As discussed in Sects. 3.2 and 3.3, the fact that the K-band emission of the latter two stars is unresolved in our GRAVITY observations (as expected for sources with purely photospheric emission at those wavelengths) does not mean that their dust inner sizes are relatively small. On the contrary, stars with IR excesses starting at wavelengths longer than those corresponding to the K-band are expected to have r_in values significantly larger than the sources that show nIR excess and are interferometrically resolved at these wavelengths. As a consequence, a robust analysis of the size–luminosity relation should combine continuum interferometry at the K-band and at longer wavelengths (e.g., with MIDI and MATISSE at VLTI, Menu et al. 2015; Kobus et al. 2019) of a more complete sample of HAeBes belonging to both Group JH and Group K. Such an analysis is likely to add many more HBes with L_* > 10³ L_⊙ and large r_in values. Thus, the apparent trend currently showing undersized HBes is most probably affected by the fact that the current sample includes a comparatively large number of such sources that have an excess starting at short, nIR wavelengths. In other words, a more complete sample of HBes – mostly with IR excesses starting at longer wavelengths – may reveal that the undersized sources are mainly exceptions.

4.3 Variability of the inner dust size

More than one estimate of r_in is available in the literature for 19 objects in Table 3, from which 10 (53%) show differences above errorbars (see Cols. 2 and 4 in that table and the corresponding caption). To establish whether the cause of the different r_in values provided for a given source is actually a temporal change of the inner dust size or is due to other reasons, such as the use of different instrumentation, baselines, and uv coverage, or slight methodological differences, is out of the scope of this work (detailed analysis in, e.g., Kobus et al. 2020, and references therein).

Instead, our aim is to remark that such differences alone cannot explain the observed ≥ ± 0.5 dex scatter in r_in for a given L_*. Indeed, the variations in r_in are smaller by a factor 2 (< ±0.2 dex changes, considering deviations from the average value) in all stars except for HD 179218 and MWC 297. Multi-epoch K-band interferometry of HD 179218 (aka MWC 614) provides r_in values of 3 ± 1 mas with the IOTA instrument (Millan-Gabet et al. 2001), 24 ± 1 mas with AMBER (combined with additional high-resolution instrumentation; Kluska et al. 2018), and 8 ± 1 mas with GRAVITY (GRAVITY Collaboration 2019). Thus, the largest variation is of a factor ~ 8 and, considering deviations from the average value, such changes represent ~ ± 0.4 dex. Concerning MWC 297, this star has been observed multiple times with different interferometers, as indicatedin the references provided in Table 3. The smallest reported value of r_in is 1.68 mas (Eisner et al. 2004). This value is adopted here as a lower limit, considering that r_in can indeed reach more than 12 mas (Millan-Gabet et al. 2001) passing through 2–3 mas (Weigelt et al. 2011; Kraus et al. 2008; Malbet et al. 2007).

In short, apart from a couple of exceptions, so far the data on hand do not support that the observed ~ ± 0.5 dex scatter in the size–luminosity relation can be fully explained from differences in the estimates of the r_in values – regardless of the cause of such differences.

4.4 Considering how distance may affect the size–luminosity correlation

Vinković & Jurkić (2007) showed that the nIR interferometric visibilities of young stars with luminosities below and above 10³ L_⊙ are intrinsically different regardless of the distance, supporting previous studies suggesting that the associated circumstellar sizes differ between these two groups (see Sect. 4 and Monnier & Millan-Gabet 2002; Monnier et al. 2005). However, Vinković & Jurkić (2007) also removed the dependence with the stellar luminosity in their study. Here, we take a critical look at the possible influence of the distance driving the size–luminosity correlation.

Spurious correlations are the ones for which the observed slope, intercept, and scatter are influenced by a common dependence of both parameters compared on a third one. An established way of potentially obtaining a spurious correlation is through arithmetic transformations of both parameters using a common factor (Pearson 1897). Both terms compared in the size–luminosity diagram, L_* and r_in, are arithmetically derived using the distance to the sources, d. In fact, the derivation of absolute luminosities involves a multiplication by a factor d², and the transformation from angular to spatial radii also requires multiplying by a factor d. Therefore, our following analysis on the potential sensitivity of the size–luminosity correlation on the distance is justified.

The partial correlation technique is best suited to deal with this type of statistical issue (e.g., Wall & Jenkins 2003). It measures the degree of association between two variables (L_* and r_in) assessing the potential dependence on a third one (d) that may be the underlying reason causing the correlation. The first row of Table 5 lists the probability of false correlation (p) and the linear correlation coefficient (ρ) quantifying the strength of the size–luminosity correlation for the whole HAeBe sample in Table 3. These quantities have been derived both from standard statistics (left columns, i.e., those that do not consider underlying dependencies) and from partial correlations (right columns, i.e., those after removing the effect of potential common dependencies on the distance). Although the p-value is smaller than the typical threshold to establish a correlation (<0.05) in both cases, the probability of false correlation inferred from partial correlations is orders of magnitude larger than from standard statistics. Similarly, the linear correlation coefficient is above the typical threshold (>0.5) from standard statistics, but not from partial correlations. The previous results suggest that although the relation between log L_* and log r_in is a proper (non-linear) correlation, this is influenced by a common dependence of both parameters compared on the distance to the sources.

The same statistical tests were then carried out considering sub-samples located at narrower ranges of distances. The corresponding ranges and results are indicated in the following rows of Table 5. An additional sub-sample of TTs with nIR interferometric sizes compiled by Pinte et al. (2008) has also been considered (rows 5 and 6 in Table 5). The r_in values listed in their Table 2 have been updated using Gaia EDR3 distances, with all TTs being closer than 200 pc. The comparison between the numbersresulting from standard statistics and from partial correlations indicates that although now the distance does not play a role (the values at the left- and at the right-side columns are similar), there is no statistically significant correlation between log L_* and log r_in for any of these sub-samples (the typical cut-offs for p and ρ are not reached). The only exception is the sub-sample made by combining the TTs and the HAeBes closer than 250 pc. The comparison between the statistical quantities in the last row of Table 5 shows that when these sources are taken together there is a statistically significant correlation that is not affected by underlying dependencies with the distance.

Figure 5 serves to illustrate the case graphically. This is the same as Fig. 3 but the stars have now been color-coded to indicate the distance ranges where they are located, following the same bins as in Table 5. Direct visual inspection supports the conclusions based on the statistical analysis included above: from the current sample, it cannot be concluded that there is a significant correlation between log L_* and log r_in once sub-samples located in narrow ranges of distance are selected. In fact, the stars in each distance bin show roughly horizontal distributions with similar values of r_in within errorbars, despite of the fact that they have stellar luminosities different by orders of magnitude. The only exception is again the sub-sample made of TTs and close HAeBes (green symbols in Fig. 5). These are visually correlated, mainly falling within the predictions of the optically thin-MA model.

The previous discussion poses the question about whether the slope, intercept and scatter of the observed size–luminosity correlation for HAeBes carry any physical significance, or if those are fundamentally driven by the distance. In the following, we argue that even when the size–luminosity correlation is indeed affected by underlying dependencies with the distance, the correlation is most likely to be physically meaningful.

The top and bottom panels of Fig. 6 show the stellar luminosities and inner dust radii of the sampled HAeBes and TTs versus their Gaia EDR3 distances. For d < 250 pc (log d < 2.4), there is a scatter of points both for L_* and r_in. Thus, for nearby sources we are probing a relatively broad range of stellar luminosities (10^−0.5 < L_*∕L_⊙ < 10^2.5) and inner dust radii (0.02 < r_in/au < 1.8) without any dependence on the distance. Precisely because of the lack of such underlying dependencies, the partial correlation statistics confirmed the result from standard statistics: there is a correlation between L_* and r_in for TTs and close HAeBes for which the distance does not play a role. In contrast, from d > 250 pc there is a rising trend linking both L_* and r_in with d. These underlying trends explain why the partial correlations statistics warn that the correlation between L_* and r_in is affected by the distance. We note that although the arithmetic to calculate both L_* and r_in involves the use of d for all sources, it is the presence of underlying correlations between both parameters with d what matters, regardless of the arithmetic. Therefore, in order to assess whether the size–luminosity relation is physically relevant or not, we need first to understand in detail what causes the underlying correlations of both parameters with the distance.

The trend linking L_* and d is common for most HAeBes known to date, as shown by their corresponding parameters from GD21 (overplotted in gray in the top panel of Fig. 6). The reason explaining such a trend is double. Firstly, high-mass, luminous sources are less frequent and evolve faster than lower-mass stars, for which a larger space volume is necessary to detect them. Secondly, low luminosity sources cannot be observed at long distances due to current observational detection limits. This is indicated by the dashed line, above which only the bright sources with an apparent bolometric magnitude <12 are located. On the other hand, the dashed line in the bottom panel of Fig. 6 roughly represents the smallest angular radius that can be resolved with the VLTI at ~ 2 μm, ~ 0.5 mas. It is tempting to associate the lower bound of the trend linking r_in and d to such angular resolution limits. However, if this were the case, we should see mainly upper limits of r_in (i.e., unresolved sources) populating the region near the dashed line, but based on the published data, this is not observed. Instead, the trend linking r_in and d seems to be a by-product of the above relation between L_* and d: high luminosity (distant) sources have intrinsically larger inner disk sizes compared to low-luminosity (nearby) sources. Therefore, even though improving our detection and resolution limits would remove the underlying dependencies with the distance, it would only add additional data for distant, low luminosity sources. But low-luminosity stars are already well probed locally, and thus the currently observed size–luminosity correlation is not expected to change significantly. This conclusion involves two assumptions: Firstly, it is assumed that the r_in values reported for high luminosity stars are representative. If there is a large amount of unpublished upper limits of r_in coinciding with our current resolution limits then the actual size–luminosity diagram differs from the current one. Secondly, it is assumed that the environmentand boundary conditions do not have a critical effect on the inner dust radii. If this was somehow not the case, then the L_* –r_in correlation could eventually change once distant, low-luminosity stars – potentially under different conditions than the ones nearby – are probed.

In summary, the size–luminosity correlation can be statistically confirmed only for TTs and nearby HAeBes, when both sub-samples are taken together, but the overall correlation is affected by underlying dependencies of both L_* and r_in on d. These observed dependencies are caused by a combination of detection limits and the lack of luminous sources in the near vicinity. Improving the detection – and resolution – limits in the future will only add additional low-luminosity sources to the size–luminosity diagram. Given that these are already well probed from nearby stars, it is expected that the whole size–luminosity correlation will remain essentially the same in the future. In this sense, and regardless of the statistical warning, the size–luminosity correlation can in principle be interpreted physically despite the underlying dependencies with the distance. However, we note that the best scenario for studying the size–luminosity relation involves removing the underlying dependencies by considering the widest possible range of stellar luminosities in a narrow enough range of distances. Improving current detection limits with future upgrades such as GRAVITY+³ will facilitate such studies.

Fig. 5 Size-luminosity relation where the stars have been color-coded to indicate different distance ranges as shown in the legend. Interferometrically resolved TTs listed in Pinte et al. (2008) are included in the plot after correction with Gaia EDR3 distances (green open circles), all of them located between 100 and 200 pc from us. The triangles, solid, and dashed lines represent upper and lower limits and possible trends for the dust destruction radius as in Fig. 3.

Fig. 6 Relation with the distance of the stellar luminosity (top) and the inner dust radius (bottom) for the HAeBes (solid circles and triangles to indicate upper and lower limits), and TTs (open circles) with interferometric estimates of rin discussed in the text. The gray symbols in the top panel represent most HAeBes known to date from the data in GD21. The regions below the dashed lines indicate apparent bolometric magnitudes >12 and angular radii <0.5 mas.

5 Conclusions

This work broadens the sample of HAeBes with K-band spectro-interferometric data by adding five HBes observed with GRAVITY/VLTI. Qualitative information has been provided with respect to the Brγ line, but this work has been mainly focused on the analysis of the adjacent continuum. Simple parametrical models have been fitted to the two spatially resolved stars (V590 Mon and DG Cir), providing geometrical estimates of their inner dust distributions. Upper limits for the sizes of the inner dust regions and photospheres are provided for the unresolved stars (PDS 281, HD 94509, and HD 141926).

The previous sample was selected because of the extreme properties with regard to the presence or absence of innermost gas and accreting modes, also covering both Group I and Group II SED shapes. Additional information, including Gaia EDR3 distances and self-consistently obtained stellar and circumstellar properties, has been compiled for these and for all HAeBes with similar interferometric results available in the literature. Based on the most complete dataset concerning K-band interferometric sizes, a reassessment of the size–luminosity relation has been made, leading to the following conclusions:

• We confirm that the location of most HAes is consistent with the predictions of the optically-thin inner disk-MA model inthe size–luminosity diagram, whereas a significant fraction of the more luminous HBes are “undersized” and fall within the range predicted by the classical-BL model with an optically thick gas disk that reaches the central star;

• Based onthe nUV Balmer excesses, the Hα luminosities, the accretion luminosities, or the mass accretion rates, there is no evidence supporting the hypothesis that the shielding effect caused by innermost gas is the general reason explaining the comparatively smaller inner dust sizes of some HAeBes. Underlying relations between the three latter parameters and the stellar luminosity must be removed before any comparison is carried out. The few stars for which direct information on the presence of innermost atomic gas inferred from spectro-interferometry is available do not support the above mentioned hypothesis either;

• Similarly, the more recent suggestion that stars with Group I SEDs have inner dust sizes larger than Group II sources is not supported from our data for any luminosity bin. In fact, the opposite trend is observed for the luminous stars with L_* 1000 L_⊙, althoughthis sub-sample is biased towards Group I stars, and additional luminous Group II sources are necessary to test such a possible trend;

• Based on the stars for which multi-epoch interferometric data is available, the reported variations of theinner dust sizes cannot explain alone the observed scatter in the size–luminosity diagram (≥± 0.5 dex for a given stellar luminosity).

• The size–luminosity correlation is statistically significant only for the sub-sample made of TTs and low-luminosity HAes located at <250 pc, but the overall correlation is influenced by a common dependence of both parameters compared on the distance to the sources. However, although future improvements in our observational capabilities will break the dependencies with the distance, the overall size–luminosity relation will most probably remain essentially the same. Thus, the currently observed size–luminosity correlation can in principle be interpreted physically in spite of the mentioned dependencies with the distance.

Our study is based on an investigation of the main proposed scenarios potentially explaining the different observed positions of HAeBes in the size–luminosity diagram. No associated observational trend capable of confirming any such hypotheses has been found. Thus, a general rationale for the distribution of the stars in the size–luminosity correlation is still lacking.

Acknowledgements

Based on observations collected at the European Southern Observatory under ESO programme 0102.C-0576. The authors acknowledge the anonymous referee for the useful comments, which have served to improve the manuscript. PMA, IM, and JGD acknowledge the Government of Comunidad Autónoma de Madrid (Spain) for funding this research through a ‘Talento’ Fellowship (2016-T1/TIC-1890, PI I. Mendigutía). The research of IM, JGD, and BM is also partially funded by the Spanish “Ministerio de Ciencia, Innovación y Universidades” through the national project “On the Rocks II” (PGC2018-101950-B-100; PI E. Villaver). EK is funded by the STFC (ST/P00041X/1). MV acknowledges the STARRY project, which has received funding from the European Union’s Horizon 2020 research and innovation programme under MSCA ITN-EID grant agreement No 676036. The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement 730890 (OPTICON). This research has made use of the Jean-Marie Mariotti Center LITpro⁴ service co-developed by CRAL, IPAG and LAGRANGE. PMA acknowledges M. Tallon in particular for his support on LITpro, and A.C. Carciofi and R.G. Vieira for their contribution to the observing proposal leading to these results.

Appendix A Log of the observations and interferometric results

Table A.1 includes the target stars, their coordinates and K-magnitudes, the corresponding calibrators and magnitudes, observing dates, average values of the seeing, coherence times, and the telescopes used. Figures A.1, A.2, A.3, A.4, and A.5 show the uv coverage for each target star; Figs. A.6, A.7, A.8, A.9, and A.10 their observed fluxes, squared visibilities, and differential phases from top to bottom for the different baselines. For Fig. A.11, the closure phases per star and baseline are given.

Fig. A.1 uv coverage of the two V590 Mon observations (noted as Obs. A and Obs. B).

Fig. A.2 uv coverage of the two PDS 281 observations (noted as Obs. A and Obs. B).

Fig. A.3 uv coverage of the two HD 94509 observations (noted as Obs. A and Obs. B).

Fig. A.4 uv coverage of the two DGCir observations (noted as Obs. A and Obs. B).

Fig. A.5 uv coverage of the two HD 141926 observations (noted as Obs. A and Obs. B). This is the only star for which the 1.8 m ATs were used instead of the 8.2 m UTs (Sect. 2).

Fig. A.6 Composition of six triplets of panels corresponding to each V590 Mon observation baseline. For each baseline is shown the normalized flux (in red; the same for all baselines), squared visibilities (in blue) and differential phases (in green) in a wavelength range centered in Brγ. Typical uncertainties are indicated in each panel, representing the standard deviation of the corresponding parameter in the adjacentcontinuum of the line.

Fig. A.7 Composition of six triplets of panels corresponding to each PDS 281 observation baseline (see the caption of Fig. A.6).

Fig. A.8 Composition of six triplets of panels corresponding to each HD 94509 observation baseline (see the caption of Fig. A.6).

Fig. A.9 Composition of six triplets of panels corresponding to each DG Cir observation baseline (see the caption of Fig. A.6).

Fig. A.10 Composition of six triplets of panels corresponding to each HD 141926 observation baseline (see the caption of Fig. A.6).

Fig. A.11 Closure phases of the stars in our sample around Brγ line, which center is indicated with the vertical line. Different colors indicate different triplets, as indicated.

References

All Tables

All Figures

	Fig. 1 Squared visibilities measured for the continuum adjacent to Brγ in V590 Mon as a function of the spatial frequency (solid dots with errorbars). The best fit is indicated with a solid line and corresponds to a Gaussian disk on top of the central star. Two additional fits corresponding to ring and uniform disk models are plotted with dashed and dot-dashed lines.
In the text

Fig. 2 Squared visibilities measured for the continuum adjacent to Brγ in DG Ciras a function of the spatial frequency (solid dots with errorbars). The best fit is indicated with a solid line and corresponds to a Gaussian disk on top of the central star, and a background contributing ~ 10%. Two additional fits corresponding to ring and uniform disk models are not plotted as they are practically indistinguishable from the previous.

In the text

Fig. 3 Size-luminosity relation based on data compiled in Tables 3 and 4 for HAeBes. Each panel emphasizes a different aspect, as indicated in the legends. Top left: HAes and HBes; the five with new GRAVITY data presented in this paper are highlighted in red. Top right: the HAeBes with Balmer excess measurements above and below the median value (0.07 magnitudes). Bottom left: the HAeBes with spectro-interferometric measurements indicating that the Brγ/Hα emitting region is smaller than the dust emitting region and the ones without such evidence. Bottom right: Group I and Group II sources based on the SED shape. (All panels): triangles indicate upper and lower limits, and the solid and dashed lines indicate the expected boundaries of rin for an optically thin-MA and optically thick-BL inner gaseous disks, respectively (see text for details).

In the text

	Fig. 4 Inner dust size versus the Hα luminosity, the accretion luminosity, and the mass accretion rate normalized to the stellar luminosity. Linear correlation coefficients (ρ) in each panel quantify the lack of correlation. Upper and lower limits of rin are indicated with triangles.
In the text

Fig. 5 Size-luminosity relation where the stars have been color-coded to indicate different distance ranges as shown in the legend. Interferometrically resolved TTs listed in Pinte et al. (2008) are included in the plot after correction with Gaia EDR3 distances (green open circles), all of them located between 100 and 200 pc from us. The triangles, solid, and dashed lines represent upper and lower limits and possible trends for the dust destruction radius as in Fig. 3.

In the text

Fig. 6 Relation with the distance of the stellar luminosity (top) and the inner dust radius (bottom) for the HAeBes (solid circles and triangles to indicate upper and lower limits), and TTs (open circles) with interferometric estimates of rin discussed in the text. The gray symbols in the top panel represent most HAeBes known to date from the data in GD21. The regions below the dashed lines indicate apparent bolometric magnitudes >12 and angular radii <0.5 mas.

In the text

	Fig. A.1 uv coverage of the two V590 Mon observations (noted as Obs. A and Obs. B).
In the text

	Fig. A.2 uv coverage of the two PDS 281 observations (noted as Obs. A and Obs. B).
In the text

	Fig. A.3 uv coverage of the two HD 94509 observations (noted as Obs. A and Obs. B).
In the text

	Fig. A.4 uv coverage of the two DGCir observations (noted as Obs. A and Obs. B).
In the text

	Fig. A.5 uv coverage of the two HD 141926 observations (noted as Obs. A and Obs. B). This is the only star for which the 1.8 m ATs were used instead of the 8.2 m UTs (Sect. 2).
In the text

Fig. A.6 Composition of six triplets of panels corresponding to each V590 Mon observation baseline. For each baseline is shown the normalized flux (in red; the same for all baselines), squared visibilities (in blue) and differential phases (in green) in a wavelength range centered in Brγ. Typical uncertainties are indicated in each panel, representing the standard deviation of the corresponding parameter in the adjacentcontinuum of the line.

In the text

	Fig. A.7 Composition of six triplets of panels corresponding to each PDS 281 observation baseline (see the caption of Fig. A.6).
In the text

	Fig. A.8 Composition of six triplets of panels corresponding to each HD 94509 observation baseline (see the caption of Fig. A.6).
In the text

	Fig. A.9 Composition of six triplets of panels corresponding to each DG Cir observation baseline (see the caption of Fig. A.6).
In the text

	Fig. A.10 Composition of six triplets of panels corresponding to each HD 141926 observation baseline (see the caption of Fig. A.6).
In the text

	Fig. A.11 Closure phases of the stars in our sample around Brγ line, which center is indicated with the vertical line. Different colors indicate different triplets, as indicated.
In the text