© ESO 2021

1 Introduction

Massive stars (>8 M_⊙) are among the most influential objects in galaxies. Their birth, evolution, and death as supernovae not only affect their immediate vicinity but contribute significantly to the dynamical structure, chemical structure, and evolution of their host galaxies. However, it remains unclear as to what sets the conditions for the formation of a high-mass star, and how the processes of accretion and ejection manifest in the dynamic environments of massive star formation. The main theoretical challenge has been to find a mechanism that can sustain the mass accretion towards the central protostar in the presence of high radiation pressure, which halts the infalling matter (e.g. Wolfire & Cassinelli 1987; Kuiper & Hosokawa 2018). To explain this, the theory is converging towards the formation of massive stars via accretion discs, similarly to their low-mass counterparts (e.g. Kuiper et al. 2010, 2011; Haemmerlé et al. 2017), or via multi-directional mass accretion (Goddi et al. 2020).

The search for circumstellar discs around massive young stellar objects (MYSOs) and determination of the disc properties (i.e. size, mass, infall rates) are essential. Detecting and characterising such discs has been extremely challenging, mainly because of the scarcity and highly embedded nature of their host environments, and the large distance (typically several kiloparsecs) to massive starforming regions. It is mostly thanks to high angular resolution and infrared (IR) interferometry that we are now finally in a position to start revealing and exploring these relatively rare objects in more detail. On large scales (300–2000 au), millimetre line observations (ALMA, SMA, VLA) reveal Keplerian-like disc structures in cold (T ~ 50 K) material (Ilee et al. 2016, 2018a; Johnston et al. 2015, 2020), while substructures down to ~45 au were recently traced by ALMA (e.g. Beuther et al. 2017; Maud et al. 2019). At smaller scales (down to a few au), where the accretion onto the star takes place, the disc is traced using near- and mid-IR (hot/warm) emission (e.g. Boley et al. 2013). Direct evidence of such hot discs is more scarce – there is only one imaged case presented 10 yr ago (Kraus et al. 2010) – because of the extremely challenging nature of the required observations. Therefore, accessing and studying a sample of hot discs around MYSOs at astronomical unit (au) scales is of great importance. In addition, observations of the inner rim of such discs is of particular interest, as this is where interaction between the disc and the star takes place (e.g. magnetospheric accretion, boundary layer accretion, disc-wind; e.g. Mendigutía 2020), and is eventually where planet formation occurs. Beltrán & de Wit (2016) provides a thorough review of discs around luminous YSOs and the processes that shape them. Infrared interferometers trace the size scales and material temperatures of accretion discs around MYSOs uniquely, as well as binarity, highly collimated jets, and stellar winds.

Determining the binary properties of high-mass young stars became particularly relevant when it was found that the binarity fraction of OB-type populations is close to 100% (Chini et al. 2012). Binarity is known to significantly affect the evolution and fate of massive stars (Sana et al. 2012). Multiple systems of massive stars with separations as large as several hundreds of au are mostly predicted by numerical simulations as a result of fragmentation processes during the collapse phase (Myers et al. 2013), while closer binaries (<100 au) may instead form during an accretion disc fragmentation (Meyer et al. 2018) or via orbital decay during internal (e.g. capture in competitive accretion, magnetic braking; Bonnell & Bate 2005; Lund & Bonnell 2018) or external interactions (e.g. with other stars; Bate et al. 2002). However, despite these theoretical findings, reproducible quantitative predictions of massive binary properties are currently lacking. Observational studies of massive binaries and multiple systems in the pre-main sequence (PMS) phase are necessary to inform and distinguish between different scenarios of their formation.

To date, only a few studies have been dedicated to the multiplicity of MYSOs. Pomohaci et al. (2019) performed a search for wide (>1000 au) binary companions of 32 MYSOs, concluding that the total multiplicity fraction of MYSOs may be nearly 100%, and reported mass ratios higher than 0.5, in agreement with studies on intermediate-mass Herbig AeBe PMS stars (Wheelwright et al. 2010). Direct detections of high-mass binaries covering tighter separations (<500 au) were obtained serendipitously. In particular, only a handful of such protobinary systems are known, with separations ranging from some tens to hundreds of au (PDS 27: 30 au, PDS 37: 48 au, V921 Sco: 45 au, NGC 7538 IRS1: 430 au, IRAS 17216-3801: 170 au, IRAS 07299-1651: 180 au; Koumpia et al. 2019; Kraus et al. 2012, 2017; Beuther et al. 2017; Zhang et al. 2019).

Here, we present the first interferometric survey in K-band of six massive YSOs using the unique spatial capabilities of AMBER (2 MYSOs) and GRAVITY (4 MYSOs) on the Very Large Telescope Interferometer (VLTI). The VLTI can achieve angular resolutions down to 1.7 mas with AMBER and GRAVITY on the four 8.2-m Unit Telescopes (UTs), bringing the inner regions of embedded high mass stars within reach. We focus on the characteristic size and geometry of the hot dust traced via the 2.2 μm continuum emission, the ionised gas traced via the hydrogen recombination emission (Brγ), and finally the binarity of MYSOs at milli-arcsecond scales. In Sect. 2, we describe our target selection and the interferometric observations (GRAVITY and AMBER) along with the data-reduction process and the interferometric observables. In Sect. 3, we spatially resolve the hot innermost parts of the MYSOs, and we trace the dust and gas components down to a few au from the central star. In particular, we apply simple geometrical models to fit the interferometric observables and we constrain the size and the geometry of the 2.2 μm continuum emission. In the same section, we constrain the size and the geometry of the Brγ emission towards the MYSOs. In Sect. 4, we investigate our results with respect to the size–luminosity relationship and assess the statistics on the MYSO binarity at mas separations (2–300 au). Lastly, we discuss our results and summarise our findings in Sects. 5 and 6, respectively.

2 Observations and data reduction

2.1 Target selection

Conducting interferometry in the near-IR with the VLTI, and in particular AMBER and GRAVITY, requires relatively bright targets. In the past, most MYSOs were considered faint for VLTI observations on a single-field mode, as the brightest typically have 2MASS K-band magnitudes between 7 and 8, close to the sensitivity limits of the instruments.

To identify program stars, we exploit the largest and most complete sample of MYSOs to date: the Red MSX Source (RMS¹) survey (see Urquhart et al. 2011; Lumsden et al. 2013). This is an unbiased survey of MYSOs throughout the Galaxy and is based on the Midcourse Space Experiment (MSX) which conducted a mid-IR survey of the Galactic plane. The candidate MYSOs have been subject to an extensive multi-wavelength campaign to confirm their MYSO nature, and establish their kinematic distances and bolometric luminosities. The RMS survey also establishes kinematic distances and bolometric luminosities. The selection of targets was bound to near-IR bright MYSOs with 2MASS K-band magnitudes <8, with the stipulation that they must be observable from Paranal, resulting in a sample of 29 MYSOs. Of the objects selected, approximately 50% are not located near a suitable optical guide star, which is a requirement for the wave-front correction by the adaptive optics system MACAO located in the Coudé focus of each UT. In this study, we present the analysis and results of six MYSOs (brightest in the K-band) comprising the largest sample of MYSOs with near-IR interferometric observations to date. These objects are massive (8 M_⊙ < M < 15 M_⊙) and are representative of the unbiased MYSO survey they are drawn from in terms of typical luminosity (L_bol ~ 10⁴ L_⊙) at the average distance of 3.4 kpc. The program stars with successful observations are listed for the whole sample in Table 1. The associated calibrators with their corresponding spectral type, K-magnitude, and size expressed as the diameter of a uniform disc are presented in Table A.1. The calibrators were observed to calibrate the atmospheric transfer function, and were selected based on their single nature, their comparable brightness to the science target, their size (small enough to be unresolved), and their proximity to their associated science target. They were also used during the telluric correction process after removing the observed absorption or emission features around the emission lines of interest (e.g. Brγ).

2.2 GRAVITY observations

Four of our targets, G282.2988, G287.3716, G301.8147, and G034.8211, were observed using the GRAVITY instrument (GRAVITY Collaboration 2017; Eisenhauer et al. 2011) on the four 8.2-m UTs, which operates in the K-band. The observed spectral setup delivered interferometric observables of six baselines in both low (fringe tracker channel) and medium (science channel) spectral resolution simultaneously. The observed projected baseline lengths, B, range between ~40 m and 130 m, corresponding to angular resolutions λ∕2B between ~5.7 mas and 1.7 mas at 2.2 μm, which at the average distance (3.4 kpc) of the present sample of MYSOs corresponds to a physical spatial scale of ~6–20 au. The uv-plane for each MYSO is presented in Fig. A.1.

The interferometric observables were recorded on six baselines simultaneously on the fringe tracker (FT) and the science channel (SC). The SC records the interferometric observables at a medium spectral resolution of R ~ 500 over the entire near-IR K-band window (1.99–2.45 μm), which corresponds to a velocity resolution of 600 km s⁻¹. Typical individual integration times were between 5 and 30 s. The technical overview of the observations including the integration times and atmospheric conditions (i.e. coherence time, seeing) is given in Table A.2. For the reduction and calibration of the observations, the GRAVITY standard pipeline recipes (as provided by ESO, version 1.1.2) were used with their default parameters.

2.3 AMBER observations

Two of the sources in our sample (G233.8306, G231.7986) were observed using the three beam combining instrument AMBER (Petrov et al. 2007) at the VLTI. AMBER was the first generation beam-combiner of the VLTI that used to operate in the near-IR H and K bands as spectro-interferometer, combining beams from three telescopes at three different spectral resolutions (35, 1500, 12 000), and it was decommissioned in 2017. The selected configuration for G231 and G233 used the 8.2 m aperture of three UTs (UT 2, 3, and 4), delivering projected baselines of between ~45 m and ~88 m. The uv-plane is presented in Fig. A.2. The fringes were recorded spatially, that is, each detector integration contains fringe pattern, also called an interferogram. The fringe patterns for each of the three baselines are detected co-spatially (incoherent addition of the light) and are therefore encoded non-redundantly (i.e. with distinct fringe pattern frequencies) such that they can be recuperated by means of a Fourier transform of an interferogram (see Robbe-Dubois et al. 2007).

For this program, fringes could be obtained with AMBER setup in its low spectral resolution mode which has R = 35. No use could be made of the FINITO fringe-tracker as the deeply embedded sources are too faint for the H-band atmospheric window. Each observation of a program star was bracketed by a calibrator star. The meteorological conditions of the observing night are reflected in Table A.3 by the coherence time (τ₀) and seeing measurements from the Differential Image Motion Monitor (DIMM) on Paranal. The values in Table A.3 indicate worse than average Paranal conditions and it should be clear that they were not favourable to obtain high-quality spectro-interferometric measurements. For the further analysis of this dataset we focus on the visibilities extracted for the continuum. The raw AMBER dataset was reduced using the amdlib data reduction package (version 3.0.9 Tatulli et al. 2007; Chelli et al. 2009). During the data reduction we selected the best 20% of the observed frames, as has been previously suggested to provide robust visibilities (Malbet et al. 2007).

Fig. 1 Normalised spectra of the MYSOs observed with GRAVITY around the 2.2 μm continuum level. The wavelength coverage contains the Brγ which is seen in emission towards all sources, and the CO bandheads which are not detected towards G301. G301 also shows some weak emission within the 2.05–2.1μm range (Fe II, He I). The narrow absorption features are spectral artefacts caused by telluric correction.

2.4 Observational results

For all six MYSOs of our sample (GRAVITY and AMBER) we were able to measure the visibilities of the 2.2 μm continuum. In addition, for the four sources observed with GRAVITY, we were able to further extract information on the spectra, the visibilities around the Brγ emission, and the closure phases of multiple telescope triangles. We note that the GRAVITY dataset (4 sources) is superior to AMBER (2 sources) with respect to spectral resolution (500 versus 35) and uv-coverage (4 UTs versus 3 UTs), which allowed the extraction of the additional information.

2.4.1 K-band spectra

We present new near-IR spectra of G282.2988, G287.3716, G301.8147, and G034.8211 as observed with GRAVITY (Fig. 1). The observed wavelength coverage includes the Brγ hydrogen recombination line emission at 2.167 μm, which is detected in emission in all sources. The GRAVITY spectra also cover the CO bandheads, which are found in emission towards all sources but G301 (Fig. 1). We note that the CO bandheads were previously reported by Ilee et al. (2014) towards G282.2988 and G287.3716. In our study, we do not focus on this molecular emission.

For each source, we fit a Gaussian distribution to the line profile of the Brγ emission and derive its line-to-continuum ratio at line peak (L/C), and the full width at half maximum (FWHM). G301.8147 shows the strongest Brγ emission among the sources, with L/C at line peak of 1.29. G282.2988 is characterised by a Brγ L/C at line peak of 1.17, followed by G287.3716 with a L/C at line peak of 1.12. G034.8211 is the source with the weakest Brγ emission, and in particular it is characterised by a L/C at line peak of 1.06. The wavelength calibration uncertainties and the low spectral resolution do now allow reliable measurements of peak velocities.

We note that for all sources, the FWHM ranged between ~550 and 596 km s⁻¹, which is at the spectral resolution limits of the instrumental observing mode (~600 km s⁻¹), and therefore Brγ is spectrally unresolved. Interestingly, the only source of the sample without signs of the CO bandheads (G301.8147) shows the strongest Brγ emission. Increasing accretion rates (e.g. stronger Brγ; see, Mendigutía et al. 2011) or evolutionary stage (Cooper 2013; Cooper et al. 2013) are known to influence the detection rate of the CO bandheads towards MYSOs (Ilee et al. 2018b).

2.4.2 Visibilities and phases

We extracted the visibilities, differential, and closure phases (CP) from the GRAVITY and AMBER observations around the 2.2 μm continuum and Brγ emission. We note that AMBER combines three telescopes instead of four resulting in a single closure phase. In general, the visibility values are lower for larger angular geometries. In practice, when comparing the geometries of sources directly by means of visibilities, the baseline position angle comes into play. Relative brightness distributions of distinct emitting regions also affect the contrast levels per baseline length. Phases can be used to assess the degree of symmetry of the emitting region. Variations in differential phases between the continuum and a line emission indicate differences in the photocentre of the two emissions. On the other hand, a non-axisymmetric brightness distribution of emission will result in non-zero closure phases. The differential phase of the Brγ with respect to the continuum is typically 0° within the error bars (see Figs. A.3–A.6), with possibly the only exception for the longest baselines towards G301. In this work, the differential phases are not taken into further consideration. We focus on the visibilities and closure phases of the continuum for all six MYSOs and of the Brγ emission for four MYSOs (GRAVITY).

2.4.3 2.2 μm continuum

All the sources in the sample (Table 1) are fully or partially resolved, showing continuum visibilities V_cont < 1 even at the shortest baselines (~40 m). The measured visibilities generally decrease with increasing baseline length (see also Table A.2). Before proceeding with a detailed geometric modelling of the emission, we extract the measured visibilities of the continuum as averaged over a bandwidth of ~0.03 μm bluewards and ~0.03 μm redwards of the Brγ emission². The measured visibilities, in addition to the angular size of an emitting region, are also linked to the position angle of the baselines and geometry (inclination) of the disc (e.g. G034.8211).

In addition, we can assess the degree of symmetry or asymmetry of the emitting region by extracting the closure phases. For near axisymmetric brightness distributions, the closure phases are always ~0° (or 180°), while any other measured value is indicative of a skewed intensity distribution. G034.8211 is the only MYSO in the sample that is characterised by asymmetric continuum emission at both the smallest (~1.7 mas) and intermediate scales (~3 mas) (Table 2).

2.4.4 Brγ emission

The interferometric observables of GRAVITY as a function of wavelength of the Brγ emission and the continuum around it for the four MYSOs are presented in Figs. A.3–A.6. A summary of the findings per source is given below.

G282.2988: The Brγ emission shows an increase in visibility of 1.3–6% (highest difference at longest baseline; ~130 m) with respect to the continuum, while the associated closure phases at longer baselines are ~0.5–1° indicating a mostly symmetric emission of the ionised gas.

G301.8147: The Brγ emission is symmetric (CP~0°) and shows an increase in visibility between 22% and 47% (highest difference at short baselines; ~50 m).

G034.8211: The visibility of the Brγ emission follows the visibility of the continuum at the long baselines (90 m, 130 m) and is higher than the continuum by ~2–12% for the rest of the baselines. In addition, the Brγ emission shows closure phases up to ~6–10° following those of the overall continuum emission.

G287.3716: The Brγ visibility increases by 7–10% at two short baselines (~50 m) and shows no significant changes for the rest of baselines. The observed closure phases are <4° following the symmetric nature of the continuum (we note the large errors of 2–5°).

We conclude that for all four MYSOs in the GRAVITY sample, Brγ shows a similar or a higher visibility value compared to the continuum for different baselines, which is beyond the associated errors and therefore can be attributed to real geometrical effects. Hence, the observed changes in visibilities indicate that the ionised region traced via the Brγ emission is of comparable size or is more compact compared to the continuum-emitting region. For G034.8211 and G301.8147, the Brγ and continuum also show some variations in differential (i.e. differences in the photocentre of the continuum and that of the line-emitting region) and closure phases (i.e. asymmetries, Fig. A.3). The spatial distribution of the Brγ emission appears to follow the morphology of the continuum for the symmetric sources (observed variations are within uncertainties and systematic errors), and it shows a more symmetric signature compared to the continuum for G034.8211.

3 Modelling the brightness distribution

To determine the size of the hot dust (2.2 μm continuum) and ionised gas (Brγ) emission, we adopt simple geometries to fit the observed visibilities of the sample of six MYSOs. We model the visibilities using the fitting software LITpro³ (Tallon-Bosc et al. 2008).

3.1 Size estimations: Inner disc

Observing and measuring the innermost radius of discs (<10 au) towards MYSOs is very important as this is the region where the interaction between the disc and the central star is most prominent, with material directly feeding the central star (e.g. via magnetospheric or boundary-layer accretion). The present dataset gives us access to scales of only a few au at the typical distances of the current sample of MYSOs (~3.4 kpc). In addition, the K-band continuum flux can be mostly attributed to the thermal hot (~1500 K) dust emission, which is also the temperature at which silicate (Si) grains typically sublimate (Kessler-Silacci et al. 2007), while these graphite (C) grains sublimate at higher temperatures (2000 K; Baskin & Laor 2018). The size of the 2.2 μm continuum emission has been traditionally used to directly probe the dust sublimation radius due to the radiation from the host star towardsT Tauris and Herbigs (e.g. Monnier & Millan-Gabet 2002). Here, we present direct measurements of the K-band continuum size towards a sample of MYSOs, and we place MYSOs with K-band measurements in the context of a size–luminosity diagram for the first time (Sect. 4; for M-band see Grellmann et al. 2011).

To determine the characteristic size of the 2.2 μm continuum emission, we translate the observed visibilities of each MYSO to angular sizes. To do so, we apply three simple geometric models of a pre-defined brightness distribution, a Gaussian, a uniform disc, and a ring, and fit the observed to the predicted visibility curves.

We ran a grid of sizes between 0.2 mas and 10 mas. Figure 2 presents the observed visibilities obtained with GRAVITY and AMBER respectively, overplotted with the best-fit modelled visibilities (lowest reduced χ²) for all MYSOs in our sample. The estimated sizes per source for the different brightness distributions are presented in Table 3. For four out of the six sources in our sample, the addition of a background component (or halo) was required to explain the observed low visibility values (<0.8) at the shorter baselines. Lastly, for the sources in our sample for which pole-on geometries resulted in a reduced χ² > 10, an additional centred point source was introduced in the fitting process. The flux weights of the additional components (1 for point source; 3 forbackground) are also presented in Table 3.

We report the 2.2 μm emission sizes to be between 1.7 mas and 6.2 mas with typical errors of between 1 and 6%. Depending on the adoptedbrightness distribution, the size estimations of a single source can vary by up to a factor of ~1.8 (e.g. G034.8211), while the reduced χ²s are mostly below 10. Overall, the lower the reduced χ², the smaller the size variations among the models. For one MYSO, G034.8211, the shape of the visibility curve does not show a gradual decrease in visibility with increasing baseline (wiggly shape, Fig. 2 a), and can only be fitted by assuming a flattened brightness distribution (reduced χ² < 30 versus >80). The best fit could be achieved for a flatten ratio of ~1.5 (ratio between major and minor axes) regardless of the initial assumption of a disc, a Gaussian, or a ring distribution. The position angle is between 20 and 34 degrees depending on the adopted brightness distribution. G282.2988 was the only MYSO for which a binary brightness distribution was necessary to improve the goodness of the fit (reduced χ² ~ 2 versus 12; see Sect. 4.2.2). We additionally confirm the estimated sizes of the 2.2 μm continuum emission of the four MYSOs observed with GRAVITY by performing a model-independent image reconstruction where the closure phases are also taken into account (Appendix B). The image reconstruction of G034.8211 revealed a similar elongated brightness distribution to the one retrieved from geometric modelling.

Fig. 2 Modelled visibilities of the 2.2 μm continuum emission as function of baseline assuming a single disc brightness distribution overplotted with the observed ones as obtained with GRAVITY for (a) G034.8211, (b) G301.8147, and (c) G287.3816 as fitted using LITpro (Sect. 3.1). (d, e) Same as before but for observations obtained with AMBER towards G233 and G231, after fitting a Gaussian, a uniform disc, and a ring. The background flux contribution is found to be 16% for G233 and ~9% for G231.

3.2 Size of the Brγ line-emitting region

The origin of the Brγ emission observed in the close vicinity of YSOs and MYSOs (Mendigutía et al. 2011; Murakawa et al. 2013; Davies et al. 2010; Kraus et al. 2008; Pomohaci et al. 2017; Lumsden et al. 2012) can be explained by theory (e.g. jets, magnetospheric accretion, disc wind; Ferreira 1997; Tambovtseva et al. 2014, 2016) and can be further tested by observationally constraining the spatial origin of the emission compared to the size of the hot dust (2.2 μm emission). The GRAVITY dataset allows us to estimate the size of the ionised gas by measuring visibilities at the central spectral channel of the Brγ emission towards four MYSOs and to directly compare it to that of the continuum emission. We note that the low spectral resolution of this dataset (~600 km s⁻¹) prevents us from studying a possible stellar wind, disc-wind, or a jet origin for the Brγ emission kinematically (typical FHWM of 100–200 km s⁻¹; Bunn et al. 1995).

To estimate the size of the ionised emission, we used the measured calibrated visibilities at the spectral channel that corresponds to the peak of the Brγ line profile. Before modelling the visibilities of any line emission, we need to correct for the continuum visibility and flux contributions (see also, Malbet et al. 2007). In particular, we applied Eq. (1), $$V_{\text{line+cont}}=\frac{V_{\text{cont}}\times F_{\text{cont}}+V_{\text{line}}\times F_{\text{line}}}{F_{\text{cont}}+F_{\text{line}}},$$(1)

where the total visibility (V_line+cont) is a function of continuum and line visibilities (V) and fluxes (F), and solved for V_line.

We derive the size of the Brγ emission by fitting the interferometric observables with a simple geometrical model. In particular, we fit three different brightness distributions, a Gaussian, a uniform disc, and a ring, for a range of sizes (0.2–10 mas), similarly to how we treated the 2.2 μm continuum emission. Figure 3 shows an example (G034) of the observed visibilities of the continuum and the Brγ overplotted with the corresponding best fit of the continuum. The estimated sizes per source for a given brightness distribution are presented in Table 3. The best fit of each adopted distribution demonstrates that the ionised gas systematically originates from a similar or smaller region (up to ~10%) compared to the continuum for three out of four MYSOs observed with GRAVITY.

Arguably, underlying photospheric Brγ absorption is known to affect the measured continuum level and Brγ, and therefore the size estimates of the Brγ emission could be overestimated (Eq. (1)). Such contributions may become more significant for G301 and G282, for which the resolved emission accounts for up to ~50% and 10%, respectively (Table 3), and less for G034 and G287, for which the stellar contribution is unresolved⁴. We note that, although such contributions may affect the absolute measurements of the Brγ sizes (we therefore provide upper limits), the qualitative findings of the smaller Brγ sizes with respect to the continuum emission will remain. A detailed study of those effects is beyond the scope of this study.

Fig. 3 Modelled visibilities and closure phases of the 2.2 μm continuum after assuming a flattened disc, overplotted with those observed towards G034.8211. The observed visibilities of the Brγ emission arealso overplotted in red, demonstrating the co-planar origin of the continuum and the ionised gas.

4 MYSOs in a broader context

4.1 Size–luminosity diagram

To proceed with the size–luminosity relation, we adopt the sizes of the ring brightness distribution determined in Sect. 3.1. Even though a ring does not always reflect the best fit of the interferometric observables, it is preferred for the purpose of this section for consistency with other studies of this nature. In particular, a ring is generally assumed to reflect the size of the inner disc where dust sublimates (see also, Monnier & Millan-Gabet 2002). In addition, the adopted brightness distribution can be supported by theory, because the hot dust emission is expected to stem from a narrow disc annulus at distances very close to the star (see also Fig. 4 and Kraus et al. 2010; Stecklum et al. 2021). The morphology and properties of the inner rim around low- and intermediate-mass YSOs has been the topic of multiple studies (Isella & Natta 2005; Tannirkulam et al. 2007; Kama et al. 2009; McClure et al. 2013). Here, we directly compare the measured ring sizes of the five MYSOs in our sample (excluding G282 which can be best fitted as a binary) with their associated luminosities (Table 1; taken from Mottram et al. 2011; Guzmán-Díaz et al. 2021; Wichittanakom et al. 2020). We investigate the location of MYSOs in the size–luminosity diagram with respect to the dust sublimation radius as predicted from theory. In particular, we investigate the dust sublimation radius as predicted by (a) oblique heating of an optically thick flat disc (classical disc; Hillenbrand et al. 1992; Millan-Gabet et al. 2001; Monnier et al. 2005), (b) direct disc heating from the star in the presence of an optically thin cavity ignoring backwarming effects (Tuthill et al. 2001), and (c) direct disc heating from the star in the presence of an optically thin cavity, but this time taking into account backwarming effects from the hot dust (self-irradiation; Dullemond et al. 2001). Lastly, we compare our results with findings from studies of the MYSOs’ less massive counterparts (Herbig AeBes, T Tauris).

Fig. 4 Modelimage of the 2.2 μm continuum emission towards M8EIR (based on Frost et al. 2021) convolved at 1.7 mas resolution (GRAVITY/VLTI). The distance to M8EIR is 1.3 kpc.

4.1.1 Massive young stellar objects

By plotting the measured 2.2 μm sizes of the present sample of MYSOs (excluding G282 because of its binarity) as a function of their stellar luminosity (Fig. 5a), it becomes apparent that there is a large scatter of inner disc radii at the luminosity range (1.1 × 10⁴ L_⊙–2.4 × 10⁴ L_⊙) of this sample. To better explain the observed sizes with respect to stellar luminosities, we overplot thepredicted dust sublimation radius of a disc with an optically thin inner cavity ($R_{\text{s}}\propto L_{*}^{1/2}$; Tuthill et al. 2001)for a range of dust sublimation temperatures (1000 K–1500 K Kessler-Silacci et al. 2007; Boley et al. 2013). With the exception of G231, the rest of the MYSOs appear to follow what is predicted by the models. G231 is the only source in our sample, which appears significantly smaller and its size cannot be explained by the optically thin scenario (Fig. 5a).

To investigate the size–luminosity relation towards a larger sample of MYSOs, in addition to our direct measurements of K-band inner radii, we overplot the other two MYSOs from the literature with K-band interferometric observations and available 2.2 μm continuum size measurements (IRAS 13481-6124, NGC 2024 IRS2, Kraus et al. 2010; GRAVITY Collaboration 2020a). Moreover, we overplot indirect measurements of the inner disc radii of eight MYSOs (Fig. 5b) derived by Frost et al. (2021). These latter authors performed advanced radiative transfer modelling to simultaneously fit high-angular-resolution interferometric observations at mid-infrared (MIDI/VLTI), images (VISIR/VLT, COMICS/Subaru), and spectral energy distributions towards their MYSO sample. A good fit of the MIDI visibilities of five of those additional MYSOs at the shorter wavelengths could only be achieved when their inner radius was set to be significantly larger than the sublimation radius predicted from the optically thin models. To verify this finding, we extracted the modelled images at the wavelength of interest (K-band; 2.2 μm) and convolved them to match the spatial resolution of GRAVITY/VLTI (1.7 mas). The measured size of the modelled 2.2 μm continuum emission for all eight sources is in agreement (within 0.5–2 au, therefore with the errors reported in Frost et al. 2021) with the inner radius predicted by models based on longer wavelength observations (e.g. M8EIR; Fig. 4). When these MYSOs are placed in a size–luminosity diagram, as we show in Fig. 5, self-irradiation (blue shade) cannot explain the sizes of their inner radii, but direct K-band interferometric measurements are necessary to independently confirm this finding.

The combined sample of MYSOs starts revealing a trend of increasing inner disc radius with increasing luminosity, and in particular the sizes scale with the square root of the stellar luminosity. We further discuss the location of MYSOs with respect to models in Sect. 5.1.

4.1.2 Herbig AeBes and T Tauris

To directly compare the size–luminosity relation of the MYSOs to YSOs of lower mass, we collect and overplot the measured sizes of the 2.2 μm emission of known T Tauris and Herbig AeBes in Fig. 5c. The initial luminosities and inner sizes are obtained from Millan-Gabet et al. (2007), Monnier & Millan-Gabet (2002), and Pinte et al. (2008, and references therein), the most recent GRAVITY YSO survey (GRAVITY Collaboration 2019), and the five new Herbig measurements presented in Marcos-Arenal et al. (2021). The luminosities and sizes of all objects are scaled taking into account the new distances obtained using Gaia EDR3 parallaxes (Guzmán-Díaz et al. 2021), while no binarity at the traced scales was reported for these sources.

In Fig. 5c, we present the size–luminosity diagram with respect to the dust sublimation radius as predicted by three different disc models for dust sublimation temperatures of 1000 K and 1500 K. Figure 5 shows that with the exception ofa few Herbig AeBes which appear smaller than the rest (classical disc regime), the observed inner radii of most YSOs followthe distribution of the dust sublimation radii of an optically thin discs, increasing with the square root of the luminosity. This trend can be seen within a wide range of luminosities. A more detailed discussion on the sizes of different classes of objects with respect to different models is presented in Sect. 5.1.

4.2 Binarity

To inform and differentiate among the theories of high-mass binary formation, studies need to provide observational information on the frequency, separation, and mass ratios of binary MYSOs. The present work targets separations on scales of a few au to a few hundred au, and provides a bridge between adaptive-optics-assisted imaging targeting 600 au to 10 000 au separations (Pomohaci et al. 2019) and a high-resolution spectroscopic survey of MYSOs which indirectly probes the closest separations at sub-au scales (Shenton et al., in prep.; for massive young stars see also, Apai et al. 2007).

4.2.1 Our sample with respect to binarity

In this study, we undertake the first companion search at milli-arcsecond scales towards a sample of MYSOs. We investigate whether our sample of six MYSOs at 2.2 μm with VLTI (GRAVITY or AMBER) contains close binaries at mas separations covered by the field of view and angular resolution of the present interferometric K-band observations. The VLTI-UT configuration at 2.2 μm is sensitive to binary separations of 0.5 milli-arcsecond up to the single telescope diffraction limit (~70 mas). For a typical distance of 3.4 kpc this corresponds to 1.7 au and 240 au, respectively. Any wider pairs should already be traced by the NaCo imaging survey. We note that GRAVITY is sensitive to companion detections of a maximum magnitude difference ΔK ~ 5 mag within a 3–50 mas range of separations (see also, Sanchez-Bermudez et al. 2017), similar to the one achieved by Pomohaci et al. (2019) in the 1–3″ separation range. For the given sensitivity limits, we would not be able to detect subsolar mass companions.

With the exception of G034.8211, all objects were also part of a survey of 32 MYSOs using adaptive-optics-assisted high-resolution K-band imaging with NaCo on the 8.2 m VLT (Pomohaci et al. 2019). Three objects in our sample (G301.8147, G287.3716, and G282.2988) were found to be part of a multiple with wide companion separations of between 1.8″ and 2.9″, while the other two objects showed no multiplicity. We note that 30% of the full sample of 32 objects was observed to be binary with physical separations of between 400 and 46 000 au (0.6″–3.1″).

Fig. 5 (a) Near-IR size–luminosity diagram of our observed sample of MYSOs. The measured inner ring size Rin of the observed MYSOs (GRAVITY, AMBER; based on Table 3) is plotted with respect to the stellar luminosity. The shaded area represents the dust sublimation radius Rs as predicted by the presence of a disc with an optically thin cavity for a range of temperatures. The adopted dust sublimation temperatures are 1000 K (solid line) and 1500 K (dotted line). (b) Same as before but overplotted with two MYSOs with K-band measurements from the literature (filled black stars; IRAS 13481-6124, NGC 2024 IRS2, Kraus et al. 2010; GRAVITY Collaboration 2020a) and with indirect measurements (star symbols; modelled sizes presented in Frost et al. 2021). The open and filled stars represent the sources for which the modelled inner radius was found to be larger than or equal to Rs respectively. The dust sublimation radius as predicted by models of different dust temperatures with (blue area) and without (yellow area) backwarming effects is also overplotted. (c) Same as before but including the available measurements of Herbig AeBes and T Tauris from the literature. The dust sublimation radius as predicted by the classical scenario of a flat optically thick disc is also overplotted.

Fig. 6 Modelled visibilities and closure phases of G282.2988 as fitted using a binary model in LITpro. The best fit reveals a component separation of 23.7 mas (consistent with Koumpia et al. 2019) and primary to secondary flux ratio of ~10.

4.2.2 Method used to find binaries

To investigate binarity among the MYSOs in our sample, we proceeded by evaluating the addition of an off-centre source component in the geometrical models presented in Sect. 3 with respect to the resulting χ². Our models revealed that single brightness distributions are sufficient (χ² < 3 for four objects with pole-on geometries) to reproduce the interferometric observables for most MYSOs. Modelling G233 and G231 as single sources resulted in areduced χ² of ~2, which, in combination with the limited number of visibility measurements, makes the inclusion of an off-centre point source unjustified. For G301 and G287, which come with a higher number of visibility or closure phase measurements, a reduced χ² as low as 1.8–3 could be achieved without adding an off-centre source. Introducing an off-centre component for G034 did not result in an improvement of the fit (χ² > 10).

G282.2988 (known also as PDS 37) is the only object for which both an elongated disc and the addition of an off-centred source improved the fit significantly (Fig. 6), reducing the χ² from 12 to 2–3. The source was recently identified as binary at shorter H-band wavelengths using PIONIER/VLTI (PDS37; Koumpia et al. 2019). The K-band modelling reveals a companion at 23.7 ± 0.1 mas with a primary-to-secondary flux ratio of 10.4 ± 1.2 and a position angle PA of 260 ± 2°. The fainter secondary companion is found to be approximately five times more extended than the primary, resulting in a χ² of 3 compared to ~70 when a point source is considered. Although the brightness is a good proxy for mass at long wavelengths, the emission at 2.2 μm in such embedded environments is still affected by significant extinction. The presence of a large disc surrounding the fainter object may suggestthat the companion is more massive but embedded. Both the reported separation and the PA are consistent with Koumpia et al. (2019), making the binary nature of G282 a favourable geometry over that of an elongated disc.

4.2.3 Binary fraction

In our sample of six MYSOs (AMBER and GRAVITY), we find a MYSO binary fraction of 17 ± 15% in K-band. Two of the sources in our sample are observed only with AMBER, and therefore are limited to only three measured visibilities for baselines of between ~40 m and 80 m. We note that the interferometric binary signal for G282.2988 becomes apparent for baselines >80 m (<2.8 mas; ~10 au), and would have been flagged as a non-binary source in our AMBER observations as the angular resolution provided by AMBER on UTs would not have been sufficient to resolve it. Therefore, this limitation should be considered for the two sources of our sample observed with AMBER. G231.7986 (PDS 27) in particular, which is found to be a single brightness distribution based on AMBER observations, was found to show strong binary interferometric signatures in PIONIER H-band observations (Koumpia et al. 2019). PDS 27 is the only source in our sample presented in GRAVITY Collaboration (2019). In this latter study, PDS 27 shows closure phase variations of up to a few degrees, but the authors argue that there is no need to model a binary companion. We note that the size of PDS 27 as obtained with AMBER (1.8 mas) is consistent with the size obtained with GRAVITY (1.66 mas GRAVITY Collaboration 2019). The fact that a binary was not detected in the case of PDS 27 with AMBER or GRAVITY could be due to the different fluxes (i.e. different masses, evolutionary stages) of the candidate companions in K- or H-band. If we take into account a more homogeneous observational dataset by focusing on the GRAVITY sample alone, the MYSO binary fraction is 25 ± 21%. A more thorough discussion on the statistics and a comparison with the literature are presented in Sect. 5.3.

5 Discussion

We discuss the relation between the luminosity of the central star and its associated near-IR size towards MYSOs. In addition, we discuss our findings on MYSO binarity at milli-arcsecond scales (2–300 au) for this sample of 6 MYSOs.

5.1 The size–luminosity relation

We observe a large scatter of about an order of magnitude in the near-IR size at similar luminosities for the classified Herbig Be stars and MYSOs (10⁴ L_⊙ < L_* < 10⁵ L_⊙). Numerous factors have been explored to explain the different sizes of the K-band continuum emission (i.e. radius of the inner rim) at a certain luminosity, both in terms of predictions and observations. The predicted dust sublimation sizes were found to be affected by the size, opacity, and composition of grains (e.g. Monnier & Millan-Gabet 2002), the nature of the accretion (e.g. turbulence Kuchner & Lecar 2002), and photoevaporation (Danchi et al. 2001). Observations on the other hand can introduce uncertainties on the interpretation of the luminosity (e.g. multiplicity; Hartmann et al. 1993) and near-IR size measurements (i.e. geometry of the brightness distribution). More recently, Marcos-Arenal et al. (2021) performed an investigation of the size–luminosity correlation and the observed size scatter towards Herbig AeBes. Here, we investigate the influence of backwarming effects and accretion on the location of MYSOs on the size–luminosity diagram.

5.1.1 Backwarming effects

It is known that for geometries where the vertical height of the inner rim is not significantly smaller than the radius of the inner rim, the backwarming by the circumstellar hot dust is not negligible, resulting in an increase in the dust sublimation radius (self-irradiation; Dullemond et al. 2001). Here, we explore this mechanism with respect to the location of the different classes of objects in the size–luminosity diagram. The addition of backwarming results in a larger sublimation radius compared to the dust sublimation radius predicted by the classical flat disc with an optically thick inner region of gas or the disc with an optically thin cavity of gas without backwarming effects (Fig. 5c).

The sizes of the T Tauri stars in Fig. 5c can be better explained with the dust sublimation radius predicted by a disc with an optically thin cavity, where the inner rim is directly heated by the stellar radiation, and taking into account backwarming. The Herbig Aes and the more luminous Herbig Bes follow the size distribution as predicted by an optically thin disc, but this time neglecting the backwarming effects. This finding indicates that the small grains are a more prominent backwarming heating mechanism for low-mass stars but that they become less important with increasing mass. The situation is less clear for the regime of the more luminous Herbig Be stars and MYSOs. Although most of the sizes of the inner rim can be explained by the optically thin disc with or without the backwarming effect, five of the MYSOs with indirect measurements are larger and cannot be explained by any of the models, while one MYSO and four of the Herbig Be stars appear undersized. The Herbig exceptions instead follow the size distribution expected by the classical scenario of a flat optically thick disc with oblique heating (Hillenbrand et al. 1992).

The two most luminous MYSOs in our sample (G034 and G301) have a larger K-band size (i.e. inner radius) than the rest, and are located in the lower temperature regime (1000 K) of an optically thin disc when backwarming effects are neglected, or the higher temperature regime (1500 K) of an optically thin disc if backwarming effects are included (Fig. 5b). Given that a temperature of 1000 K is very low for dust sublimation to occur, we consider the scenario of backwarming effects becoming more important for these two MYSOs more likely. In contrast, the inner disc sizes of G287 and G233 follow the distribution of the dust sublimation predicted by an optically thin disc assuming a dust sublimation temperature of 1500 K and neglecting backwarming effects. G231 appears to have an undersized inner radius compared to the rest of the MYSOs which cannot be explained by either of the two optically thin scenarios, while it appears oversized compared to the predictions of an optically thick disc.

We conclude that introducing backwarming effects in the discs surrounding YSOs appears to have a prominent role in explaining their location on the size–luminosity diagram, but there is no clear trend or explanation as to why this mechanism seems to be more prominent for some objects than others. Differences in dust grain compositions and more complex inner disc geometries (e.g. a flat, exposed inner rim) may contribute to the differences in backwarming contributions. However, a detailed exploration of those effects would require a dedicated modelling study.

5.1.2 Evolution and accretion

To investigate the possible influence of the evolutionary status on the different locations of the MYSOs on the size–luminosity diagram, we retrieved the available evolutionary class of G034 and G231 from Cooper (2013) and Cooper et al. (2013). The classification scheme presented in Cooper (2013) is based on the morphology of the near-IR spectra and defines three types (with their own subtypes) of (M)YSOs: Type I, II, and III. Sources of Type I are the youngest (redder among subtypes) and show strong H₂ emission and no ionised lines. Type III sources are the oldest (most blue subtype), and they show strong H I lines, prominent fluorescent Fe II emission at 1.6878 μm emission, and, if any, weak H₂ emission. Based on that scheme, G034 was classified as a class II object, while G231 was classified as a more evolved, IIIb type. Further near-IR spectra (see Fig. C.1) of G287, G231, G233, G282, and G301 were obtained using IRIS-2 on the AATin a wider-reaching project related to RMS between 2006 and 2008. The spectral resolution in the H and K bands was ~2400, allowing clearer identification of weak lines than the data presented in Cooper (2013) and Cooper et al. (2013). Four of the objects, G231, G233, G282, and G301, show strong Br series, indicating a class III designation in Cooper’s evolutionary scheme. Of these, only G282 shows evidence of H₂ emission, which makes it class IIIa, while the other three MYSOs are class IIIb. The class IIIs also show evidence fora stronger ionising continuum through the presence of fluorescent FeII; in particular the IIIbs show both the 1.688 and 2.089 micron lines. G233, G282, and G301 also show [FeII] emission, showing that shocked gas in an outflow type process is still present. G287, like G034, can be classified as Type II given the stronger H₂ emission and much weaker HI emission.

Following the same reasoning as for Herbigs, we would expect the less evolved objects with higher accretion rates, namely G034 and G287, to follow a regime of smaller inner disc radii in the size–luminosity diagram, and the remaining, more evolved MYSOs to follow a regime of larger sizes. This is in contrast to what we observe for four out of the six sources, with G287 and G301 being the only two MYSOs following that expectation. Based on these observations, it appears that similar to Herbigs, the different evolutionary stages alone cannot explain the observed discrepancy in the measured 2.2 μm sizes.

In addition, three MYSOs with indirect measurements (Frost et al. 2021) follow the same relation as our more luminous MYSOs (optically thin disc with backwarming effects), while five of those objects are systematically larger and cannot be explained by the models. The larger sizes of the inner radius compared to the dust sublimation radius for those MYSOs are attributed to increasing inner holes with age due to photoevaporation or the presence of binary companions (Frost et al. 2019, 2021). Figure 5c shows that the distribution of the locations of the oversized Herbigs and MYSOs appears to follow a trend, although this cannot be explained from the current models. We note that even if one assumes a dust sublimation temperature of 2000 K as an upper limit, the predicted sizes of the inner radius would get smaller increasing the discrepancy even more. To verify that the observed trend is realand investigate its origin in more detail, direct measurements of the inner sizes (traced in the K-band) of the sample presented in Frost et al. (2021) are necessary.

To explain the observed scatter in sizes at given luminosities (Fig. 5), we also investigate the impact of accretion rates. High accretion rates may result in optically thick gaseous environments, shielding the dust in the inner radius. In that context, objects with higher accretion rates would be characterised by a smaller inner radius compared to the rest at a given luminosity, therefore deviating from the general optically thin trend. This argument could be used to explain the one undersized MYSO and the undersized Herbigs (Muzerolle et al. 2004). Both in theory and observations, the low-mass star formation is characterised by accretion rates of 10⁻⁹–10⁻⁷ M_⊙ yr⁻¹ (T Tauris regime, Ingleby et al. 2013; Hartmann et al. 1998), while high-mass star formation generally requires accretion rates that are three to four orders of magnitude higher (Hosokawa & Omukai 2009). We note that star formation is characterised by variable rather than steady accretion (Vorobyov 2009). The observed discrepancy in the size–luminosity diagram for massive objects of similar luminosities could be indicative of less evolved massive young stars that are naturally characterised by higher accretion rates. If we extend this argument to Herbigs, less evolved objects could then explain their location in the regime where sublimation occurs for the optically thick disc scenario. Marcos-Arenal et al. (2021) investigated the size–luminosity distribution of HAeBes with respect to the nUV Balmer excesses, the Hα and accretion luminosities, and the mass accretion rates, and they did not report a clear trend. In addition, the presence of emission from the CO bandheads, which were previously found to originate from an inner gaseous disc (Ilee et al. 2014; GRAVITY Collaboration 2020a), could act as a shielding mechanism for the dust, allowing it to survive at distances closer to the central star and therefore resulting in a smaller inner radius of the dusty disc (i.e. smaller 2.2 μm size). Based on the near-IR spectra of this sample, the presence of CO could explain the small sizes of G287 and G233, but not the large size of G304, while similarly the absence of this molecular emission could explain the large size of G301 but not the small size of G231.

In conclusion, we find that the dust inner rim radius of MYSOs, when directly measured via the 2.2 μm emission, does not show a clear trend with respect to the stellar luminosity. When MYSOs are treated as a class and compared to low luminosity T Tauri, Herbig Ae, and most of the Herbig Be stars, then we observe a general trend of increasing inner rim radius with the square root of the stellar luminosity. This finding is suggestive of a universal trend in the observed size–luminosity diagram, indicating that the sizes of inner regions of discs around young stars scale with luminosity independently of the stellar mass, and are consistent with the dust sublimation radius predicted by models.

5.2 Origin of the Brγ emission

The VLTI/GRAVITY observations indicate that the Brγ-emitting region is similar to or smaller in size than the region where the hot dust resides. Figure 7 demonstrates this by means of the ratio of the Brγ to the 2.2 μm continuum sizes as a function of source luminosity. We would like to put this result into context.

Heating and ionisation by shocks and radiation occurs in the gaseous structures that make up the accretion environment in young stars. The Brγ transition is therefore a prime diagnostic in young, embedded stars and is often resolved in spectro-interferometric observations at 100 m baselines. In accreting, low-mass stars, the Brγ transition is compact and traces the magnetospheric accretion columns (GRAVITY Collaboration 2020b; Bouvier et al. 2020). In this case, the emission is located well within the dust sublimation radius of the disc. Stepping up in mass however, the magnetosphere becomes less important as field strengths decrease while at the same time the star is hotter. As a result, the Brγ emission in Herbig Ae stars is no longer restricted to the magnetosphere but is observed to be more extended, albeit still smaller than the dust sublimation radius (Garcia Lopez et al. 2015; Caratti o Garatti et al. 2015; Mendigutía et al. 2015; Ellerbroek et al. 2015; Kurosawa et al. 2016). Whether the emission is restricted to the protoplanetary disc of the HAe star or subtends a larger angle remains unclear. In some late B-type PMS stars, where the line emission region is smaller than the dust continuum, a disc-wind originating in the gaseous parts of the inner (0.2 au) disc is favoured (Kreplin et al. 2018). In mid-B type stars, Brγ and the dust occupy a similar emission region (Hone et al. 2019).

A break with the ‘Brγ smaller than dust continuum’ trend is found in the early B-type Herbig Be star MWC 297 (17 M_⊙ Vioque et al. 2018), for which the ionised emission was found to be 40% larger than the hot dust continuum (Malbet et al. 2007), while the kinematics are consistent with that of a disc wind at scales of a few au (Hone et al. 2017; Weigelt et al. 2011). For the B 1.5IV star MWC 297, having settled on the ZAMS, the generation and ionisation of a disc wind extending beyond the disc’s hot dust is perhaps not surprising. On the other hand, what could be considered surprising is the observed situation in MYSOs, where the ionised emission is systematically more compact than the dust emission.

MYSO examples for which Brγ could be spectro-astrometically mapped by means of Integral Field Unit observations and closure phases indicate a bipolar geometry at high (500 km/s) velocities. Both W33A VLA1 (Davies et al. 2010) and IRAS13481-6124 (Caratti o Garatti et al. 2016) demonstrate the origin of Brγ in fast, collimated jets and/or in collimated winds at the base of the jets. This picture is extended to deeper embedded sources, where fast ionised jets can be mapped in the radio, as done in the MYSO sources Ceph HW2, GGD 27, and G345.4938 (Curiel et al. 2006; Masqué et al. 2015; Guzmán et al. 2016). High shock velocities would likely destroy any molecule within, contrasting this outflow component from any molecular disc emission. Notably, at high accretion rates (M_acc > = 10⁻³ M_⊙ yr⁻¹), the accreting MYSO is expected to be bloated (Hosokawa et al. 2010) and therefore cool, which can arguably prevent direct ionisation of the disc. However, Simon et al. (1983) showed that winds of MYSOs can be so dense that hydrogen is collisionally excited to its n = 2 state, whichmakes its ionisation from cooler stars possible (see also, Koumpia et al. 2020; Drew et al. 1998; Drew 1998). Indeed, a bloated star could also explain the narrow single line profiles of Brγ observed in MYSOs in various studies (e.g. Pomohaci et al. 2017), which is in contrast to the relatively broader and double-peaked lines predicted by disc models around hot main sequence stars (Sim et al. 2005). Indeed, a slightly bloated object would have a larger disc (inner) radius and thus lower rotational velocities (i.e. narrow line profiles).

Compact Brγ emission such as that discussed here, can also be the result of shocks from a jet or a disc-wind. Exceptionally, disc instabilities that allow accretion to proceed will allow the formation of a very compact gaseous disc which can be shock-ionised close to the star (G345.4938, Guzmán et al. 2020). Jets were recently traced via radio thermal emission and found to be abundant in high-mass star formation (up to 84%, Purser et al. 2021). In our sample of MYSOs, we find that Brγ emission originates from a smaller area, but co-planar to that of the continuum, and therefore if a jet is the underlying mechanism, we most likely trace the base of the jet.

We surmise that ionised disc emission does not require the star to become hot enough for example to settle on the main sequence once the accretion rate goes down. Brγ being significantly smaller than the hot dust continuum (approximately the sublimation radius) could also trace jet emission in MYSOs, or rather the base of the jet in a magneto-centrifugal disc-wind. The relatively small size (~3–10 au) of the Brγ-emitting region measured here roughly matches the MYSO jet collimation region, which is usually located at the Alfven radius, from a few to several tens of au from the source (see e.g. Fig. 13 in Kölligan & Kuiper 2018; Staff et al. 2019), depending on the stellar mass and age. Although the full details are not yet known, with our GRAVITY findings we can constrain the geometry so that axi-symmetric models are to be favoured.

Fig. 7 Measured ratio of the Brγ-emitting size (RBrγ) over the size of the 2.2 μm continuum emission (Rcont) plotted as a function of luminosity. The Brγ emission originates from a smaller region compared to the dust continuum (Table 3). The different colours correspond to the different adopted brightness distribution: Gaussian (blue), ring (red), or disc (black).

5.3 Binarity

In the current study we investigated binarity in a sample of six MYSOs (8.6–15.4 M_⊙) targeting separations of between ~2 au and 300 au (assuming a distance of 3.4 kpc) and find that one object out of the six can be better modelled with a binary geometric model (binary fraction: 17${}_{-17}^{+21}$% with a 70% confidence interval). Here, we compare our results to those observed for a wide range of mass, evolutionary status, and targeted separation, and to those predicted by theory (for a thorough review on multiplicity, see Duchêne & Kraus 2013).

5.3.1 Observed binary fractions

Our aim here is to compare our findings for MYSOs with similar parameters for different classes of objects. To do so it is important to look at the statistics tracing a similar range of separations.

Kraus et al. (2011) targeted a sample of low-mass YSOs (0.25 M_⊙–2.5 M_⊙) at separations of 3–5000 au and reported a binary fraction of between 63% and 76% (decreasing with mass). This fraction drops to 55% at separations comparable to those of our study (up to ~300 au). Connelley et al. (2008) investigated binarity in embedded Class I solar-type objects targeting separations of between 50 and 4500 au and reported a binary fraction as high as 55%, dropping down to 15% for more evolved Class I objects (for classification of low-mass protostars, see Lada 1987). Once again, if we focus on separations up to ~300 au, the Class I binary fraction drops to ~13%. When one moves to the main sequence solar-type stars, Raghavan et al. (2010) reported a binary fraction of 44%, which drops down to 28% for separations of between ~2 and 300 au. For further statistics on low-mass Class 0 protostars see Chen et al. (2013) and Tobin et al. (2016).

Moving on to more massive objects, Baines et al. (2006), studied the intermediate-mass Herbig Ae/Be stars reporting a binary fraction of ~70% at 50–750 au separations, and an increasing binary fraction with increasing mass. The non-coverage of the 2–50 au separations prevents a direct comparison with our sample, but at first instance the binary fraction is at least two times higher than what we report. Pomohaci et al. (2019) targeted a sample of MYSOs at wide separations (from 600 au up to ten of thousands of au) and report a multiplicity fraction of 31%. Although the traced separations are two orders of magnitude larger than those of our study, the reported statistics are similar and within the errors. GRAVITY Collaboration (2018) studied 16 massive young stars (Trapezium), and report a decrease of companions (<30%) at separations of between 1 and 100 au, which is in alignment with our findings (for spectroscopic close massive binaries this fraction is as low as 12%, Apai et al. 2007). We note that the multiplicity of the Trapezium which is probably very dynamically evolved (Allison et al. 2009), which may explain the low fraction of binaries at 1–100 au separations.

Lastly, massive main sequence stars (O- type) are reported to have a fraction of 53% at 2–200 au separations (Sana et al. 2014), which is more than a factor of two higher than what we find. More recently, Frost et al. (in prep.) investigate binarity towards a sample of ~40 B-type MS stars using PIONIER on the VLTI (H-band observations) and find a very high binary fraction at ~2–180 mas separations (sub-au to a few hundred au).

When examining the binary fractions at face value for similar ranges of separation (Fig. 8, Table 4), the MYSO binary fraction is at least a factor of two to three lower than the low-mass T Tauris and high-mass main sequence O stars, but is similar to what is found towards less evolved Class I objects and wider MYSO binaries. We note that the large difference in stellar mass between Class I low-mass objects and MYSOs introduces limitations regarding the natural separations of the binary components, and therefore exploring the same ranges of separation may point to different stages of dynamical processes in the evolution ofthose objects. Massive stars are more commonly found in binaries with separations of up to a few hundred au when they are in a more evolved main sequence phase (~53%) compared to their forming and young stages (<30%). Therefore, at first glance, the observed statistics suggest an increase in the massive binary fraction with evolution, which contradicts both observational findings for low-mass objects and theoretical predictions (e.g. Reipurth et al. 2014).

We investigate this finding further, and take a closer look at the specific observations and techniques used in each study. The direct comparison of various samples and methods is limited mainly because the statistics are based on non-uniform observations and techniques. As also demonstrated in Table 4, the targeted studies are characterised by different sample sizes, mass ranges, and sensitivity limits. The present study in particular is based on a sample of only six objects, and as a consequencethe associated statistics suffer from significant uncertainty. In addition, our observations do not provide a uniform threshold for companion detection within the interferometric field of view. The limited uv-coverage combined tothe specific S/N for each source produces a rather separation-dependent contrast threshold for detection (for quantitative studies on those effects see, Absil et al. 2011; Davies et al. 2018). We note that the observed uv-coverage alone could result in missing ~50% of companions at the smallest scales (~0.5–4 mas).

The studies based on K-band observations are similar in detection sensitivity (ΔK ~ 4–5) with the exception of the study of solar-type MS stars (ΔK ~ 2.5), which can partially explain the observed drop in their binary fraction as they move from the pre-main sequence phase. In Table 4 we can see that although the sample presented by Sana et al. (2014) targets very similar separations to those of our study, we trace a mass range which goes up to the lowest limit of the mass traced in that study. Taking this into consideration is more sensible when comparing our observed MYSO binary fraction with the fraction observed for MS B-type stars which cover similar mass ranges to our study. We note that the separations traced by Frost et al. (in prep.) cover our separation ranges. In addition, both studies on MS OB stars are characterised by similar sensitivity in H-band, while our study uses K-band observations. The different ages and filters between the current study and those studies probably suggest that the instrumental sensitivity corresponds to a different range of physical masses of the candidate companions. Even after one takes all these limitations into consideration, it is still difficult to attribute the striking difference between the very high binary fraction of MS B-type stars and the low binary fraction of MYSOs (17%) on observational biases alone. Lastly, Oudmaijer & Parr (2010) studied massive stars in MS, and in particular a sample of Be stars and ‘normal’ B stars in K-band, targeting ~30–2400 au separations. The binary fraction of the combined sample is ~30% (29 ± 5).

In summary, to be able to provide a confidence level on the reported fractions, and confirm or discredit the observed increase in binary fraction with age, it is necessary to survey homogeneous samples using similar observational techniques and targeting stars at similar distances. Future studies should adopt an approach similar to the one presented in Fig. 8, but should make use of large samples of observations that (i) strictly trace the same primary mass range in each bin and (ii) the same companion mass ranges sampled in each bin, and (iii) are of comparable dynamical environment (i.e. no, or the same, external influence on multiplicity).

When we turn our focus to the binary fraction of MYSOs in the K-band at different scales but similar sensitivity ΔK (this study and Pomohaci et al. 2019), we observe that the binary fraction is about the same (within the errors, 17${}_{-17}^{+21}$% to 31 ± 8%) and can be probably considered flat within the entire separation range of 2–10 000 au. We note that this comparison is among the most robust we have.

Fig. 8 Binary statistics (see Table 4) for companions located between 2 au and 300 au for low- and high-mass stars of various evolutionary stages: (i) embedded: low-mass Class I vs. MYSOs (8–16 M⊙), (ii) PMS: YSOs (0.25–2.5 M⊙) vs. young OBs, (iii) MS: OBs vs. solar-type.

5.3.2 Comparison with theory

Some of the most prominent theories on binary formation are those of the core accretion and fragmentation, disc fragmentation, and stellar migration or capture. Numerical simulations of core collapse fragmentation can predict multiples with separations of several hundred au (Myers et al. 2013). Meyer et al. (2018) predicts tighter binaries with unequal components and early-stage separations of a couple of hundred au and down to less than tens of au as the system evolves through accretion disc fragmentation. Such systems are tighter than what disc fragmentation theories previously predicted (100 au–1000 au; Kratter & Matzner 2006). Tighter binaries can also be the result of capture (Bonnell & Bate 2005) or magnetic braking during accretion (Lund & Bonnell 2018), or external stellar interactions (Bate et al. 2002); we note that most of these theories start with separations of several hundred au before evolving to tighter systems. Recently, Ramírez-Tannus et al. (2021) presented evidence for an inward migration of stars as a function of cluster age, indicating that massive binaries start their lives in wide pairs (~100 au) before evolving into tighter massive binaries after about 1.5 Myr. We note that this timescale is more than an order of magnitude greater than the age of PDS 37 (~0.06 Myr; Vioque et al. 2018). Sana et al. (2017) favoured such a hardening scenario to explain the lack of close companions at birth (~12%), starting off with pairs at >0.5 au separations. Observational studies of massive binaries in the PMS phase are needed to verify, distinguish and inform the several formation theories in place (e.g. Moe & Di Stefano 2017).

Taking all this into consideration, our findings in terms of separations (<100 au) appear to be more consistent with the theoretical predictions of disc fragmentation, which, as it evolves, favours a system with a main high-mass component and accreting low-mass companions at tight separations (Meyer et al. 2018). We note that the inclination of the binary and the individual components could not be extracted from the geometric modelling alone and projection effects cannot be excluded (i.e. the physical separation of the system may be larger). Knowledge of the orbital inclination of such systems andof their circumstellar or circumbinary discs is extremely valuable for constraining and distinguishing among the proposed theories. In particular, disc fragmentation predicts mostly coplanarity among the orbital plane and the discs of the stellar systems (Kratter & Matzner 2006), which is also what was observed towards Herbig stars (Wheelwright et al. 2011). On the other hand, binary formation via capture predicts random orientations between discs and orbital plane.

6 Summary

We present the first interferometric survey in K-band of six massive YSOs. Our study increases the MYSOs with K-band interferometric measurements by a factor of four. Below, we summarise our findings on the characteristic sizes of the hot dust at 2.2 μm and ionised gas (Brγ), and on high-mass binarity at milli-arcsecond scales using spatial information:

• We spatially resolve the crucial star–disc interface in a sample of MYSOs in K-band, and finally confirm the prominence of au-scale discs in high-mass star formation with observations;

• The K-band continuum emission is spatially resolved for all MYSOs in our sample. The 2.2 μm measured characteristic size of MYSOs shows a large scatter for the given range of luminosities, but is consistent overall with the location of the inner rim (i.e. dust sublimation radius) of a disc with an optically thin cavity;

• When the inner sizes of MYSOs are compared to those of lower mass Herbig AeBe and T Tauri stars, they seem to follow a universal trend at which the sizes scale with the square root of the stellar luminosity. Such a trend indicates that similar radiative processes take place at inner regions of young stars independently of their mass;

• The measured continuum and Brγ visibilities of G034.8211 could only be fitted with a flattened or elongated disc geometry, indicative of a close to edge-on geometry. For the remaining MYSOs in our sample we find no signatures of a flattened brightness distribution in their visibility curves;

• We find that the Brγ emission is comparable to or more compact in size than the thermal emitting dusty region, and the two emissions are spatially aligned. This new finding gives credence to disc wind and disc accretion models and their ability to describe the geometry of the inner parts of MYSOs, which appear to be prominent in massive star formation;

• G282.2988 is the only MYSO in our sample of six that required a binary geometry to fit the interferometric observables. Therefore, we report a MYSO binary fraction of 17 ± 15% in the K-band at the traced scales (few au to a few hundred au). This fraction is comparable to what was previously reported at 600–10 000 au scales, indicating a flat fraction for a wide range of separations;

• Based on the present statistics, MYSO binaries at 2–300 au separations are less common than massive main sequence stars on similar scales. This finding contradicts the observational findings towards low-mass stars and the theoretical predictions of a decrease in multiplicity with evolutionary stage. However, our findings are not free from observational biases and further investigation is needed.

In this paper, we adopt a simple approach to understanding the near-IR size distribution of MYSOs with respect to other classes of lower mass YSOs for the given dust destruction predicted by three basic disc models. We find that MYSOs follow similar behaviour to that of low-luminosity T Tauri, Herbig Ae, and most of the Herbig Be stars. In particular, the 2.2 μm size of MYSOs can be directly related to the dust destruction radius predicted by an optically thin scenario where the inner rim is directly heated by the central star. A more detailed physical modelling, where the interferometric disc sizes of this sample of MYSOs are fitted simultaneously with their SEDs and near-IR spectra (covering the Brγ, NaI and COemission) is necessary to form a more detailed picture on the innermost environment of those enigmatic objects, and will be the aim of future studies.

The present study is the first attempt to address the multiplicity of MYSOs on milli-arcsecond scales, and in particular via direct spatial measurements. Multi-wavelength observations (e.g. PIONIER; H-band, MATISSE; M-, N-bands on VLTI) of large sample of MYSOs would be a great asset, providing binary statistics and helping to fully identify and characterise companions.

Acknowledgements

We would like to thank the anonymous referee for providing helpful comments and suggestions that improved the paper. E.K. is funded by the STFC (ST/P00041X/1). J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1. A.C.G. has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement no. 743029). I.M. acknowledges the Government of Comunidad Autónoma de Madrid (Spain) for funding his research through a ‘Talento’ Fellowship (2016-T1/TIC-1890). We would like to thank Jacques Kluska for stimulating discussions when preparing the paper. This publication is based on observations collected at the European Southern Observatory under ESO programme(s) 0102.C-0838, 0103.C-0459 (GRAVITY), and 092.C-0064 (AMBER). This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France. The GRAVITY data reduction was undertaken using ARC3, part of the High Performance Computing facilities at the University of Leeds, UK. This research has made use of the AMBER data reduction package of the Jean-Marie Mariotti Center (http://www.jmmc.fr/amberdrs).

Appendix A Observations

Fig. A.1 uv-plane coverage of VLTI/GRAVITY of G282.2988, G287.3716, G301.8147, and G034.8211. North is taken as position angle PA=0° and east as PA=90°.

Fig. A.2 uv-plane coverage of VLTI/AMBER of G233.8306, and G231.7986. North is taken as position angle PA=0° and East as PA=90°.

Fig. A.3 Relative flux, visibility, differential and closure phase as a function of wavelength around the Brγ emission (and continuum around it) towards G034 using GRAVITY on the UTs. The baseline length, and the position angle for all observed six baselines and four triplets of the UT configuration are also given.

Fig. A.4 Same as Fig. A.3 but for G282.

Fig. A.5 Same as Fig. A.3 but for G301.

Fig. A.6 Same as Fig. A.3 but for G287.

Appendix B Image reconstruction

In addition to the simple geometric modelling of the visibilities, we investigate asymmetries (making use of the closure phases) of the brightness distribution of the 2.2 μm continuum emission. To this end we perform model-independent image reconstruction towards a sample of MYSOs for the first time, and in particulartowards the sources observed with GRAVITY. The quality of the image reconstruction increases with: (i) the ratio between the emitting region and the angular resolution, (ii) the brightness of the source, and (iii) the coverage of the uv-plane. In our case, we are limited by all three factors (e.g. Figure A.1), and therefore the resulting images are highly influenced by the beam shape of the observations and cannot trace detailed structures of the emission. Nevertheless the image-reconstruction algorithms make it possible to independently measure the size of the 2.2 μm continuum emission and to make use of the closure phases to trace possible asymmetries, and to investigate structures that simple geometrical model fitting may have missed (See Sect. 3.1, Table 3).

Figure B.1 presents the resulting image reconstruction using WISARD (Meimon et al. 2004, 2008), which is implemented in JMMC OImaging tool⁵ towards G282.2988, G287.3716, G301.8147, and G034.8211. For direct comparison with the geometric modelling, we present the measured sizes of a 2D Gaussian distribution of the reconstructed images in Table 3. The measured reconstructed sizes of G034, G301, and G287 differ by 0.2%, 3%, and 8% respectively, compared to the sizes obtained by geometric modelling after applying a Gaussian brightness distribution. This difference is up to 27% for G282. We attribute the observed asymmetries towards G287.3716 and G301.8147 to image reconstruction artefacts, because they cannot be justified by the small observed closure phases of the sources.

One source of our sample, G034.8211, appears to be the only source with a ratio between its size and the angular resolution higher than 1. The measured ratio of G034.8211 is 1.1, while for the other three sources is between 0.67 and 0.93. Therefore, G034.8211 makes a strong case to investigate the image reconstruction further. Wanting to limit the algorithm dependency on the reconstructed images, we proceed to a thorough comparison of WISARD with the Multi-aperture ImageReconstruction Algorithm (MIRA; Thiébaut 2008) and BSMEM (Baron & Young 2008). A detailed description of the initialisation of the image reconstruction algorithms is presented in Koumpia et al. (2020).

As we see in Figure B.1 all three algorithms produce very similar elongated disc structures and show a striking similarity with what the geometric modelling revealed for this MYSO. We note that the reconstructed image of G034 using the BSMEM algorithm reveals some asymmetric emission within the inner ~3mas, in the form of a ‘midplane shadow’. The observed asymmetry may depict the observed changes in closure phases of the continuum emission towards the source (up to 20 degrees at the smallest scales; Table 2). We note that similar midplane shadows have been observed at larger scales (hundreds of au) towards a number of low mass T Tauris sources in the near-IR (e.g. IM Lup, RXJ 1615, MY Lup, DoAr 25 Avenhaus et al. 2018; Garufi et al. 2020), and have been attributed to opacity effects and scattered stellar light. Lastly, such asymmetries can be also attributed to an unresolved binary component at <10 au separations, or ongoing fragmentation being the cause of an observed asymmetry in the brightness distribution.

In conclusion, the image reconstruction results in measured sizes which are consistent within 10% (for three out of four sources) with the sizes resulting from the geometric modelling presented in Sect. 3 after applying a Gaussian brightness distribution.

Fig. B.1 Top: Image reconstruction of the 2.2 μm continuum towards G282.2988, G287.3716, G301.8147, and G034.8211 using WISARD. Bottom: Image reconstruction of the 2.2 μm emission towards G034.8211 using three different algorithms (WISARD, BSMEM, and MIRA). In addition, we present the image resulted independently by applying a simple geometric model (bottom left). The white contours represent the reconstructed image of the calibrator star for evaluation of the achieved spatial resolution.

Appendix C Near-IR spectra

Figure C.1 presents the near-IR spectra of G287, G231, G233, G282, and G301. These spectra were obtained using IRIS-2 on the AAT between 2006 and 2008. The spectral resolution in the H and K bands was ~2400.

Fig. C.1 Near-IR (left column: H-band, right column: K-band) spectra of G233, G282, G287, and G301 (The spectra of G034 and G231 are published in Cooper 2013; Cooper et al. 2013))

References

All Tables

All Figures

Fig. 1 Normalised spectra of the MYSOs observed with GRAVITY around the 2.2 μm continuum level. The wavelength coverage contains the Brγ which is seen in emission towards all sources, and the CO bandheads which are not detected towards G301. G301 also shows some weak emission within the 2.05–2.1μm range (Fe II, He I). The narrow absorption features are spectral artefacts caused by telluric correction.

In the text

Fig. 2 Modelled visibilities of the 2.2 μm continuum emission as function of baseline assuming a single disc brightness distribution overplotted with the observed ones as obtained with GRAVITY for (a) G034.8211, (b) G301.8147, and (c) G287.3816 as fitted using LITpro (Sect. 3.1). (d, e) Same as before but for observations obtained with AMBER towards G233 and G231, after fitting a Gaussian, a uniform disc, and a ring. The background flux contribution is found to be 16% for G233 and ~9% for G231.

In the text

	Fig. 3 Modelled visibilities and closure phases of the 2.2 μm continuum after assuming a flattened disc, overplotted with those observed towards G034.8211. The observed visibilities of the Brγ emission arealso overplotted in red, demonstrating the co-planar origin of the continuum and the ionised gas.
In the text

	Fig. 4 Modelimage of the 2.2 μm continuum emission towards M8EIR (based on Frost et al. 2021) convolved at 1.7 mas resolution (GRAVITY/VLTI). The distance to M8EIR is 1.3 kpc.
In the text

Fig. 5 (a) Near-IR size–luminosity diagram of our observed sample of MYSOs. The measured inner ring size Rin of the observed MYSOs (GRAVITY, AMBER; based on Table 3) is plotted with respect to the stellar luminosity. The shaded area represents the dust sublimation radius Rs as predicted by the presence of a disc with an optically thin cavity for a range of temperatures. The adopted dust sublimation temperatures are 1000 K (solid line) and 1500 K (dotted line). (b) Same as before but overplotted with two MYSOs with K-band measurements from the literature (filled black stars; IRAS 13481-6124, NGC 2024 IRS2, Kraus et al. 2010; GRAVITY Collaboration 2020a) and with indirect measurements (star symbols; modelled sizes presented in Frost et al. 2021). The open and filled stars represent the sources for which the modelled inner radius was found to be larger than or equal to Rs respectively. The dust sublimation radius as predicted by models of different dust temperatures with (blue area) and without (yellow area) backwarming effects is also overplotted. (c) Same as before but including the available measurements of Herbig AeBes and T Tauris from the literature. The dust sublimation radius as predicted by the classical scenario of a flat optically thick disc is also overplotted.

In the text

	Fig. 6 Modelled visibilities and closure phases of G282.2988 as fitted using a binary model in LITpro. The best fit reveals a component separation of 23.7 mas (consistent with Koumpia et al. 2019) and primary to secondary flux ratio of ~10.
In the text

	Fig. 7 Measured ratio of the Brγ-emitting size (RBrγ) over the size of the 2.2 μm continuum emission (Rcont) plotted as a function of luminosity. The Brγ emission originates from a smaller region compared to the dust continuum (Table 3). The different colours correspond to the different adopted brightness distribution: Gaussian (blue), ring (red), or disc (black).
In the text

	Fig. 8 Binary statistics (see Table 4) for companions located between 2 au and 300 au for low- and high-mass stars of various evolutionary stages: (i) embedded: low-mass Class I vs. MYSOs (8–16 M⊙), (ii) PMS: YSOs (0.25–2.5 M⊙) vs. young OBs, (iii) MS: OBs vs. solar-type.
In the text

	Fig. A.1 uv-plane coverage of VLTI/GRAVITY of G282.2988, G287.3716, G301.8147, and G034.8211. North is taken as position angle PA=0° and east as PA=90°.
In the text

	Fig. A.2 uv-plane coverage of VLTI/AMBER of G233.8306, and G231.7986. North is taken as position angle PA=0° and East as PA=90°.
In the text

	Fig. A.3 Relative flux, visibility, differential and closure phase as a function of wavelength around the Brγ emission (and continuum around it) towards G034 using GRAVITY on the UTs. The baseline length, and the position angle for all observed six baselines and four triplets of the UT configuration are also given.
In the text

	Fig. A.4 Same as Fig. A.3 but for G282.
In the text

	Fig. A.5 Same as Fig. A.3 but for G301.
In the text

	Fig. A.6 Same as Fig. A.3 but for G287.
In the text

Fig. B.1 Top: Image reconstruction of the 2.2 μm continuum towards G282.2988, G287.3716, G301.8147, and G034.8211 using WISARD. Bottom: Image reconstruction of the 2.2 μm emission towards G034.8211 using three different algorithms (WISARD, BSMEM, and MIRA). In addition, we present the image resulted independently by applying a simple geometric model (bottom left). The white contours represent the reconstructed image of the calibrator star for evaluation of the achieved spatial resolution.

In the text

	Fig. C.1 Near-IR (left column: H-band, right column: K-band) spectra of G233, G282, G287, and G301 (The spectra of G034 and G231 are published in Cooper 2013; Cooper et al. 2013))
In the text