© The Authors 2024

1. Introduction

Ever since it has been unanimously established that a significant fraction of stars, regardless of their mass, evolve in pairs or multiples, the study of star formation has taken a major turning point. Yet, there is not a consensus on the physical processes and mechanisms that settle multiplicity, and it is of prime importance to include multiple formation into star formation studies, in order to comprehensively understand star formation as a whole. The studies of young populations of low-mass stars aim to quantify the intrinsic physical properties of multiple star systems and the distribution of binary parameters accurately, and to unravel the complex processes and history of binary formation. The presence and abundance of these low-mass (proto-)stars in our close environment (within 1 kpc), in clusters or as field stars and visible in the optical light, has made it possible to converge on reliable statistics for this mass regime (see for instance: Duchêne et al. 2007; Connelley et al. 2008; Raghavan et al. 2010; Kraus et al. 2011; Duchêne & Kraus 2013; Ward-Duong et al. 2015; Elliott et al. 2015; Tobin et al. 2016, 2022; Elliott & Bayo 2016; Moe & Di Stefano 2017; Kounkel et al. 2019; Zúñiga-Fernández et al. 2021). These studies converge towards a fraction of 50% of solar-type stars in the field with a companion, with an increased incidence in star-forming regions, indicating that many low-mass stars are born in binaries. These statistics have been derived from a wide range of observing techniques, ranging from spectroscopy, interferometry, sparse aperture masking, and adaptive optics (AO) assisted imaging, and thus they provide a comprehensive view of binarity from the innermost to the outermost environments (from 0.1 au to a couple dozen 10⁴ au).

The question of the formation of massive binaries that populate the upper part of the Hertzprung-Russell (HR) diagram (M_i, primary > 8 M_⊙) is all the more crucial as massive stars are almost entirely observed while on or already off the zero-age main sequence (MS) phase, as a result of their rapid evolution. As such, most surveys about multiplicity refer to fully formed MS O and B stars, for which nearly 100% are known to belong to a multiple system with a high tendency towards close (≲1 au) orbits (Sana & Evans 2011; Sana et al. 2012, 2013, 2014; Dunstall et al. 2015; Moe & Di Stefano 2017). Despite the observational obstacles to image young massive stars, teams have combined efforts to address binarity in the younger stages of massive stars (Apai et al. 2007; Pomohaci et al. 2019; Koumpia et al. 2019, 2021; Shenton et al. 2023). Binarity in the high-mass regime is known to significantly impact the evolution and fate of the systems, and as such constitutes an important axis of study for clarifying the origin of gravitational waves and understanding the dynamical and chemical structure of our Galaxy. Nonetheless, even with multiplicity being a ubiquitous feature of massive stars, binarity rates and orbits derived from MS populations do not necessarily reflect the primordial binary properties as the latter is impacted by secular evolution and dynamical processes (see e.g. the review of Offner et al. 2023). An important point to address concerns the physical processes and mechanisms that settle multiplicity.

Numerous theories try to explain how binaries and multiples are formed. They include disk fragmentation, core fragmentation, and capture in the protostellar phase (Tohline 2002; Kratter et al. 2008; Krumholz et al. 2009; Meyer et al. 2018). These proposed pathways do not solely apply to massive binaries, but rather explain the possible pairing mechanisms regardless the mass of the primary. Dynamical capture occurs when, in dense stellar environments where close interactions are frequent, an object initially born single and unbound to any other system may capture passing protostars to form a multiple system (Tohline 2002). This process is expected to take place later on, during the dissolution of the star-forming region (Fujii & Portegies Zwart 2011; Parker & Meyer 2014).

The two most prominent theories involve fragmentation within the gravitationally collapsing parent cloud. However, the timescales and scales at which fragmentation occurs is still very uncertain. Fragmentation of a turbulent pre-stellar core occurs when density inhomogeneities caused by turbulence create clumps above the local Jeans mass that collapse faster than the background main core, ultimately forming a multiple system with separations 100 au ≲ a ≲ 10 000 au (Goodwin et al. 2004; Offner et al. 2010; Myers et al. 2013). In disk fragmentation theories, multiplicity is presumably generated by an accretion disk that is subject to strong gravitational instabilities causing subsequent fragmentation (Bonnell & Bate 1994). Numerical simulations involving this latter process have been efficient in forming companions around massive protostars (Kratter & Matzner 2006; Krumholz & Bonnell 2007; Oliva & Kuiper 2020; Mignon-Risse et al. 2023). The observations of fragmenting accretion disks suggest that this mechanism tends to form rather close systems with separations a ≲ 100 au, and provides an ideal framework for inward migration that may result in spectroscopic binaries (Sana et al. 2017; Meyer et al. 2018).

In order to gain a better understanding of the formation of massive multiples as well as to predict their future evolution, we must have a clear description of the multiplicity properties, the mass ratio, and the initial separation distribution of forming massive stars (Sana et al. 2012). The massive young stellar objects (MYSO) phase is a short-lived evolutionary stage towards the end of the hot core phase in which the protostar experiences accretion events necessary to grow in mass. Evidence of circumstellar accretion disks is, for instance, observed in Kraus et al. (2010) and Frost et al. (2019) as well as bipolar outflows (Cooper et al. 2013) and a surrounding dusty envelope, which is opaque at visible wavelengths (Wheelwright et al. 2012). Keplerian-like disks that could be associated with accretion in larger scales are also reported in various studies (Johnston et al. 2015; Zapata et al. 2019; Maud et al. 2019). Unlike low-mass protostars, the large distances at which they are located, their high dust extinction, and rarity has made systematic studies challenging. It is only recently with the advent of high angular resolution and high-contrast observations assisted with AO that the study of multiplicity among small populations of MYSOs has been made possible, slowly lifting the veil on the confusion that reigns as to their formation. We can, for example, quote the studies of Pomohaci et al. (2019) whose work probed for wide (> 1000 au) multiple companions of 32 MYSOs, assessing a field binary fraction of 31 ± 8%, or the recent study of Koumpia et al. (2021) that reported a fraction of close binaries as low as 17% 17 ± 15% (2 mas ≲ a ≲ 300 mas) among a sample of six MYSOs. Our present study aims to complement the pilot study led by Pomohaci et al. (2019), using a smaller sample of MYSOs with a slightly longer wavelength filter (3.8 μm), ensuring the detection of even younger companions that might not have been visible in the K-band images. Expanding these studies is crucial to find a consensus on the multiplicity and the initial orbital parameters in order to determine which of the scenarios presented above is the most likely to form multiple massive systems. These tabulated values over the full stellar mass and separation range can serve to feed numerical and analytical models or as a comparison for massive star formation.

Here, we study the environment around thirteen young massive stars using mid-infrared (MIR) wavelengths, with the ultimate goal to address their multiplicity. We investigate the presence of potential companions in wide orbits with AO-supported observations at 3.8 μm. In Sect. 2, we present the observations and the reduction routine. Section 3 deals with the data analysis, ranging from source detection, point-spread function (PSF) fitting, detection limits, and chance projection probability. The results per object are described in Sect. 4. The multiplicity and companion fractions and the implications of our findings in the context of star formation are discussed in Sect. 5. We summarise the results and conclude in Sect. 6.

2. Observations and data reduction

Our target list of 13 MYSOs is taken from the Leeds Red MSX Source (RMS) survey, where MSX is the Midcourse Space Experiment telescope (Lumsden et al. 2013). This survey aimed to search for MYSOs across the Milky Way. Colour cuts were used to find the initial 2000 candidate MYSOs from the survey data, following the criteria of Lumsden et al. (2002). Follow-up observations at different wavelengths were performed to exclude other objects that satisfied the colour requirements such as planetary nebulae (PNe), UCHII regions and evolved stars to allow the robust identification of the MYSOs. The survey ultimately found ∼800 (candidate) massive young stellar objects and H II regions within 18 kpc and L > 2 × 10⁴ L_⊙ (Lumsden et al. 2013). The target selection from this database for our sample was based on bolometric luminosity, a 4 kpc distance cut, with declinations (δ < 10°) accessible to ESO’s Very Large Telescope facility.

Our observations aimed to detect putative companions to a sample of MYSOs. In order to do so, we chose to probe multiplicity in the mid-IR wavelength region at the highest possible angular resolution. This wavelength choice was motivated by minimising the effects of dust extinction and the envelope’s dust emission. Depending on the source, dust emission from the envelope dominates the spectral energy distribution (SED) from about 5 microns onwards. Hence, L′ was chosen as the optimal wavelength to be able to detect companions to the central MYSO that are up to a factor of 1000 fainter.

We selected MYSOs with L′ magnitudes comprised between 7 and 10. All 13 objects are classified as MYSOs in the RMS catalogue and a summary is provided in Table 1. Direct adaptive-optics (AO) assisted imaging is to date the most relevant method to explore the wide-orbit environments around massive protostars. To study the multiplicity and the environment of the objects listed in Table 1 therefore, we made use of the capabilities the NAOS (Rousset et al. 2003, Nasmyth Adaptive Optics System) adaptive optics system integrated with the CONICA (Lenzen et al. 2003, COudé Near IR Camera) NIR imaging camera (or NACO), an instrument decommissioned in 2019 from ESO’s facilities.

The L′-band observations were acquired in service mode over 9 different nights spread between November 26, 2012 and January 27, 2013. The seeing was variable during these nights, ranging from 0.68″ to 1.42″ in the visible, meeting the requirements for good quality products. For the purpose of our study we used the camera mode L27 with the narrow band filter Lp (3.8 μm), providing a medium resolution of ∼27 mas pix⁻¹ and resulting in a field of view (FoV) of 28″ × 28″. The pixel size of 0.027″ corresponds to 54 au given the average distance of 2 kpc of our sources. However, the survey probes the separation range from the minimum achieved full width at half-maximum (FWHM, calculated from the detected new sources) of ∼160 mas up to the full FoV of 28 × 28 arcsec. For a typical distance of ∼2 kpc, this translates in separations ranging from ∼320 au to 56 000 au. The detection limit of low-mass companions ranged up to L′ = 16, given the high extinction of the targeted MYSOs. The average L′ magnitude of the targets is 5.1 (see Table 1). With K < 10.5 mag, all the targets could be used as their own reference sources for wavefront sensing (the magnitude limit for wavefront sensing being K < 13 mag). We used the infrared wave-front sensor (IR WFS) with the JHK dichroic that passes 100% of the NIR light to the NAOS IR WFS and forwards all the L′ light to CONICA onto the science detector. The observations were obtained in short-exposure cube mode that stores each individual co-add image instead of a median-combined one. It allows for frame selection with integration times of less than 0.2 s.

In order to get rid of the background contamination and obtain a precise bad pixel map, each target was observed using a jitter sequence (a dozen images) with a fixed sky offset. This moves the telescope alternatively between object and sky positions within a predefined box of 20″ × 20″, allowing the sky from all the observations to be estimated. For each source, a total of two cycles were observed, each cycle consisting of a sequence of object and sky observations following the jitter pattern.

We reduced the data initially using the standard ESO NACO pipeline recipes (public version 4.4.11). However, the recipes were found to be unsuitable for proper estimation of the background using dedicated sky frames. We opted instead to custom-made our own reduction routine for the pre-processing of the raw images. We downloaded the processed calibration files from the ESO archive. All science and sky cubes were bias subtracted, flat-fields corrected and bad pixels masked. Given the large number of frames per observation, we performed a quick test throughout the cube (both for object and sky) to identify and remove bad frames (based on the cross-correlation of each frame to the median of all frames in a given cube). This led to the removal of ∼3% of all images on average, with mainly the first 3 frames’ quality being below average, likely due to the start of the adaptive optics loop or instrumental effects. For each of the 2 cycles, we median-combined together sets of 4 consecutive object (resp. sky) images, leading to 2 separate object (resp. sky) stacked images. Each object image was sky subtracted by its corresponding combined sky image. We aligned the two reduced images using the offset information stored in the header and the master image was obtained by stacking the two images. The reduced images are presented in Figs. A.1–A.3 and B.1.

3. Data analysis

We have imaged a total of 13 young massive stars for which we aim to identify possible companions between 300 and ∼60 000 au and characterise the thermal infrared environment. To achieve this purpose, we follow the same data analysis procedure for each system. Our analysis consists of; (1) measuring the background, (2) identifying all possible point sources and (3) determining the astrometry and photometry of these detected objects. We provide all the details about the analysis in this section and discuss the sensitivity of our observations.

3.1. Background and noise estimation

It is essential to derive an accurate estimate of the background emission level to measure the photometry of our targets. Similarly, the significance of source detection and photometric errors are affected by the background noise. We use the sigma clipping technique that removes the pixels that peak above a specified sigma level from the median and recalculates the statistics until convergence is reached. We fixed σ at 3 and the median filter is calculated over each image following the procedure written by Bradley et al. (2022). However, the background level may vary across the image. To probe background level variation across the detector, we generated for each science observation a 2D image of the background. We computed sigma-clipped statistics on a grid to create a low-resolution background image. The meshes have a size of 50 × 50 pix (1.3″ × 1.3″), which is small enough to catch any background variations and larger than the size of the sources. We obtained the final background map by interpolating the low-resolution image. The background map was subsequently subtracted from the co-added images to generate the final reduced one.

3.2. Companion detection

The identification of astronomical sources in the images is a precondition to performing photometry and astrometry measurements. In order to identify physical companions in the frames, we first need to pinpoint all possible point sources. A first visual approach and inspection of the reduced images confirm the presence of nearby and rather bright sources cohabiting with fainter point sources. We developed our own routine to numerically identify star-like objects in the FoV around each central MYSO. The code mainly uses the PYTHON affiliated package PHOTUTILS (Bradley et al. 2022) of ASTROPY that provides an implementation of the DAOFIND algorithm. The latter is designed to detect astronomical objects in crowded and non-crowded background-subtracted images. A similar approach has been used in Bodensteiner et al. (2020), for example, and has proven its efficiency. We initialised the star finder to search for objects with FWHMs of around three pixels and peak at least 5σ above the background signal (background map). We first set standards on the roundness such that only star-like objects are returned and blob-like-shaped objects are discarded. We visually examined the images and looked for any missed detection. After inspection, it seems that no object was left behind, but given the possible extended nature of MYSOs and keeping in mind that most of the companions might be young forming stars, this parameter was later on left free in order to allow the detection of any exotic shape. In practice, no additional objects were found when relaxing the roundness parameter. The returned parameters contain the results of the search, including the source locations and peak amplitude. These parameters feed the PSF fitting code that performs photometry measurements on the detected sources, which we explain in the following subsection.

3.3. PSF fitting

We follow the PSF-fitting methods as detailed in Bodensteiner et al. (2020) and Rainot et al. (2020). From the companion detection routine, we obtain accurate centroid positions for each peak detected in the image. The low source density in the images allows simple aperture photometry as a first estimation. We build a set of aperture objects (typically a circular aperture with a radius of three pixels) from the peak locations and perform the photometry taking into account the overlap of apertures on the pixel grid. The fractional overlap of the aperture is calculated and the aperture sum for each detected peak is derived accordingly. However, this method is not exact and some flux contribution outside the extraction area might be missed when summing the photons. We therefore also used a second method consisting of PSF-fitting, that is, modelling the radial shape of each star in the image. Simultaneous PSF-fitting photometry better handles stars with close neighbours, which is the case in some of our images (see Sect. 4) for the achieved spatial resolution. Indeed, some companions appear to be in such pairs. Figure 1 shows an example of an overlapping aperture originally detected as a single source without the DAOFIND method. For such systems, the PSF-fitting method derives more accurate fluxes. To generate the PSF, we use the high signal-to-noise primary star and use the derived aperture photometry as a first instance. None of the central stars appears to be resolved into a close binary. This PSF model was fit to each source individually whose positions are the fixed centroids derived from the source detection procedure (see Sect. 3.2). In this way, we obtain estimates of the total flux and flux errors.

Fig. 1. Zoom over the northern companion system for target G268.3957, showing two distinct components. The white solid lines is the 4-pix aperture applied in the aperture photometry method and the markers show the position of the centroids as found with DAOFIND. For reference, north is up and east is left.

3.4. Photometry and astrometry

The fluxes and positions derived from both methods, aperture and PSF fitting, are consistent within error bars for isolated sources, i.e. non-blended sources. In the case of close binary systems, derived centroids are unchanged and consistent with the PSF method. However, the photometry is equally balanced between the two components when modelling the brightness with a PSF. For the remainder of the paper, we adopt all fluxes and source locations derived via the PSF-fitting method.

The fitted fluxes of the newly detected sources are converted into L′ mag relative to the brightness of the primary star, for which we have an estimation of the L′ mag (see Table 1). No photometric calibrators were contained in our observations thus we derive relative photometry with respect to the central star. For every single point source that is detected, we calculate the angular separation from the primary, first in detector coordinates that we then convert into angular units, applying the pixel scale of 0.027″ pix⁻¹, as stated in the NACO user manual¹ and confirmed in the header file of the targets. The physical distances and position angles are straightforwardly computed using the distances reported in the RMS Catalogue (Lumsden et al. 2013). We additionally checked the Gaia eDR3 (Bailer-Jones et al. 2021) distances for every single source. Massive young stars are expected to become visible once they have reached or are close to the birthline. Given their strong IR emission rising from their circumstellar material, the young high-mass objects are often opaque to Gaia bands, reducing our chances to find reliable parallaxes. Among our sample of 13 MYSOs, we find reliable Gaia distances for three objects (G265.1438, G268.3957 and G318.0489) whose coordinates matched the RMS catalogue within 2 arcsec. For these three targets, the Gaia distances are consistent with the RMS distances. We use the RMS distances for the remainder of the paper.

3.5. Sensitivity limits

We determine our contrast limits following the methods of Pomohaci et al. (2019) and Rainot et al. (2020). The sensitivity of the observations is determined for each target in terms of flux ratio with respect to the central object. We computed the contrast limits by injecting 300 artificial stars into each image. The brightness varies from a thousandth to a hundredth of the central object and the simulated objects are randomly distributed across the image. Each artificial star is a 2D Gaussian source that is built following with the same FWHM as the PSF. We then run the same procedure as for the companion detection, and set the minimum flux at which the simulated source was detected at the 5-σ level above background noise as the limiting flux. The results of a single iteration are shown in Appendix C. Since the artificial sources have random brightness and are placed at random distances across the image, we repeat the process over 20 iterations and the final detection limits recorded in Table 1 are issued following the bootstrapping procedure. Figure 2 shows the detection limits as a function of the L′ mag of the primary. For each system, we display both the contrast limit measured in the real images and the one derived from the injection of artificial sources. The contrast varies from 7 to 1 depending on the L′ mag of the central star but is roughly constant in terms of the apparent companion magnitude $L_{\text{comp}}^{\prime }~8$, showing that the images are rather background limited. The noise does become indistinguishable from background noise at an angular separation of about 4″ from the central star.

Fig. 2. Detection limits for each observed field. The star symbol denotes the Δmag of the faintest source detected in the NACO image (if any) while the dot symbol refers to the Δmag of the faintest detected artificial source, injected into the image. The minimum angular separation at which companions are found is about 2.7″ which translates to ∼5000 au at a typical distance of 2 kpc. However, the faintest companion of each field is found at a typical separation of 8.7″ (∼17 000 au).

4. Results

We have imaged 13 MYSOs in the mid-IR at a spatial resolution of ∼160 mas in order to probe the presence of companions. We identified a total of eight multiple systems from 13 targets, with a total of 41 newly detected mid-IR companions candidates, summed over all fields. We test the probability of physical association of the new sources to the central object in Sect. 4.1. We also describe below each of the 13 sources together with the results of our analysis.

For four of them, we compare our findings to those of previous AO-assisted observations, published in Pomohaci et al. (2019). For a thorough comparison with their study, we quickly present the similarities and differences with regard to the sample size, specific observations and techniques used in their study. Their work aimed at probing the widest companions (with a physical separation range of 400−46 000 au) around 32 MYSOs classified in the RMS Catalogue using K-band observations with NACO. Their observations were carried out roughly three years later than the L′-band ones, and provide K-band detections and their relative binary parameters. Despite a shorter wavelength filter and narrower FoV that might have lead to a higher spatial resolution probe (120 mas to 7 arcsec), they only achieve a full coverage for distances up to 3 arcsec away from the main source, as a result of technical issues during the observations that affected the bottom left quadrant of the detector. They reach an average limiting magnitude of 14 and are able to probe sources 5.5 mag fainter than the primary source, in the K band. They infer the chance projection following a similar method but using the 2MASS point-source catalogue in the K band to derive the background density of sources.

Further archival and important notes on individual sources are provided in Appendix D. The list of detected companion candidates together with their properties are given in Table 2. We point out that the sources are arbitrarily classified and are not sorted according to their proximity to the primary star or their brightness. We simply follow the order returned to us by the PSF-fitting procedure.

4.1. Eliminating spurious associations

An important step in confirming the nature of the systems is to assess whether the detected companions share a physical connection with the central object. The analysis of common proper motion and orbital motion would unambiguously discard that a visual pair is just a chance alignment. However, with only one observation in hand, we lack sufficient temporal coverage to detect whether the companion candidates are co-moving. To overcome the lack of multi-epoch observations, we employ statistical methods to determine the likelihood of a given source being bound to the central object. We follow the procedure described in Rainot et al. (2020), Pauwels et al. (2023).

At first glance, the NACO images are not densely populated and candidate companions are not uniformly distributed across the observed field. The denser the field, the greater the chance of spurious associations. Quantitatively, we estimate the potential contamination from fore- and background sources given the target’s galactic coordinates and the limiting magnitudes of the observations. As surveys and catalogues in the L′-band do not exist, we simulate a collection of background stars by running the Besançon model of stellar population synthesis of the Galaxy (Robin et al. 2003) in the L-band, from 0 to 150 kpc and spanning a surface of 1 deg² in the direction of the object. The model predicts between 0.2 and 4 stars in a 28″ × 28″ area, depending on the magnitude properties of each source. We define the probability of spurious association (P_spur) of a given companion candidate i as the probability that, given the local number density of stars (Σ(L_i)) as bright as companion i, at least one source is found by chance within a circle of radius ρ_i (i.e. separation of companion i). For consistency, we use a Monte Carlo simulation to generate 1000 populations of stars featuring these properties (Σ(L ≤ L_i)) and we sum the iterations where at least one source is spotted at ρ ≤ ρ_i. The final P_spur is given by the result of the sum, divided by the total number of simulations.

The results are provided in the last column of Table 2. We consider P_spur = 20% as a limit to disprove any physical association (Pomohaci et al. 2019). All of the 41 detected sources but two have P_spur ≤ 8% suggesting that they most likely are real companions. Companion candidate B of G310.0135 and A of G268.3957 have respectively over 30% and 69% chance of being galactic contaminants (see boldface values in Table 2).

Among the stars around which no source is detected, we perform a statistical test to find out if the non-detections are in agreement with the P_spur estimates of the other sources with detected companions. We proceed in the same way by simulating a background stellar density and using the L′-band detection limits to calculate the probability that a source is detected in a given field. We find probabilities ranging from 5−23% to find a source in the five fields where our tools failed in detecting sources, showing the real lack of L′-bright companions around these MYSOs or that NACO is not sensitive enough to detect fainter companions.

4.2. MYSOs within a distance of 2 kpc

4.2.1. G203.3166+02.0564

G203.3166+02.0564 (hereafter G203.3166) also known as AFGL989, NGC 2264 IRS1 or Allen’s object (Allen 1972) is one of the massive members of NGC 2264, a young (Venuti et al. 2019, ∼3 Myr) star forming region located at 760 pc (Dahm et al. 2008). G203.3166 has been the subject of several observation campaigns, ranging from NIR to sub-mm, we gather the most prominent results in Appendix D (Thompson et al. 1998; Schreyer et al. 2003; de Wit et al. 2009; Grellmann et al. 2011). Six companions were discovered in 1998 thanks to K-band images (Thompson et al. 1998) and two additional ones were reported in Schreyer et al. (2003).

Our NACO L′-band images allow separations between 120 and 15 000 au to be explored for this source. We confirm the presence of these companions, tagged from A to G, and we clearly detect the binary (A and B objects, in Table 2 and Fig. A.1) in the south-west (SW) direction, that was described as the extended “object 8” and proposed to be a binary in Schreyer et al. (2003). We also detect 2 additional objects, H and I, that were not previously referenced in the literature. The geometric arrangement of companions in line, on one side of the central object is intriguing. Whether they could be knots in a jet could be further investigated upon the acquisition of new data. Thompson et al. (1998) had noticed this distribution of objects and suggested that it shows evidence of a triggered star formation. We discuss this scenario in the next section.

4.2.2. G232.6207+00.9959

G232.6207+00.9959 (hereafter G232.6207) lies at a distance of about 1.68 kpc as estimated by Reid et al. (2009) using the trigonometric parallax measurements of the 12 GHz methanol masers. The associated IRAS point source is IRAS 07299−1651 for which we find multi-wavelength studies in the literature (see Appendix D).

The image of G232.6207 shows a bright nebulosity in the L′-band surrounding the targeted source, which could be associated with the outflows detected at longer IR wavelength, for example with SOFIA (De Buizer et al. 2017). With NACO we probe separations ranging from ∼270 − 33 000 au. We detect seven sources around this young star at separation ranges between ∼8000 to ∼25 000 au, with flux ratios ranging from ∼5% to almost 70%. All of these candidate companions are statistically likely to be bound to the central object, or at the slightest fault belong to the same cluster. Pomohaci et al. (2019) report the detection of two candidate companions in the K-band with a magnitude contrast of 5.0 and 2.6 but these are likely to be merely visual binaries given their high probability of spurious associations. Being located at ∼6.4″ and 6.7″ from the central source, these detections are most likely different stars than the ones we find. Objects G and E are aligned with the direction of the jet axis. Thus, one cannot ignore the possibility that they could be knots of material that have formed due to outflow processes (e.g. Bally et al. 1995). Spectroscopic study or time-series imaging could assist in determining the definite nature of these sources, but that is outside the scope of this work.

4.2.3. G263.7434+00.1161

G263.7434+00.1161 (hereafter G263.7434) has been classified as a candidate MYSO and belongs to the Vela Molecular Ridge (VMR; Mottram et al. 2007). Urquhart et al. (2007) report that the protostar has a spherical and unresolved morphology. The 10.4 μm TIMMI2 image displayed in Mottram et al. (2007) shows a single object while the 2MASS K-band images are slightly more crowded.

With a distance of 0.7 kpc, NACO can detect companions from ∼110 to 14 000 au. No source apart from the central one is detected in the L′-band images. We found a 5σ detection limit of Δ = 6.8^m.

4.2.4. G263.7759–00.4281

G263.7759−00.4281 (hereafter G263.7759) is better known under the name IRAS 08448−4343 and is, like G263.7434, located in the VMR. The protostar has a luminosity equivalent to that of a main sequence spectral type B2 (Navarete et al. 2015) and is located at a distance of ∼0.7 kpc, which means that any companion in the sensitivity of NACO can be detected between ∼110 and 10 000 au. The resolved 20 μm images of Wheelwright et al. (2012), as well as the H₂ continuum subtracted images (Giannini et al. 2005), agree on the extended emission of the central source, showing clear signs of a structure likely to be a bipolar jet.

The extended emission is also noticeable in our L′-band images when performing the source detection code. We relaxed the roundness parameter giving more flexibility to the algorithm to detect the central source. NACO’s FoV allows the detection of companions between 110 and 14 000 au. We detect three additional objects, a binary system at ∼10″ in the SW direction and a single object at ∼3″ pointing to the north-west (NW) direction. After visual inspection of the L, M, N images of Giannini et al. (2005), the position of object 40 seems to coincide with that of the binary. Object 40 is a sub-complex of six resolved peaks in the K-band, characterised as embedded protostars, and proposed to be the driving sources of the outflow. Such a sub-cluster is not resolved in our images, apart from the binary-like system. Nonetheless, none of the two sources cited above report the detection of the NW object. Given its direction that seems aligned with the H₂ jet detected by Giannini et al. (2005), we cannot rule out the hypothesis that this object is a knot of ejected material from the central object.

We point out the extended shape of our detected object C, closer to the primary than the double detection AB. The peak emission is broader than the PSF and can be of a different origin. Most MYSOs are expected to be extended and all the sources in the sample are classified as MYSOs. Being even more extended than the others would imply that this companion is particularly embedded, hence of a younger age. Given the relatively low spurious association probability, it may also be a young and massive forming star whose central object has brightened to the point where its surrounding material (likely envelope and disk) is better illuminated than the other sources. In the latter case, such radiation could be associated with accretion. Finally, object C could also be fortuitously inclined, producing a broader emission, which could also be a distance effect since it is particularly close. However, it seems that other equally close sources (e.g. G265.1438) do not show any extended geometry beyond the PSF.

4.2.5. G265.1438+01.4548

G265.1438+01.4548 (hereafter G265.1438) is a MYSO and member of the massive star-forming region RCW 36 located at 0.7 kpc, in Cloud C of the VMR (Ellerbroek et al. 2013). Unlike the rather void 20 μm images of the central source (Wheelwright et al. 2012), we find no fewer than eleven objects in the FoV (110 − 14 000 au), making G265.1438 the source for which the most sources are found. The sources are located along the NW to SE direction, usually above ∼10 mas around the central source, with the exception of sources G and F which are located about 7 mas in the orthogonal direction (see Appendix A). We find source J, which is slightly brighter than the MYSO, at the edge of the FoV. Chance alignment is low (< 2.3%) so one could rightfully assume the physical association with the central source. However, looking back at Ellerbroek et al. (2013)’s images, the region seems rather crowded, with a lot of Class 0/I/II objects. The latter authors and Wheelwright et al. (2012) report the presence of an outflow. No more details about the direction of the outflow can be found in the literature but the arrangement of the detected objects in the L′-band images can be indicative of a NW to SE direction. Constraining the origin of the outflow, the true nature of the detected objects (real stars or knots of ejected material) and their physical association would require more observations.

4.2.6. G268.3957–00.4842

G268.3957−00.4842 (hereafter G268.3957) is an MYSO candidate located at approximately ∼0.7 kpc, a rather short distance allowing for source detection from 110 − 14 000 au with NACO. The NACO L-band image shows the presence of five companion candidates, including a relatively NIR bright binary system (resolved as Objects D and E) in the north vicinity of the central object. All but one (object A in Table 2, which we discard in the remainder of the discussion) have spurious associations lower than 1.1%, pointing to the fact that they are most likely physical companions. The companions have projected separations ranging from ∼1350 au to ∼5800 au for the binary system. The projected separation of the two components forming the binary amounts to 0.15″ which translates to ∼107 au. All companions have also been detected in the NACO K-band study of Pomohaci et al. (2019) with angular separations matching our findings with a difference of less than 1% for objects B and C and of about 7% for the components of the binary (objects D and E). However, they find a chance projection of the binary greater than 48% and conclude that the system is rather a fore- or back-ground multiple which is inconsistent with our derived probabilities of less than 1%.

4.2.7. G269.1586–01.1383A

G269.1586−01.1383A (hereafter G269.1586) is classified as an MYSO candidate and associated with an H II region, approximately 7″ to the north-east (NE; Urquhart et al. 2007). Located at a distance of 0.7 kpc, NACO is sensitive to source detection ranging from 110 − 14 000 au. We detect two additional relatively bright objects, both at around ∼8″ (which translates to nearly 6000 au) from the protostar and with low spurious association probabilities. Both objects are respectively located towards the NE and NW directions from the central protostar, as shown on the Fig. A.1. The newly detected sources are approximately ∼32% and ∼77% as bright as the primary.

4.3. MYSOs above a distance of 2 kpc

4.3.1. G194.9349–01.2224

G194.9349−01.2224 (hereafter G194.9349, also known as RAFGL 5186) is an MYSO located 2.0 kpc away (Lumsden et al. 2013). G194.9349 has been part of NIR and FIR spectroscopic campaigns (see Appendix D for further notes Cooper et al. 2013; Cunningham et al. 2018) but we could not find NIR nor MIR imaging archives in the literature. We did not find indications for a companion star in our L′-band images and the probed separations of 320 − 40 000 au (see Appendix A). We derive a shallow detection limit of Δ = 2.0^m.

4.3.2. G254.0548–00.096108

G254.0548−00.096108 (hereafter G254.0548) has also been classified as an MYSO located at a distance of 2.75 kpc. This protostar is among the faintest of our sample, with a L′ = 7.2 mag. Similarly to G263.7434, a 10.4 μm the image shows a single central object while the 2MASS K-band image indicated a more crowded region (Mottram et al. 2007).

No signature of a companion star was found in the NACO L′−band images (see Appendix A), for the probed separation ranging from ∼440 to 54,500 au. We derive a 5σ detection limit of Δ = 2.2^m.

4.3.3. G305.2017+00.2072A

G305.2017+00.2072A (hereafter G305.2017) is located at a distance of ∼4 kpc and is thus the most distant MYSO of our sample. With such a large physical distance, NACO probes for separations between 640 and 56 000 au. We focus here on the multiplicity of G305.2017 but additional notes are provided in Appendix D as the protostar is a broadly studied object and its geometry is well constrained (Frost et al. 2019; Krishnan et al. 2017).

Liu et al. (2019) review all the detections surrounding G305.2017, as of 2019, for a large wavelength coverage, ranging from sub-mm to NIR. Their high-spatial-resolution observations in the MIR (8.8−24.5 μm) show an emission component in the NE of the massive protostar lying at about 1″ (∼4000 au), that they name G305B2. However, their study shows that this emission is no longer visible in shorter wavelength filters, at 3.8 μm, for example, which is our band of interest. Their MIR observations probing a wider environment also show strong emissions towards the east of the central object, that peak around ∼30 μm and that they call G305A. Our high-angular resolution NACO observations look over similar separation ranges and have the capacity to resolve both objects. Nonetheless, images do not show the presence of a NE-4000 au-orbiting companion and rather confirm the absence of such NIR emissions corresponding to G305B2 or G305A. Pomohaci et al. (2019) do not mention the presence of any additional sources near G305.2017 in their NACO K-band images. We report the detection of two new companions (L′-band), with projected separations of respectively ∼1.5″ (∼6000 au) and ∼8.8″ (∼35 000 au) and flux ratios of around 17%±1%. Both objects have relatively low spurious associations resulting from a low stellar surface density in this direction. None of these two objects is consistent with the detections presented by Liu et al. (2019). The fact that two additional objects are detected at longer wavelengths might indicate that the emitting source is in a much younger and more embedded stage (likely a hot core phase) than the central G305.2017.

4.3.4. G310.0135+00.3892

G310.0135+00.3892 (hereafter G310.0135) has been extensively studied under the name of IRAS 13481−6124. The object is classified as a high-mass YSO of an equivalent ZAMS spectral type O9. The central object has an estimated mass of 20 M_⊙ (Grave & Kumar 2009). It is one of the brightest sources in our sample, with L′ = 2.3^m. We provide some extra important facts about this system in Appendix D, as they do not directly impact our results but give a comprehensive and consistent picture of the system. However, it is worth comparing our multiplicity results to those of Pomohaci et al. (2019) in the near-IR K-band.

Located at 3.25 kpc, L′-bright sources with separations ranging from 520 to 60 000 au should be detectable. Pomohaci et al. (2019) detected a companion located at a distance of 2.56″ which coincides, within errorbars, with our detected L′-band source at ∼2.6″. Using the 2MASS catalogue in the K-band, they derived a chance projection of 16% yielding a likely binding connection between the central star and the companion. We derive a chance projection of only 0.11% in the L′-band given the low background surface number density of L-bright stars along the line of sight. Pomohaci et al. (2019) infer a mass of 7.6 M_⊙ when using foreground extinction. We detect a second object lying at ∼16″ in the NW direction. The presumed companion is most likely a background object (P_spur = 69% and is discarded for the rest of the analysis. Wheelwright et al. (2012) present 20 μm images of G310.0135, showing a resolved and slightly extended object in the NE/SW direction that matches the direction of the outflow modelled in Kraus et al. (2010). The contours of the MIR emission also seem to show emission in the NE direction of the MYSO, approximately 2−3″ away. However, the authors do not comment on this localised emission. Several studies probed the astronomical-unit scales by means of high-angular resolution interferometry: none of which reported the presence of an inner companion. It seems that the companion discovered by Pomohaci et al. (2019) whose detection is confirmed here, is the closest object (∼8400 au) orbiting the MYSO.

4.3.5. G314.3197+00.1125

G314.3197+00.1125 (hereafter G314.3197) is an MYSO located at a distance of ∼3.6 kpc. Both the 20 μm and 10.4 μm images show a single object, without any MIR bright companion in the ∼10″ surroundings (Mottram et al. 2007; Wheelwright et al. 2012). Three NIR-bright sources are seen in the 2MASS K-band images.

We do not detect any extra source in the 28″ × 28″ L′-band NACO images (∼580−50 000 au). We determine companion detection limits on a 5σ level of Δ = 3.4^m.

4.3.6. G318.0489+00.0854B

G318.0489+00.0854B (hereafter G318.0489) is an MYSO associated with two H II regions whose radio counterparts have been established by Urquhart et al. (2007), Navarete et al. (2015). G318.0489 is the faintest object in our sample with L′ = 7.7^m. The 20 μm image from Wheelwright et al. (2012) shows a single resolved object with a slight cometary-like morphology towards the east direction.

Within the probed separation range (540−65 000 au), no companion source was detected around G318.0489 in the L′-band NACO images. We place a 5σ detection limit of Δ = 2.7^m.

5. Discussion

Our sample composed of 13 MYSOs offers an exclusive opportunity to study multiplicity among massive young stars for separations ranging from ∼300−50 000 au. In this section, we aim to compare our multiplicity results to those found for a wide range of masses (∼0.05−25 M_⊙), evolutionary stages (from (M)YSO to PMS/ZAMS and MS) and probed separations (0.1 − 10⁴ au). We further discuss the implications for massive multiple formation.

5.1. Multiplicity of MYSOs

We have identified 8 potential multiples (1 binary, 2 triples, 1 quadruple and 4 higher-order multiple systems) from 13 targets. We rejected 2 detected sources that have a high probability of chance projection (P_spur > 20%: G310.0135B and G268.3957A). This results in a master sample of 39 physical companions, detected around 13 MYSOs. Our newly discovered companions comprise a significant minority of all known systems: only a few of the detected companions have already been reported in the literature. A compilation of the measured physical properties of the multiple systems is provided in Table 2.

Multiplicity and companion fractions (MF, CF) can be computed in different ways but we adopt here the quantities recently presented in the review by Offner et al. (2023). For consistency, we have defined two ranges of physical separations in the following analysis. We derived the physical separations of the multiple systems using the distances (kinematic, parallax and H₂O masers) shown in Table 1. We provide a global MF and CF for the whole separation range probed in this study, 600−35 000 au. For this parameter space, we consider 39 companion detections. We define another range as 600−10 000 au, resulting in a lower sample of 31 companion detections. However, the number of multiples remains unchanged, consequently, the frequency of multiples is identical for both subgroups.

We estimate a raw multiplicity frequency of MF_{600 − 35 000} = MF_{600 − 10 000} = 62 ± 13% and a companion fraction of CF₆₀₀₋ 35 000 = 3.0 ± 0.5 for the complete dataset and CF_{600 − 10 000} = 2.4 ± 0.4 for the restricted group, agreeing with each other within the errors. The uncertainties on the MF follow binomial statistics while the uncertainties on the CF follow Poissons statistics, as described in Sana et al. (2014) and Moe & Di Stefano (2017). Our survey probes the multiplicity of massive stars at their youngest stage and complements well-existing surveys which focus on closer orbits, from spectroscopic to interferometric regimes. Comparing our findings derived from a group of MYSOs with those of different classes of objects and evolutionary stages is not a trivial exercise and has some limitations. Hence, it is important to consider statistics tracing similar ranges of separations, but also the specific observations and techniques used in each study. One should keep in mind that the statistics are highly affected by the homogeneity of the samples, the sample size and the sensitivity limits. As a result, the statistics associated with low-mass stars campaigns are much more robust than for their massive counterparts. In fact, the samples contain several hundreds or even thousands of low-mass objects while the studies on young massive stars contain at most a few tens of objects. It is for these various reasons that we decided to split our results into two different separation ranges, for more consistency and a fairer comparison with the existing studies, especially the one of Pomohaci et al. (2019), for which the comparison is among the most robust ones we can refer to. Figure 3 gathers a handful of multiplicity results for different populations of young stars, in clusters and field, targeting different mass ranges and separations.

Fig. 3. Multiplicity frequency for young stars (low- to high-mass stars) for different spatial scales: from 0.1 to several hundreds of thousands of au. The horizontal bars of different colours indicate the separation intervals specific to each study. The different symbols refer to a category of young stars (see legends).

For similar separation ranges (600−10 000 au), the previous survey of Pomohaci et al. (2019) reports MF_{600 − 10 000} = 31 ± 8%, lower than our results despite a larger sample size and average distance (3.1 kpc). Unlike their conclusions, we notice a significant difference in the MF over distance. Among the seven sources lying at distances shorter than 2 kpc, six are multiples yielding MF_{< 2 kpc} = 85 ± 13% while we measure MF_≥2 kpc = 33 ± 19% for sources further away than 2 kpc. The distance cut we used here seems to show that at further distances, the fainter companions are not detected. A similar comparison can be drawn for the MF as a function of the bolometric luminosity. We count nine MYSOs whose bolometric luminosity is under 10 000 L_⊙ and four brighter than 10 000 L_⊙, yielding MF_{< 10 000 L⊙} = 55 ± 17% and MF_{≥10 000 L⊙} = 75 ± 22% respectively. For both luminosity ranges, the MFs are consistent with each other within the errors. Keeping the comparison strictly to the MF regardless of the distance towards the source, a recent study lead by R. Shenton shows notable consistency with our derived MF, for similar separation probes (see Fig. 3). Using K-band data from VVV, UKIRT and UKIDSS, they report a MF varying from 54 to 70%, for a population of 402 MYSOs and spanning separations of 10³ − 10⁵ au with a mean physical separation of companions of ∼18 000 au (Shenton et al. 2023). The similarity with our work is all the more encouraging given their much higher sample size. The main asset of their analysis is the deeper limiting magnitude at which they can detect companions, however, their spatial resolution stays below the NACO’s spatial resolution presented in this work or in Pomohaci et al. (2019): the inner ∼1 to 1.5″ of each MYSO is essentially a blank spot and only the brightest companions between 1.5 and 2″ are detected, with a peak of overall detected companions between 3 and 6″.

Zooming in at the interferometric scales, the results of Koumpia et al. (2021) whose sample was composed of six MYSOs with masses ranging from 8.6 M_⊙ to 15.4 M_⊙, are somewhat smaller than ours but in agreement with those of Pomohaci et al. (2019). They report a rather low MF of ${17}_{-17}^{+21}$%. These statistics are in stark contrast with those of GRAVITY Collaboration (2018) and Bordier et al. (2022) who report an MF of 100% for relatively small samples of pre-main-sequence massive stars (4 and 6 stars respectively), for similar traced separations. The latter results are consistent with the more evolved populations of massive stars (Sana et al. 2014, see also Fig. 3) but raise a major question about the evolution of multiplicity between the youngest and more evolved stages. Dynamic rearrangements within dense clusters must be at play to shape the final structures of the systems and could point towards an agreement with the competitive accretion models which require the stars to interact via mass or momentum exchanges in order to shape the final structure of the cluster. This point is further discussed in Sect. 5.2.

Moving on to spectroscopic separations, Kounkel et al. (2019) and Zúñiga-Fernández et al. (2021) studied the low-mass regime of young stars, the former focusing on a sample of Class II/III pre-main-sequence T Tauri stars and the latter on a sample of 410 nearby YSOs (< 200 pc). They respectively derive a multiplicity fraction of 30 ± 6% (averaged across their three groups) and 33 ± 3%. We also mention the results of Shenton et al. (2023) who investigate close binarity of 29 YSOs by the means of high-resolution spectroscopy with X-shooter/VLT and find a multiplicity occurrence varying from 35 to 55%. All these values are in alignment with the findings of Tobin et al. (2018, 2022) that derived similar protostellar multiplicity statistics (30 ± 3% and 38 ± 7%), considering separations as small as 20 au and up to 10⁴ au (ALMA observations), within their sample of young low-mass stars in Perseus and Orion star-forming regions, respectively. Meanwhile, the close binary fraction within a sample of Herbig Ae/Be protostars, whose masses lie between low-mass YSOs and the MYSOs we study here, inferred by Apai et al. (2007) reaches ${20}_{-0}^{+10}$%. This is within the errors of the above-mentioned results, but considerably lower than the wide binary fraction derived by Wheelwright et al. (2010) for the same class of protostars (74 ± 6%).

Figure 3 does not show a clear trend in the multiplicity statistics of young stars, independently of their mass and environment. It seems that two trends are emerging: a significant number of studies (Apai et al. 2007; Elliott et al. 2015; Duchêne et al. 2018; Tobin et al. 2018, 2022; Kounkel et al. 2019; Koumpia et al. 2021; Zúñiga-Fernández et al. 2021) have recently reviewed the multiplicities of several samples and whatever the mass and separation explored, the multiplicity seems to vary between 20 and 50% (bottom part of the figure). The other trend appears for multiplicities above 50%. We note that no spectroscopic study unambiguously derives a multiplicity fraction above 50%. Assuming that this class of stars is observable with the current technical means and that they meet the requirements of instruments’ sensitivity, the low spectroscopic multiplicity found among MYSOs is questionable. One could naturally conclude that 1 Myr old massive stars do not have any close companions. The present study shows that companions can be much fainter than the primary, which might also be the case in the spectroscopic regime. The most massive star could be two or three orders of magnitude brighter than the faint companion in which case the primary would overrule the spectroscopic emission and hide the spectroscopic features from the companion, preventing one to proceed with a double spectroscopic binary analysis. Similarly, a low-mass object would not significantly affect the orbit of a higher-mass primary, falling off the radial velocities detection limits. In the scenario in which there is a dearth of spectroscopic companions as it was shown in M17 (Sana et al. 2017; Ramírez-Tannus et al. 2021), the companions lying at intermediate separations (GRAVITY Collaboration 2018; Bordier et al. 2022, as observed by) can eventually experience an inward migration on Myr timescales, as theorised by Sana et al. (2017). On the other hand, most interferometric studies, whether conducted with ALMA or VLTI, converge towards rather high multiplicities in the 1−200 au scales, (Wheelwright et al. 2010; GRAVITY Collaboration 2018; Bordier et al. 2022), with the exception of the study conducted by Koumpia et al. (2021), for which the MF is lower but the evolutionary stage differs slightly as their objects are all MYSOs. These latter results have to be taken with more caution as the probed stars tend to have already reached the birthline. Here again, the physical distances to the object are an important source of error and we recall that unlike young low-mass stars, the majority of massive stars are found at large distances. Thus, the average distance of the Pomohaci et al. (2019) study for instance, is 3.3 kpc while that of our study is slightly more than 2 kpc. In contrast, Tobin et al. (2018)’s studies are settled for regions at 200 pc from the Sun, a factor of 10 closer.

We further investigate our findings for a specific class of objects: massive stars. To do so, we relax the evolutionary stage parameter and compare our results to those of different studies (Apai et al. 2007; Sana et al. 2014; GRAVITY Collaboration 2018; Pomohaci et al. 2019; Koumpia et al. 2021; Rainot et al. 2022; Bordier et al. 2022). The results are shown in Fig. 4. We keep the comparison consistent by providing statistics tracing similar separation ranges. From 600 to 10 000 au we note that the multiplicity frequency is consistent across all the evolutionary stages and within error bars, with the exception of the study conducted by Pomohaci et al. (2019) who found a lower fraction, almost two times lower, for the same probed separations. However, looking at the companion rate panel, the frequency of companions at wide separations seems to drop across the evolution, while it is higher when studies go down to the spectroscopic regime. Indeed, the SMASH survey that assesses the companion fraction among main-sequence O stars shows that the overall companion fraction (resolved+unresolved) is more than a factor of two higher than the statistics for the sole group of resolved companions (Sana et al. 2014). It may indicate that the presence of close companions is preponderant when deriving the companion rates for main-sequence massive stars, while it does not seem to be the case for young massive stars. For instance, the spectroscopic statistics derived by Apai et al. (2007) show a rather low limit for the binary fraction. This trend is consistent with the idea that massive star formation does not result in a tight, close binary system, but rather forms at large distances and experiences dynamical effects over time (Ramírez-Tannus et al. 2021). Overall, we note that our derived statistics are in better agreement with those derived for older ages than the statistics derived for similarly aged populations (Pomohaci et al. 2019; Koumpia et al. 2021). It is clear from Figs. 3 and 4 that the multiplicity statistics are strongly dependent on the mass of the targeted sources.

Fig. 4. Multiplicity statistics, including the MF (right) and CF (left) for massive stars and different evolutionary stages. The different colours and symbols stand for the different probe separation ranges, as stated in the legend box.

We calculated the physical separations of the multiple systems using the RMS kinematic distances, displayed in Table 1. Figure 5 shows the projected physical separation of the detected companions per system. The size of the markers refers to the brightness ratio in relation to the primary. The projected separation is highly dependent on the physical distance towards the primary star. The closest companion detected is G268.3957B, at around 1350 au of its primary. The plot suggests a preferred separation range for companion stars, which likely indicates the separation at which they were formed given the young age of the central object. We clearly observe that the vast majority of the companions (80%) are detected within 10 000 au around their central MYSO and at least one companion per system is observed in this parameter space. For a given system, the nature of the detected companion is rather heterogeneous, faint companions are found at a separation of ∼5000 au as well as brighter companions. Figure 5 shows the powerful nature and effectiveness of AO-imaging for detecting new and faint companions, down to a factor of about 500 times fainter than the massive central YSO in the present study (Δ_max = 6.8^m for G203.3166-I at a distance of ∼6800 au).

Fig. 5. Companion separation per system with given their L′ magnitude and distance in kpc. The colours indicate different systems and the marker size scales with the magnitude difference in relation to the primary. The vertical grey line indicates the 10 000 au separation that we used as a cut earlier for a fairer comparison with other studies. Small variations in the y-direction are to avoid overlapping markers.

In view of the results, there is no clear correlation between the separation and the type of companion that is formed in the vicinity of the massive primary. It is important to remind that we observe the L′ emission in the images, which includes for instance the dust emission, whose radiation comes from the natal envelope and the circumstellar disk. Thus, L-band magnitudes cannot be trivially used as a proxy for stellar mass companions as it is commonly the case for K-band magnitudes (Pomohaci et al. 2019; Rainot et al. 2020; Reggiani et al. 2022; Pauwels et al. 2023). In fact, we lack crucial information such as a mass–luminosity relation for the L′ band, an accurate distance, and both the extinction (interstellar and circumstellar) and embeddedness of the systems. Masses are commonly estimated using comparison with evolutionary tracks (Bordier et al. 2022) or directly from the observed spectral type of the primary, which strongly depends on the extinction correction (Beltrán & de Wit 2016). All these unknowns could be hypothesised, such as assuming that the extinction towards the primary spreads for the whole system, but this is less likely to be valid for distant companions (as it is the case here) since these companions may not be obscured by the same amount of dust as the primary and such an assumption would result in an overcorrection. Combined with the uncertainty about the age of the systems, the derived masses would include large error bars which would not allow us to have a clear view as to the real mass of the formed companions or their spectral type. In light of the present results, the only conclusive statement we can make about the mass of the companions is that they appear less massive than the central massive protostar, assuming that the companions are shrouded by the same or lesser amount of dust than the primary.

That being said, companion separation alone allows us to discuss the most likely scenarios for the formation of such systems, which we highlight in the next section.

5.2. Implications for star formation models

As predicted by the initial mass function (IMF), massive stars (M_i, primary > 8 M_⊙) make up < 1% of the total stellar population, meaning that many MYSOs are most likely the most massive objects locally. Observationally, regions of massive star formation in our Galaxy are located in more remote regions than those of low-mass star formation, which has the effect of considerably reducing the level of detail that we can probe and rendering their study arduous. Unlike low-mass star formation, obtaining a clear view of the formation process of massive stars with a high level of detail is difficult. As such, defining the frontier of similarities and differences between low-mass and high-mass star formation is not a trivial task, with an added complexity when multiplicity is at play. We mainly rely on (3D) hydrodynamic and magneto-hydrodynamic (MHD) simulations of collapsing molecular clouds to understand all the physical processes involved (Krumholz & Bonnell 2007; Myers et al. 2013; Meyer et al. 2018; Oliva & Kuiper 2020; Mignon-Risse et al. 2021, 2023). Because massive binary formation involves many complex processes, various observational campaigns aiming at probing different regions around MYSOs as well as studying the substructures reinforced the numerical works and yield a better understanding of individual effects at one time.

Several theories may explain the origin of massive binaries. A complete review of the current models and modes of fragmentation is provided in Offner et al. (2023). Here we provide the three main and independent scenarios for multiples formation and discuss the most likely channels for our observations. Owing to the probed spatial scales in this work that set a limit on multiple system separations, our results can only place constraints on formation mechanisms at core scales (> 1000 au).

The gravitational fragmentation of massive cloud cores or filaments into two Jean’s unstable over-densities, providing a sufficient reservoir of gas, may explain the formation of wide binaries (> 1000 au). In that specific scenario, the main triggers of fragmentation are rotation and turbulence which create density and velocity asymmetries. The Jeans-radius being of the order of 10 000 au, two wide and bound components can be generated. Some characteristics may be imprinted during formation that are still detectable when catching young multiple systems right after the birth parent cloud has dispersed. The initial separation is one of those predicted signatures: a lower limit is set at ∼100 au, strictly fixed by the effects of magnetic support, angular momentum and tidal forces (Lee & Hennebelle 2018), and 0.1 pc is proposed as an upper limit as two fragments have nearly no chance to be initially bound at such separations and will inevitably break off. Turbulent fragmentation also tends to produce randomly distributed alignments of stellar spins, accretion disks, and protostellar outflows (Offner et al. 2016; Bate 2018; Lee et al. 2019). Lee et al. (2019)’s numerical simulations confirmed the existence of such wide embedded and accreting multiples for a short amount of time, before they migrate towards closer orbits, at a couple hundreds of au. Pomohaci et al. (2019) and the present study also shows observational evidence for such wide companions. At this stage, we cannot confirm that core fragmentation alone explains the observed structures, as we would need to probe the innermost regions of MYSOs to seek for the presence of closer companions and conclude as to the other mechanisms that could simultaneously take place. We point out the case of G282.2988 around which Koumpia et al. (2019) showed the presence of an interferometrically resolved companions at ∼50 au while Pomohaci et al. (2019)’s study had already reported two wide companions at ∼5800 au and ∼10 000 au respectively. The complementary of the probed spatial scales here, shows the existence of a wide range of companions, however, given the very young age of that source, core fragmentation would fail at explaining the origin of the 50-au close companion. This brings us to the next theory that may explain the presence of closer companions.

Independently from the previous mechanism and during the monolithic collapse of a massive core, the required massive circumstellar disks which need to feed the central object to high masses are prone to gravitational instabilities and lead to fragmentation and subsequent formation of companions (Krumholz & Bonnell 2007; Kuiper et al. 2010; Meyer et al. 2018; Oliva & Kuiper 2020; Mignon-Risse et al. 2023). The main physical mechanism responsible for the creation of spiral arms where fragments are created resides in the so-called Toomre instabilities (a gravitationally unstable disk will have Toomre Q parameter < 1, Toomre 1964). In fact, it has been demonstrated that early in evolution, disks reach a sufficient mass to become self-gravitating, form spirals and fragment. An increasing amount of observational evidence and characterisation of disks around MYSOs and its substructures, including the detection of a fragment in the outskirts of the accretion disk, have been reported since the advent of high-angular resolution facilities such as ALMA and the VLTI (see for e.g. Kraus et al. 2010; Cesaroni et al. 2017; Beuther et al. 2017; Ginsburg et al. 2018; Ilee et al. 2018; Frost et al. 2019, 2021; Maud et al. 2019; Johnston et al. 2020). These freshly created fragments can grow in mass through secondary smaller disks that form locally. If they successfully overcome the fragmentation epoch imprinted with many mechanisms of destruction such as shearing, collisions, accretion bursts or being accreted by the central protostar, they can become actual companions. Ultimately, simulations show that the central massive protostar gains sufficient mass through different episodes of accretion and the orbits of the surviving fragments shrink, which provides an ideal framework for the formation of spectroscopic binaries (Kratter & Matzner 2006; Oliva & Kuiper 2020; Tokovinin & Moe 2020). Low-mass cluster member stars would thus be formed all within the size of the accretion disk, giving rise to an unstable small-N system, prone to (partial) dynamical disruption. In order to explain binaries with separations smaller than 10 au, the scenario fragmentation then migration has gained popularity since massive disks inefficiently fragment at distances smaller than 50−100 au from the central object (Krumholz et al. 2009; Kratter et al. 2010; Zhu et al. 2012). In this way, binaries can shrink their orbit by one or two orders of magnitude after they have formed at wider separations (Moe & Kratter 2018).

The observational predictions supporting disk fragmentation is undoubtedly the presence of companions within the expected disk size around MYSOs (up to 2000 au, refer to the previous references for disk sizes around massive stars). Another reliable test would be to measure the alignment of the accretion disk with that of the binary orbit. In fact, disk fragmentation would lead to a final orbit within the same plane as the accretion disk. We detect four companions within 2000 au, in three different systems (see Fig. 5). NACO offers sufficient angular resolution to probe closer regions but these companions are not found. Several possibilities may explain the observed lack of closer companions: the systems may be very young such that these companions have just formed in the outskirt of the disk and have not experienced, if any, inward migration nor any other dynamical interaction. Conversely, disk fragmentation could have been successful at forming companions within the accretion disk, but substantial torque generated by small individual secondary disk in addition to the ongoing accretion could push the companion outwards on a wider orbit, as demonstrated by Muñoz et al. (2019), in the case of binaries embedded in disks. Another explanation is that companions do not form in the 100−2000 au but rather at larger scales or in the innermost regions that are not probed in the present study, but that were detected in (Kraus et al. 2017; Koumpia et al. 2019) for example. Moreover, any companion that would have already cleared up its dusty cocoon might be opaque at the probed wavelength. Though this assumption is less likely as the central source might evolve faster than the surrounding lower-mass companions. Lastly, we cannot rule out the possibility of highly eccentric orbits (Maíz Apellániz et al. 2017; Tokovinin 2020). As a deduction, the companion spends a consequent amount of time at a large separation before reaching the periastron location, which can be much closer to the MYSO. Statistically and with only one snapshot in hand, it is most likely that we are overestimating the separation to the primary. Only high-precision astrometric measurements can provide accurate eccentricity measurements and eliminate this hypothesis. Even though such measurements may also provide key constraints on the binary formation mechanisms, they are however challenging to obtain as binaries with separations of ∼1000 au have orbital periods of several thousand years, which would translate into very tiny variations of the overall orbit by performing a typical monitoring with a baseline of a few years (Hwang et al. 2022).

Dynamical interactions and capture may occur during the dissolution of the forming regions, that is, post fragmentation, where an initially single and unbound star is paired with or attracts another nearby star or system, due to energy and angular momentum exchange. The high degree of clustering in most star-forming regions facilitates the interactions between stars, readjusting the hierarchy, multiplicity and orbital parameters of the systems formed via the two above-mentioned channels. The timescales and upper limit scales at which capture occurs are uncertain. This mechanism would explain the formation of both the tightest and widest systems. Murillo et al. (2016) demonstrate that a significant fraction of high order multiple systems are non co-eval, showing evidence for different episodes of core fragmentation and possible reconfigurations via dynamical interactions. Most of our targets do not seem to belong to particularly crowded regions, as per the derived local stellar densities in the L-band (see Sect. 3, at the exception of G265.1438). It also questions the real binding connection between the central object and the detected sources and their association to the cluster rather than to the system. At this point, it is not an easy task to assess whether all these associations shared a common origin, due to the absence of proper motion measurements.

Orbital hardening through disk-driven migration, causes the shrinkage of the inner orbits and is likely associated with the ejection of outer low-mass companions. Dynamical-friction forces are also at play and could efficiently reduce orbital separations from thousands of au to hundreds or even tens of au, in less than a few Myr (Ramírez-Tannus et al. 2021). Such inward migration is consistent with disk fragmentation simulations and produces spectroscopic binaries (Oliva & Kuiper 2020). Kozai-Lidov cycles coupled with tidal friction are yet the most efficient way to produce the tightest binaries, caused by the presence of an inclined outer tertiary (Kozai 1962; Lidov 1962; Michaely & Perets 2014; Toonen et al. 2020). All these dynamical interactions are the ultimate opportunity to shape the final multiplicity distribution observed among pre-MS and more evolved field star populations, however, they highly depend on the birth local environment, such as radiation and protostellar outflows.

On that latter point, we want to point out the case of G203.3166. Parker & Schoettler (2022) recently published a study aiming at tracing back the formation history of NGC 2264, star-forming region that hosts G203.3166. They isolate two main sub-clusters of about 2 pc radii, based on the derived age dichotomy: S Mon (∼5 Myr) and IRS 1/2 (∼2 Myr) sub-clusters. G203.3166 (IRS1 in their study) is thus the main driving force of one of the sub-clusters. They report that NGC 2264 might have experienced multiple episodes of star formation, given the different age of the sub-cluster without further iterating on the preferred pathways for their formation. Flaccomio et al. (2023) add that NGC 2264 is not dynamically relaxed and its present configuration is the result of multiple dynamical processes that are most likely still ongoing. If we focus on G203.3166 itself and our results, providing that the observed objects are physically associated with each other, if one looks closely at the disposition of companions around G203.3166, one should notice the distribution of the objects in an arc on the SW side of the MYSO. This specific arrangement would result in a high-pressure front or shock wave formed by a trigger that is most likely G203.3166 being the brightest and most luminous object in the field (Allen 1972; Thompson et al. 1998). The reported presence of an outflow further supports that hypothesis, in which G203.3166 would form first and the remaining stars are formed within the next 2000 yr, of the onset of the outflow assuming a typical outflow velocity of 10 km s⁻¹ and a distance of 4000 au. After that estimated time, most of the gas and dust reservoir needed to form new stars would have been swept away. This timescale argues strongly for the concept of triggered star formation.

6. Conclusions

This study was designed to witness one of the earliest stages in the lives of massive stars, after the formation of the star itself in order to tackle the uncertainty as to the formation process by characterising the epoch of the onset of stellar multiplicity. In this work we present the multiplicity statistics in the L′-band for a sample composed of thirteen MYSOs, from AO-assisted imaging VLT/NACO observations. The parameters from individual observations including night conditions are gathered in Table 1 and a summary of the detected companions per system is provided in Table 2. We have derived multiplicity properties of our sample within the limited angular range of 0.16″ − 28″, corresponding to projected separations of ∼120 au for the closest targets (∼0.7 kpc) to ∼60 000 au for the widest ones. These results are compared to other populations by conscientiously selecting the comparison criteria (evolutionary stage, mass, and separation range) and seeking statistical similarities or differences. Our main results can be summarised as follows:

1. We report the identification of eight multiple systems (one binary, two triples, one quadruple, and four higher-order systems), with a total of 39 companions detected from the observation of thirteen MYSOs.

2. The detected sources cover a wide contrast range of between Δ = −0.1^m and Δ = 6.9^m but the majority of them are found in the Δ = 2.5^m − 5^m range. The poor local source density in the L band in most of our fields leads to a negligible chance of spurious association apart from two candidates that we discarded in the statistics.

3. No clear-cut trend could be found in the radial location of the sources, or in their brightness ratios. However, we note that nearly 70% of the companions are to be found between 3000 − 10 000 au and that at least one companion per system is to be found within 10 000 au.

4. When no companion is found within the field, we provide a detection limit computed from the artificial injection of fake companions, demonstrating that NACO reaches the necessary capabilities to detect faint (L′∼9) companions over the entire FoV of our observations.

5. The detection of companions at different scales informs us on possible pathways for massive star formation and evolution, regardless of their separation. The presence of wide companions (≳1000 au) indicates that they rather formed via fragmentation during the initial free-fall collapse, in the earliest stages of massive star formation. For closer companions and given the large scales that the accretion disk can reach around MYSOs, fragmentation can occur later on in the outskirts of the disk and produce companions. This specific category of companion is most likely to be affected by dynamical interactions such as inward migration and meet the orbital properties measured in older populations of massive stars.

Given the data currently available, including our results and the broader literature, it is still too early to properly conclude whether the observed structures preferentially follow one formation scenario over the other. Further attempts to observe these objects at higher angular resolutions with more sensitive instruments such as the VLTI/GRAVITY will allow us to cover the missing innermost separation range and as such detect the substructures involved in the formation of the central star (accretion disk and outflows) or closer companions. In addition, observing these sources in a band with a well-known mass-luminosity relation would help in assessing the mass and spectral type of the detected companions, thus discarding some dynamical interactions that are less likely to occur for a given mass ratio. Probing the full range of separations and having proper mass estimates would assuredly help us better understand the interactions between the disk, the host star, and its companions.

Acknowledgments

E.B. gratefully acknowledges support from the studentship programme founded by the European Southern Observatory. This project has further received funding from the European Research Council under European Union’s Horizon 2020 research programme (MULTIPLES, No. 772225). This publication is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 090.C-0207(A). This paper made use of information provided on the Red MSX Source archive at http://rms.leeds.ac.uk/. This research has made use of the SIMBAD database and VizieR catalogue access tool, operated at CDS, Strasbourg, France. Part of our code made use of Astropy (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018) and of the python code PHOTUTILS. We used the internet-based NASA Astrophysics Data System for bibliographic purposes.

References

Appendix A: Overview of the observed multiples

We provide here an overview of the observed multiple systems. The left hand panels correspond to the full FoV of 28″ × 28″ and the right hand panels display a zoom-in (15 × 15 pixels) on the individual detected sources.

Fig. A.1. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.

Fig. A.2. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.

Fig. A.3. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.

Appendix B: Summary of single objects

We display the 28″ × 28″ reduced NACO L′-band images of sources for which we do not detect companions.

Fig. B.1. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.

Appendix C: Detection limits

Here we show the limiting magnitude determination by injecting artificial sources in all the images. Purple points correspond to artificially injected sources, not detected by our source detection algorithm at a 5σ level, while the light blue ones show the non-detections. For reference, we added to the plot all the detected sources (orange symbol of a star) in the real image and provided for both groups the lowest detected magnitude. The y-axis values are the differences in magnitudes between the artificial sources and the primary source, and the x-axis is the distance in arcseconds between these two. Using that method we ensure that all the sources within the sensitivity of the instrument and matching the requirements of our source detection algorithm are detected. Any missed source would either be fainter than the limiting magnitude, masked by the background variations or untrustworthy detections.

Fig. C.1. Detection limit artificial source tests for all the targets. We display only one iteration. The orange star symbols stand for the real detected companions (if any) and the plain dot symbols refer to artificially injected sources. We display the magnitude contrast of the faintest (i) source detected in the real image (orange dot and line) and (ii) artificial source (purple dot and line).

Fig. C.2. Detection limit artificial source tests for all the targets. We display only one iteration. The orange star symbols stand for the real detected companions (if any) and the plain dot symbols refer to artificially injected sources. We display the magnitude contrast of the faintest (i) source detected in the real image (orange dot and line) and (ii) artificial source (purple dot and line).

Appendix D: Notes on the previously studied MYSOs

G194.9349: This source has been the object of NIR spectroscopy Cooper et al. (2013) and far-IR observations to detect the presence of outflows and characterise the infall dynamics through the SiO(8-7) emission line on one side, and the HCO⁺/H¹³CO⁺ emission on the other hand (Cunningham et al. 2018). They report a non-detection of SiO but the presence of both HCO⁺ and H¹³CO⁺. The former non-detection could either be associated with the absence of an active outflow or a potential fossil remnant-driven outflow. The structures associated with G194.9349 cannot be confirmed at this stage and may require higher spatial resolution observations where infall motions onto the protostar and outflows can be resolved.

G203.3166: G203.3166 known under various names (NGC224 IRS1, AFGL989, Allen’s source) has been studied at different wavelengths and angular resolutions. Schreyer et al. (2003) described G203.3166 as a young B-type star, devoid of disk but surrounded by low-mass companions located in a low-density cavity and surrounded by denser gas. As part of the commissioning of the NICMOS (K−band) camera on board Hubble Space Telescope (HST), Thompson et al. (1998) report the discovery of six point objects at projected separations of 2.6″ to 4.9″ from the central source. In the L−band NACO images it seems that we observe four out of the six objects reported in Thompson et al. (1998): objects G(#6), E(#5), D(#4) and C(#2). Schreyer et al. (2003) further complete the picture with K− band (SOFI at ESO’s NTT) and mid-infrared (4.6 and 11.9 μm, TIMMI2 st ESO’s 3.6m telescope) observations. Within the field of view of our interest, they confirm the presence of the 6 objects detected by Thompson et al. (1998) and report the detection of two new objects (#7 and #8), that respectively coincide with object F and the binary system (A and B) in our images. The K−band polarisation map displayed in Figure 2, clearly shows that object #8 has its own polarisation pattern confirming the presence of circumstellar dust. Object #8 being brighter at longer wavelengths further supports the presence of dusty circumstellar material. The advent of AO-assisted images allowed us to resolve the binary system as shown in our images, which was not the case in the NICMOS/Hubble images. de Wit et al. (2009) looked at the 24.5 μm emission of G203.3166 featuring COMICS observations. In their 6″ × 6″ image, they do not detect companions.

Grellmann et al. (2011) explored the closer IR environment of G203.3166 seeking for closer companions. The interferometric observations (MIDI/VLTI) do not provide any hint of resolved companions around G203.3166, in the 30-230 au separation range. However, a model of a geometrically flat but optically thick circumstellar disk successfully reproduces the short-wavelength region of an SED and visibilities (SED taken with SWS/ISO satellite). Finally, a multi-scale analysis by Frost et al. (2021) fitting these visibilities, the 24.5 μm image profile and the SED simultaneously found that the source is likely to harbour a small dust disk (0.3M_⊙ with a maximum radius of 500R_⊙) which contains a gap-like substructure, surrounded by a dusty envelope with outflow cavities.

G232.6207: NIR and MIR observations (K−, L− and N-band) conducted by Walsh et al. (1999) and Walsh et al. (2001) have shown that the source is composed of a compact centre with a very reddened object and diffuse emission. De Buizer et al. (2017) further confirm the presence of a NIR/MIR bright compact core and extended emission, using 11.7 μm Gemini South images, Spitzer 8 μm and 31 and 37 μm SOFIA data. Both Walsh et al. (2001) and De Buizer et al. (2017) highlight the influence of a putative outflow that is driven by a radio continuum source, coinciding with the NIR/MIR source. More recently, Zhang et al. (2019) reported that high angular resolution observations with ALMA of the 1.3mm continuum and the H30α recombination line emission confirmed the central object to be a massive protobinary with apparent separation of 180 au. The massive central binary is expected to have a total mass of 18M_⊙ and a maximum period of 570 years. They also report features exhibiting stream-like structures, at the 10⁴ au scales. Their measurements show evidence for ongoing accretion at different spatial scales: from infalling streams at 1000 − 10000 au, to forming a massive binary on 100 − 1000 au scales via a circumbinary disk and finally to circumstellar accretion disks feeding each star separately on 10 au scales. Apart from the central object, they detect three other compact sources within 17000 au, two of them being low-mass protostars (0.022 and 0.014M_⊙) and lying respectively at 6″ and 9″ away from the central massive binary. The third one 9″ away (0.25M_⊙) from the binary was not retained as a protostellar because of its non-point-like core shape.

G254.0548: G263.7434 has no radio counterpart as found by Urquhart et al. (2007). They report the protostar to be unresolved with spherical morphology. It would seem that this object has not been further studied, there is no archive with images at NIR or MIR wavelengths.

G263.7434: G263.7434 has been observed at 20 μm as part as the VISIR and COMICS campaign lead by Wheelwright et al. (2012). However, they report the protostar to be unresolved and could not conduct the modelling of a putative envelope as they did for the other sources in their sample.

G263.7759: In the resolved 20 μm images, Wheelwright et al. (2012) noticed a bipolar morphology that contrasts with the circular symmetric morphology revealed by most of the resolved MYSOs in their sample. They suggest that the morphology traces the cavities of a bipolar outflow seen close to edge on. Indeed, the extension from the NW to the SE direction they found is in clear agreement with a well-collimated H₂ jet detected by (Giannini et al. 2005). The jet seems to be driven by several different sources that they detected thanks to ISAAC/VLT. The structure of the jets and the various associated knots is also well seen in the OSIRIS/SOAR narrow-band H₂ images Navarete et al. (2015). Being a resolved object, Wheelwright et al. (2012) could model the SED and intensity profile of G263.7759 with an infalling envelope of mass 15.5M_⊙ and a mass infall rate of 1.5 × 10⁻⁴M_⊙yr⁻¹.

G265.1438: G265.1438 is classified as an accreting Class 0/I YSO of about 1.1 Myr (Ellerbroek et al. 2013). The latter authors lead a multi-wavelength study of the RCW 36 region which includes NIR/MIR imaging and photometry (NTT/SOFI and Spitzer/IRAC), as well as NIR integral field spectroscopy (VLT/SINFONI) and optical spectroscopy (VLT/X-SHOOTER). They report that G265.1438 shows a highly reddened spectrum with many emission lines and outflow features. They also detect the presence of a dust shell or disk that might responsible for the severe reddening. Wheelwright et al. (2012) had already reported the possible presence of outflow and their model is compatible with an envelope of mass 15.5M_⊙ and similar mass infall rates as G263.7759, viewed at relatively low inclination (i = 32°).

G268.3957: G268.3957 was observed as part of a VISIR/VLT campaign and was spatially resolved in the MIR. Wheelwright et al. (2012) explore the nature of its circumstellar structure using radiative transfer models. G268.3957 exhibits a rather circular symmetric morphology in the MIR, a consequence of the low inclination at which it is observed (i = 30°). The spatial extension of its 20 μm emission and its SED is better reproduced with a model featuring an infalling (3.5 × 10⁻⁴M_⊙yr⁻¹) rotating envelope (9.4 × 10³M_⊙) with outflow cavities (10° opening angle). No companions were reported within the 30″ × 30″ field of view of the VISIR (20 μm) observations.

G269.1586: Wheelwright et al. (2012) observed G269.1586 as part of their VISIR/VLT campaign and report the object to be spatially resolved in the MIR. They also report a distinct localised emission at the NE in the 20 μm image, that they associate with the H II region detected by Urquhart et al. (2007), in the radio regime. They describe a cometary-like morphology for both objects. Wheelwright et al. (2012) suggest that the exotic morphology, for which they cannot model the intensity profile, may indicate that the object is experiencing a different phase of massive star formation than the other MYSOs observed in their sample, a fact that is supported by the nearby presence of the H II region. As for G268.3957, Wheelwright et al. (2012) could model the SED with an infalling envelope model. This is the most massive envelope they derive, with a mass of 32.2M_⊙ and a mass infall rate of 7.5 × 10⁻⁴M_⊙yr⁻¹.

G305.2017: Using two sets of spatially resolved data in the MIR (MIDI/VLTI and VISIR/VLT) together with an SED, Frost et al. (2019) constrained the geometry of the MYSO. Using a 2.5D radiative transfer models fit to N-band interferometric data, Q-band imaging data and an SED simultaneously, they show that the best model implies the source is surrounded by a large extended dusty envelope featuring bipolar cavities and a dusty disk. They further concluded that an inner hole is likely present in the disk and the low envelope density indicates that the protostar is at an advanced stage in its formation, similar to a transition disk phase of low-mass star formation.

G310.0135: G310.0135 is one of the first MYSOs around which a disk of scales small enough to be an accretion disk was found which has been studied using diverse multi-wavelengths techniques (Kraus et al. 2010; Ilee et al. 2013; Boley et al. 2016; Caratti o Garatti et al. 2015, 2016). The ∼17.3 au-diameter compact disk revealed by NIR K−band interferometry (Kraus et al. 2010) is in Keplerian rotation (Ilee et al. 2013). Perpendicular to the disk plane, (Kraus et al. 2010) observed a molecular outflow and two bow shocks, showing that a bipolar outflow originates from the inner regions of the MYSO. The MIR observations lead by (Wheelwright et al. 2012) that trace the warm dust show an extended emission along the direction of the outflow. Being resolved at ∼20 μm the latter authors could model the SED and intensity profile. The model predicts an infalling envelope of 9.8M_⊙ with a mass infall rate of 7.5 × 10⁻⁴M_⊙yr⁻¹, and an inclination of i = 32°, in agreement with the intermediate inclination of 45 found by Kraus et al. (2010).

G314.3197: Wheelwright et al. (2012) observed G314.3197 as part of their VISIR/COMICS campaign. The 20 μm image shows a spherically symmetric object whose circumstellar environment is unresolved at the resolution of 0.6″, most likely as a result of the great distance to the protostar.

G318.0489: G318.0489 was observed during the VISIR/COMICS observing campaign lead by Wheelwright et al. (2012). They report a slight cometary morphology emanating from the 20 μm image. Facing the inaccuracy around the distance to this source, Wheelwright et al. (2012) chose to not attempt to reproduce its intensity profile and SED as they cannot be unambiguously modelled. Navarete et al. (2015) shows the continuum-subtracted H₂ maps of G318.0489. The area surrounding the MYSO does not show any particular excess in the H₂ emission nor any excess in the continuum filter. However, the protostar is associated with a radio-quiet source.

All Tables

All Figures

	Fig. 1. Zoom over the northern companion system for target G268.3957, showing two distinct components. The white solid lines is the 4-pix aperture applied in the aperture photometry method and the markers show the position of the centroids as found with DAOFIND. For reference, north is up and east is left.
In the text

Fig. 2. Detection limits for each observed field. The star symbol denotes the Δmag of the faintest source detected in the NACO image (if any) while the dot symbol refers to the Δmag of the faintest detected artificial source, injected into the image. The minimum angular separation at which companions are found is about 2.7″ which translates to ∼5000 au at a typical distance of 2 kpc. However, the faintest companion of each field is found at a typical separation of 8.7″ (∼17 000 au).

In the text

	Fig. 3. Multiplicity frequency for young stars (low- to high-mass stars) for different spatial scales: from 0.1 to several hundreds of thousands of au. The horizontal bars of different colours indicate the separation intervals specific to each study. The different symbols refer to a category of young stars (see legends).
In the text

	Fig. 4. Multiplicity statistics, including the MF (right) and CF (left) for massive stars and different evolutionary stages. The different colours and symbols stand for the different probe separation ranges, as stated in the legend box.
In the text

Fig. 5. Companion separation per system with given their L′ magnitude and distance in kpc. The colours indicate different systems and the marker size scales with the magnitude difference in relation to the primary. The vertical grey line indicates the 10 000 au separation that we used as a cut earlier for a fairer comparison with other studies. Small variations in the y-direction are to avoid overlapping markers.

In the text

	Fig. A.1. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.
In the text

	Fig. A.2. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.
In the text

	Fig. A.3. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.
In the text

	Fig. B.1. 28″ × 28″ reduced NACO L′-band images. The Right Ascension (RA) is on the x-axis and the Declination is displayed on the y-axis. For reference, north is up and east is left.
In the text

	Fig. C.1. Detection limit artificial source tests for all the targets. We display only one iteration. The orange star symbols stand for the real detected companions (if any) and the plain dot symbols refer to artificially injected sources. We display the magnitude contrast of the faintest (i) source detected in the real image (orange dot and line) and (ii) artificial source (purple dot and line).
In the text

	Fig. C.2. Detection limit artificial source tests for all the targets. We display only one iteration. The orange star symbols stand for the real detected companions (if any) and the plain dot symbols refer to artificially injected sources. We display the magnitude contrast of the faintest (i) source detected in the real image (orange dot and line) and (ii) artificial source (purple dot and line).
In the text