We present the first search for binary companions of Massive Young Stellar Objects (MYSOs) using AO-assisted $K$ band observations, with NaCo at the VLT. We have surveyed 32 MYSOs from the RMS catalogue, probing the widest companions, with a physical separation range of 400 - 46,000 au, within the predictions of models and observations for multiplicity of MYSOs. Statistical methods are employed to discern whether these companions are physical rather than visual binaries. We find 18 physical companions around 10 target objects, amounting to a multiplicity fraction of 31$\pm$8\% and a companion fraction of 53$\pm$9\%. For similar separation and mass ratio ranges, MYSOs seem to have more companions than T Tauri or O stars, respectively. This suggests that multiplicity increases with mass and decreases with evolutionary stage. We compute very rough estimates for the mass ratios from the $K$ band magnitudes, and these appear to be generally larger than 0.5. This is inconsistent with randomly sampling the IMF, as predicted by the binary capture formation theory. Finally, we find that MYSOs with binaries do not show any different characteristics to the average MYSO in terms of luminosity, distance, outflow or disc presence.

To date, there is no explanation as to why disc-tracing CO first overtone (or `bandhead') emission is not a ubiquitous feature in low- to medium-resolution spectra of massive young stellar objects, but instead is only detected toward approximately 25 per cent of their spectra. In this paper, we investigate the hypothesis that only certain mass accretion rates result in detectable bandhead emission in the near infrared spectra of MYSOs. Using an analytic disc model combined with an LTE model of the CO emission, we find that high accretion rates ($\gtrsim 10^{-4}\,{\rm M}_{\odot}{\mathrm{yr}}^{-1}$) result in large dust sublimation radii, a larger contribution to the $K$-band continuum from hot dust at the dust sublimation radius, and therefore correspondingly lower CO emission with respect to the continuum. On the other hand, low accretion rates ($\lesssim10^{-6}\,{\rm M}_{\odot}{\mathrm{yr}}^{-1}$) result in smaller dust sublimation radii, a correspondingly smaller emitting area of CO, and thus also lower CO emission with respect to the continuum. In general, moderate accretion rates produce the most prominent, and therefore detectable, CO first overtone emission. We compare our findings to a recent near-infrared spectroscopic survey of MYSOs, finding results consistent with our hypothesis. We conclude that the detection rate of CO bandhead emission in the spectra of MYSOs could be the result of MYSOs exhibiting a range of mass accretion rates, perhaps due to the variable accretion suggested by recent multi-epoch observations of these objects.

We present medium-resolution (R~7000) near-infrared echelle spectroscopic data for 36 MYSOs drawn from the Red MSX Source (RMS) survey. This is the largest sample observed at this resolution at these wavelengths of MYSOs to date. The spectra are characterized mostly by emission from hydrogen recombination lines and accretion diagnostic lines. One MYSO shows photospheric HI absorption, a comparison with spectral standards indicates that the star is an A type star with a low surface gravity, implying that the MYSOs are probably swollen, as also suggested by evolutionary calculations. An investigation of the Br$\gamma$ line profiles shows that most are in pure emission, while 13$\pm$5 % display P Cygni profiles, indicative of outflow, while less than 8$\pm$4 % have inverse P Cygni profiles, indicative of infall. These values are comparable with investigations into the optically bright Herbig Be stars, but not with those of Herbig Ae and T Tauri stars, consistent with the notion that the more massive stars undergo accretion in a different fashion than lower mass objects which are undergoing magnetospheric accretion. Accretion luminosities and rates as derived from the Br$\gamma$ line luminosities agree with results for lower mass sources, providing tentative evidence for massive star formation theories based on scaling of low-mass scenarios. We present Br$\gamma$/Br12 line profile ratios exploiting the fact that optical depth effects can be traced as a function of Doppler shift across the lines. These show that the winds of MYSOs in this sample are nearly equally split between constant, accelerating, and decelerating velocity structures. There are no trends between the types of features we see and bolometric luminosities or near-infrared colours.