© ESO 2021

1 Introduction

Massive stars (those with masses greater than or equal to 8 M_⊙) shape the Universe at all scales, with their winds, outflows, and supernovae (SNe) affecting phenomena such as interstellar chemistry, the nature of molecular clouds, and the morphology of galactic superwinds (Leitherer 1994). However, despite the importance of massive stars, their formation is still not completely understood. While forming low-mass stars lose their natal material (envelope) faster than their disk, revealing a clear star-disk system, young massive stars remain deeply embedded. This, combined with the fact that they are numerically rare, makes them more difficult to investigate. One crucial stage of massive star formation is the massive young stellar object (MYSO). Here, the forming star is deeply embedded in its envelope, rendering it optically invisible. Energy does escape an MYSO at micron wavelengths and longer with the most luminosity emitted at around 100 μm, meaning they are IR bright. Despite the high bolometric luminosity of MYSOs, there is no ionising radiation. This implies that the star is still accreting at a high rate, making the MYSO a key stage of the massive star formation process. With current instrumentation available for the IR and (sub-)millimetre wavelengths, similar spatial resolutions (~0.01″) can be reached, allowing direct comparison of morphologies of distinct but related emitting regions and their physical processes (Beltrán & de Wit 2016).

Disks, a key component of the low-mass star formation process, have recently been detected around MYSOs (e.g. Kraus et al. 2010; Johnston et al. 2015; Ilee et al. 2018; Maud et al. 2018). Caratti o Garatti et al. (2017) detected evidence of accretion around the MYSO S255 NIRS 3, observing a radio burst after a maser and IR burst and interpreting this as an accretion event of the order 10⁻³ M_⊙ yr⁻¹. This implies that the detected disks and outflows observed around MYSOs are engaged in the accretion process, as they are in low-mass star formation. The low-mass star formation process has been assigned four distinct evolutionary stages based on the evolution of the envelope and disk (Shu & Adams 1987). However, this has not been achieved for massive star formation as massive protostellar disk evolution in particular is yet to be constrained. Now that they have been confirmed to exist, determining the nature of massive protostellar disks and their role in the accretion process for massive protostars can be investigated.

Various studies have aimed to resolve the near- and mid-IR emitting regions of MYSOs by means of IR interferometry, in particular using the Very Large Telescope and Interferometer (VLTI). For example, de Wit et al. (2010) determined the importance of the bipolar outflow cavities as a source of mid-IR emission through the study of the MYSO W33A. Alternatively, Wheelwright et al. (2012) and de Wit et al. (2009) used Q-band imaging to investigate the envelopes, allowing them to probe cooler regions of the protostellar environments. These studies made use of radiative transfer modelling in order to compare synthetic images to observed ones. These studies also considered spectral energy distributions (SEDs), which provide a multi-wavelength view of MYSOs but only provide, at best, indirect spatial information. As such, they must be combined with spatially resolved observations to provide geometrical information on the MYSOs.

Another potential source of IR photons are the warm circumstellar disks and the dust at sublimation temperature. In Frost et al. (2019) (hereafter Paper I), we combined the above techniques to provide a multi-scale view of the MYSO G305.20+0.21 (hereafter G305). By simultaneously fitting interferometric data between ~7 and 13 μm, Q-band direct imaging data at ~20 μm, and an SEDwith a 2.5D radiative transfer model, the characteristics of the MYSO were obtained. We found that the broad ~10 μ emission traces not only the cavity emission (de Wit et al. 2010) but is also sensitive to the emission originating in the disk. This is in agreement with previous work by Boley et al. (2013), who found that their data from the MID-infrared Interferometer (MIDI, VLTI) can trace both disks and outflows. Furthermore, the inner dust radius, corresponding the inner dust disk boundary, had to be over three times the sublimation radius to satisfy the observables. The SED is dominated by emission coming from the outflow cavities and envelope, constituting reprocessed emission from the star and hot inner disk regions. The fact that the disk is likely to have an inner hole and has rather low envelope density suggests that this source is more evolved in its formation and is more comparable to a transition disk phase, similar to those seen around low-mass stars.

In this paper, we apply the methodology of Paper I to an additional seven MYSOs. In total, this constitutes the largest MYSO sample to date to be subject to a multi-wavelength detailed analysis of this sort. In doing so, we obtain a detailed set of characteristics for each MYSO that can be consistently compared to each of the others, which has been lacking in the literature. In Sect. 2, we describe the observations. In Sect. 3, the methodology is briefly reiterated. In Sect. 4, the main results of the fitting process are described (an in-depth presentation of each source is provided in Appendix A). In Sect. 5, we discuss the results of the sample. Finally, we summarise and conclude in Sect. 6.

2 Observations

The sample of MYSOs consists of eight sources located at a variety of distances, listed in Table 1. The sources were chosen based on the availability of extensive observational data coverage. The types of dataset used between this paper and Paper I are the same, although the exact instrument and source of the data varies slightly, the details of which we describe here.

2.1 N-band interferometric data

An important source for the mid-IR interferometric data of this sample of MYSOs is Boley et al. (2013). These data were all collected with the MIDI instrument. MIDI is a two-telescope beam combiner that operates at N-band and deliversspectrally dispersed interferometric observables. It was fed by the VLT interferometer until its decommissioning in 2015. For all except two sources, the MIDI data used were reduced and previously presented in Boley et al. (2013). In addition to the data from Boley et al. (2013), one configuration of MIDI data from Grellmann et al. (2011) was used for one source (NGC 2264 IRS1) and the MIDI data from de Wit et al. (2010) was used for W33A. The data of W33A were observed using a similar sequence to the Boley et al. (2013) data, using the HIGHSENS mode (where correlated and total flux measurements are taken separately and the photometric observations immediately follow the interferometric measurements to produce the final visibilities) and the same prism. They used calibrators from Cohen et al. (1999) as well as a different version (1.6) of the data reduction software MIA+EWS, but the observing procedure followed was very similar to that of Boley et al. (2013). For the Grellmann et al. (2011) data an estimate of the sky background was made and subtracted over an area corresponding to the size of the mask over which the correlated spectra were extracted. The data from Grellmann et al. (2011) were taken in the SCIPHOT mode using a grism with a spectral resolution of λ/Δλ ≈ 230.

Overall the MIDI data were obtained using a homogeneous set up, with the main difference being the specific u-v coverage of each set of observations. This depends on the specific observing configuration, and the areas probed for each source are illustrated in Fig. 1. A full list of the MIDI observations used in our work can be found in Appendix B.

2.2 Q-band imaging data

This paper uses single-dish Q-band data (~20μn) obtained with the VISIR instrument at the VLT (Lagage et al. 2004) and the COMICS instrument at the Subaru telescope. VISIR is a mid-IR imager and spectrograph (mounted on Unit Telescope 3 (UT3) of the VLT) and was used in service mode to observe the sources G305 and IRAS 17216-3801 as a part of ESO run 097.C-0320(A) (PI: A. J. Frost). The observations were taken with the Q-band filter which has a central wavelength of 19.5 μm and a half-band-width of 0.4 μm, tracing cooler material than that of MIDI and therefore providing an alternative view of the protostellar environments. The data reduction and further details on this VISIR data were referenced in Paper I. COMICS observations were used at a central wavelength of 24 μm for the remaining sources, somewhat longer than the VISIR observations but still sensitive to cavity wall and envelope emission. Images were taken with a pixel size of 0.13 × 0.13 arcsec² and a field of view of approximately 40 × 30 arcsec². The COMICS data presented here were reduced as part of de Wit et al. (2009) and processed with IRAF (Tody 1986).

All Q-band observations were accompanied by PSF standards and all target objects are spatially resolved. HD 123139 was observed as PSF standard on the night of the VISIR observations with a measured FWHM of 0.48″ and α Tau was used as the PSF during the post-processing of the COMICS observations and has a FWHM of 0.34″. We performed aperture photometry on the Q-band images in a standard way comparing the flux behaviour of flux standard throughout the night. In our modelling, we reduce the 2D images to azimuthal averaged radial profiles, which are then compared to synthetic data.

2.3 Spectral energy distributions

The fluxes for the SEDs were compiled on an object-by-object basis. The Red MSX Source (RMS) Survey (Lumsden et al. 2013) database was used in addition to the literature to create an observational SED for each object. The literature was consulted to ascertain whether, as was the case for G305, other sources could be contaminating various fluxes, and if this was the case such measurements were not considered in the fitting. All the included SED fluxes are listed in Appendix C.

3 Methodology

The methods of this work closely follow those of Paper I. In short, the multi-wavelength dataset described in Sect. 2 is used toprobe the protostellar environments of MYSOs at 10 mas and 100 mas scales and as a whole. By fitting one RT model to all these observations simultaneously, the characteristics of the MYSO environment at these multiple scales are probed, providing a multi-scale view of the MYSOs.

HOCHUNK-3D (Whitney et al. 2013) is the code used to generate the RT models from which the required simulated observables are extracted. The model consists of two main geometrical components; an envelope with bipolar outflow cavities and a disk, which surround a central object. Polynomial-shaped or streamlined cavities present the two main geometrical options for the bipolar outflow cavities, and two envelope density distributions – power-law and Ulrich (Ulrich 1976) (a rotating and infalling envelope) – are explored. The disk is described by two separate density components. One is a large-grain disk component and the other a small grains disk component, which are assigned different dust types. The base stellar parameters sampled are the stellar radius, mass, and temperature. We find that it is the stellar luminosity, as opposed to the individual radius and temperature, which affect our observables, so this is discussed in later sections. The disk is modelled via two separate disk components, a large-grains disk (disk 1 – described by Model 1 from Wood et al. 2002) and a small grains disk (disk 2 – described by the interstellar medium (ISM) grain model of Kim et al. 1994). For each component the sample parameters were the inner radius, outer radius, and scale-height. In Paper I, only the large-grain disk component was included; in this work, both disk components are used. This allows for different geometrical arrangements of grains in the disk to be experimented with. The disk scale-height and the inner and outer radii are free parameters for each of the disk components. Additionally, we fit for the total disk mass and the relative contribution by each component. The scale-height exponent and the radial density exponent are fixed at values 2.2 and 1.2 respectively. These values are recommended by Whitney et al. (2013) based on the hydrostatic calculations of D’Alessio et al. (1998). The overall disk parameters were the total mass of the disk and the ratio of this disk mass between the two disk components. For one source, the spiral features of the code were also required, where the free parameters were the pitch angle of the spiral, the summation parameter (that determines the width of the arms), the fraction of mass entrained in the arms and the radius where the spiral arms begin. The cavity parameters investigated were the cavity opening angles, the cavity density exponent, and the coefficient for cavity density distribution. The envelope free parameters were the centrifugal radius, the minimum and maximum envelope radii, and the envelope infall rate.

The dust types within the models are not varied throughout the fitting process. The dust type for the large-grains disk follows a dust distribution based on models from Wood et al. (2002), with grain sizes ranging between 5 μm and 1mm. The cavity and ‘small grains’ disk share the same dust type, based on that of Kim et al. (1994) for the diffuse ISM (with maximum grain sizes of 0.01 μm). For the envelope grains, Whitney et al. (2003) produced a size distribution of grains in agreement with R_v ~ 4, which is typical of dense regions of molecular clouds (Whittet et al. 2001). To include the effects of ices in the cooler outer regions of the envelope a layer of water ice was included as the outer 5% of the grains’ surface. This final percentage was decided upon after variation testing and comparison to polarisation observations, as described in Whitney et al. (2003). Polyaromatic-hydrocarbons (PAHs) and very small grains (VSGs) can also be accounted for in the code, resulting in emission that is visible in the SED and also has effects on the other observables. If PAH emission is recorded for a source in the literature, the dust files which are identical to the above files are used, with the exception that they do not include the grain population with sizes smaller than 200 Å. These are instead calculated for using the dust distribution of Draine & Li (2007). PAHs were not detected for the source studied in Paper I, so the results of including these PAHs on our simulated observables are published here for the first time.

The interferometric observables are simulated in Fourier space. HOCHUNK-3D model images are generated at 1 μm intervals within the MIDI N-band range, the fast-Fourier transform is taken of the images, and the specific u-v point or visibility according to the position angle and baseline of the observations is extracted, corresponding to the data point obtained with MIDI. By applying this to each image, a visibility profile with wavelength is obtained. In order to simulate the Q-band images, model images are generated at the required wavelength and convolved to the telescope’s resolution. Azimuthally averaged radially profiles are then used to compare the model, source and PSF images. The SEDs were generated with IDL post-processing scripts (Whitney et al. 2013), which could then be compared with photometric data points compiled from the literature.

Fig. 1 Schematics of the suspected geometries of each of the sources with the position angles of each configuration of their MIDI data additionally shown. G305.20+0.21, the subject of Paper I, is also shown. The envelope is shown in blue, the cavity in yellow and the disk in red. Two different schematics are shown, one that applies to sources that have inclinations of less than 30 degrees. In this case, most of the cavity emission will be along the same region of the disk in the line of sight. The other schematic represents sources inclined to the point where cavity emission extends beyond the disk in the line of sight.

Fitting process

The fitting process is the same as that of Paper I, with the SED constituting the starting point for the fitting, before the high-resolution datasets were considered. The Ulrich-type envelope and cavities successfully used in de Wit et al. (2010) and Wheelwright et al. (2012) were the starting geometry for the envelope and cavity specifications in this work. The fitting process starts with the simplest geometry, which consists of only of an envelope. If a fit could not be obtained, then the model geometry becomes progressively more complex by changing outflow cavity parameters and adding disks and eventually disk substructure.

When combining a variety of observations, the balance between the fits of those observations needs to be considered. When the SED alone is utilised, a common method of uncertainty quantification has been the chi-square minimisation. The SED Fitter tool of Robitaille et al. (2007), which utilises a large grid of models and has been used to analyse the data of many protostars, employs this minimisation. Once a method is expanded past SED considerations, which as previously discussed is important to better constrain protostellar environments, the uncertainties in the fit parameters must be and have been considered differently. This is because different observations, at different wavelengths and using different techniques will all trace different scales of MYSO environments in contrasting ways. Few papers dedicated to MYSOs in the literature fit RT models to more than one set of observations, but the papers that are most similar to our approach have tackled this is different ways.

Linz et al. (2009) fit MIDI interferometric data and an SED for one MYSO. They used the aforementioned SED Fitter tool to constrain the fits on their SED and compared the model visibilities associated with a variety of acceptable models from that process in order to decide upon their best fitting model. A similar approach was taken in de Wit et al. (2010), where the SED Fitter was again used to provide some quantitative analysis of the SED fits, but the overall fit between the MIDI and SED was determined by eye. In de Wit et al. (2011), where one MIDI dataset and an SED were used, the overall fit was determined by eye. Johnston et al. (2013) modelled the SED and 2MASS image profiles of the MYSO AFGL 2591 using a Monte Carlo radiative transfer code. The reduced combined chi-square was calculated for the line profiles and their SEDs, giving equal weight to each observation within the final chi-square. They studied the constraints on their parameter space by systematically varying a few key parameters and investigating how the chi-squared varied with each parameter.

Boley et al. (2013) determined their best fits among 2D Gaussian models fit to their MIDI data by performing a grid search over their parameter space and the chi-square minimisation. The uncertainties on the fitting parameters were found using a Monte Carlo approach. This involved generating synthetic observations, which were normally distributed about the mean visibility of their observations, and repeating this to optimise. This was a viable approach given that only one observable (their MIDI data) was being used in their work. In Boley et al. (2016), a variety of different models were applied to a combined dataset of MIDI observations, T-ReCS aperture masking data and K-band interferometric data for MYSO. Gaussian modelling was again used to fit their visibilities, and the chi-square minimisation (as in Boley et al. 2013) was used to determine good fits; however a different fit found for the different wavelength datasets and no overall chi-square values were computed. Temperature-gradient disk models were also used to fit the K-band visibilities and the N-band correlated fluxes and a reduced chi-squared was calculated across the different observations; however, the weighting was not explicitly explained. Finally, radiative transfer models were employed via the use of the 2003 version of HO-CHUNK (Whitney et al. 2003). They did not attempt to reproduce the SED throughout this modelling; only the K-band visibilities and N-band fluxes were modelled. No errors are formalised on their final parameters, and the chi-square was not utilised during this portion of their analysis.

In Liu et al. (2019), the model grid of Zhang & Tan (2018) was used to fit the SEDs of a number of MYSOs. Following the method of De Buizer et al. (2017), they performed the chi-square minimisation to find five best fitting models for each of the MYSOs. Five parameters alone were varied within the fitting process in order to investigate how the SED is affected by models corresponding to different evolutionary phases that occur within the turbulent core model (Tan et al. 2014).

The aforementioned works constitute those most similar to the work on this sample, but more observations are used here and all are fit simultaneously. Hardly any of these works calculated a global chi-square for different datasets, nor did they derive formal errorbars on their parameters when RT codes were used. This appears to be the inherent compromise of combining many datasets and using a complex modelling approach. Including multiple observations and a large model parameter space provides a much better picture of a source, but formalising the combined error using these methods becomes arduous. Calculating an overall chi-square and giving each observation equal weight, as has been done in some of the works mentioned in the previous subsection, was not determined to be an appropriate approach for the work on our sample. The MIDI spectrum and VISIR profiles have dozens of points covering a very small wavelength regime, while the SED covers a wavelength regime from one to thousands of micron with a (most frequently) small number of data points. This means that there will be an obvious bias within any chi-squared value calculated for the combined datasets, in favour of the 7–13 μm and ~20 μm datasets. Another reason to not combine the chi-squared is that each observation traces very different scales so, similarly, the innermost and smallest geometrical features of the source, which are traced by the data-rich MIDI spectrum, will be unfairly focused on through the use of a combined chi-square.

Overall, it can be concluded that while the chi-squared may not be an appropriate way to provide a blanket estimate of a fit when multiple observations that trace different scales are used, it can assist in the determination of which observables are key for tracing different physical characteristics of MYSOs. In order to comment on the uncertainties involved in the fitting of the wide parameter space provided by the RT code, the chi-squared was calculated for each observation. The chi-squared was calculated separately for each configuration of the MIDI observations and then averaged. The image profile chi-squares were calculated from the radial profiles up to the point where the background noise exceeded the rms. The SED chi-squared was calculated by convolving the model SED with filter fluxes corresponding to the available data points. While the SED was the starting point for the fitting (as in de Wit et al. 2010; Linz et al. 2009) the SED fitter tool is not used. As such, the fits of these models avoid the inherent biases of the Robitaille grid of models, which are discussed at length in Robitaille (2017).

Following the above considerations and investigations into the individual chi-squares, we find that the different parameters can be traced within certain ranges. The disk inner radius could be tightly constrained within 5°, as could the cavity opening angles. Generally, significant changes were found in the MIDI fits for variations of the order 0.1 in scale-height in the two disk components. Similarly the cavity density exponent was constrained within ±0.1 and notable changes were observed in the observables when the cavity density was varied by a factor of two. Some parameter,s such as total disk mass, outer disk radius, and envelope outer radius, are not optimally traced by IR observables and as a result are poorly constrained. On average, the disk outer radius could usually be constrained to within 500 au, the disk mass could be constrained to an order of magnitude, and the envelope outer radius could be constrained within a factor of five. For sources that were particularly close or orientated such that the disk emission dominated or had a mass-dependent substructure such as a spiral, the disk mass could be better constrained, and we discuss these ranges on a source by source basis in the appendix.

4 Results

In this section, we present the results derived when applying the above methods to the sample of MYSOs. Table 2 describes summarises the key features of each MYSO, while Table 3 lists their specific parameters. An in-depth discussion of the fitting and features of each source is provided in Appendix A.

G305.20+0.21 is also included in the final sample but the results of its fitting have been previously published in Paper I. During the consideration of the full sample, we experimented with including the second disk component to see if improvements could be made to the fit of G305. Including the small grain disk component worsened the fit and, as such, the final results of Paper I are carried through to the comparative analysis made in Sect. 5. The final geometry for G305 is characterised by a disk with an inner radius of 60 au, over three times the sublimation radius of the central protostar. Similar to the pilot object, G305.20+0.21 (see Paper I), all of the sources studied with our methodology appear to share the same base geometry – a central object surrounded by a disk, bipolar outflow cavities and an envelope.

The central protostars of the models in this study cover a range of masses and luminosities. The protostellar masses were only varied from those obtained from the literature to improve the SED fits, as they had effects on the infall properties, and therefore the density distribution, of the envelopes.

The envelope properties of the MYSOs were traced predominately by the longer wavelength emission of the SED and the 20 μm emission, which is sensitive to the centrifugal radius (for closer sources) and the envelope infall rate. The larger the envelope infall rate, the larger the density within the infall radius (according to the Ulrich density distribution Ulrich 1976) and the larger the amount of material available for re-emission of the 20 μm flux. The infall rate reaches values of ~10⁻³M_⊙ yr⁻¹ for our sample. This is consistent with infall rates reported by McKee & Tan (2003) for the turbulent core model. Accretion rates of this magnitude are reported in numerical work (e.g. Hosokawa et al. 2010; Kuiper et al. 2011), which is consistent assuming that envelope infall rates translate to dynamical mass accretion rates onto the star. The envelope outer radii are all similar across the sample, but our high-resolution data are not sensitive to such large scales. All the MYSOs studied within our work have observables that can be fit with an Ulrich-type envelope, implying that all of the MYSOs are actively infalling and contain the subsequent density enhancements in the midplane that affect the 20 μm emission. Follow-up kinematic studies could confirm this to be the case following methods such as those detailed in Ohashi (2004) (e.g. searching for inverse P Cygni line profiles in the sources’ millimetre spectra). S255 IRS3 appears to be the exception in terms of envelope parameters, given its low envelope infall rate and smaller outer radius (by a factor of five). S255 is the only source confirmed to be episodically accreting so the fact that this source’s envelope infall rate is lower is particularly interesting, as it could be related to the accretion activity and the fact that the observations in this paper correspond to an accretion lull.

All the MYSOs in this sample have disks as a component of their final models. These disks are not toroids, which are 100 s of solar masses (Beltrán & de Wit 2016), but are objects that (while they could exist at larger radii) are present in the innermost (~100 au) regions of the protostellar environment and within the accreting regions. The fact that all of the sample harbour a disk, lends support to the idea that massive forming stars accrete their mass through similar mechanisms to low-mass stars.

The flaring and inner geometry of the disks had the greatest effects on the MIDI visibilities. For the majority of the model disks, the small grains disk component has a scale-height twice that of the large-grains disk. Such a ratio is in agreement with the hydrostatic calculations of D’Alessio et al. (1998) and models of observed images (Cotera et al. 2001; Wood et al. 2002). IRAS 17216-3801 has a flaring factor very close to the value most commonly observed throughout the sample of 2 at 1.75, but the fact that a binary system is likely to be having dynamical effects on this disk could explain this. Two of the remaining disks are flat (scale-heights of the disks were 0.1) and one of them required no small grains to be present in the disk to successfully fit the datasets.

Whitney et al. (2013) recommended that the distribution of the dust be split such that 20% of the grains are in disk 1 (the large-grains disk component) and the remainder are in disk 2 (the small grain component). This was based on the expected dust settling timescales observed for low-mass source. Even though the exact nature of dust amalgamation within massive protostellar environments is yet to be accurately constrained, we find that this approximation suited all of the observations well, save for those of G305. The fact that no disk 2 was required to fit the MIDI visibilities for this source was discussed as a significant element that characterises this object in a possible transition phase in its disk evolution towards full disk dispersion, as discussed in Paper I.

Inner clearing, as found for the disk of G305.20+0.21, was also found to be present in the disks of M8EIR, AFGL 2136 IRS1 and Mon R2 IRS2. The degree of clearing varied between the four sources, ranging between 2 and 4 times the sublimation radius. NGC 2264 IRS1 displays a different form of substructure. The spiral parameters of the code were utilised, and we find that including a spiral with a loose pitch angle of 30°, which contains 90% of the mass of the disk and begins and ends at 30 and 60 au, respectively, assisted in fitting the observables. These parameters essentially create substructure that manifests as a disk with a gap rather than an obvious spiral. NGC 2264 IRS1 is by far the closest MYSO of the sample and, as such, was traced at the smallest scales, so the fact that evidence is found for disk substructure only in the nearest of all sources in our sample is deemed consistent.

Of the disk parameters, one weaker constrained parameter is the disk mass. The only source where the disk mass was constrained to less than a factor of ten was NGC 2264 IRS1. This is likely due to the fact that it is near, that it has a close to face-on inclination, and that it has a mass-dependent substructure. The outer disk radius is only loosely constrained by the observables. A larger disk creates more millimetre emission, which can affect the longer wavelength regime of the SED and the maximum disk radius is also tied to the centrifugal radius, which affects the shape of the SED particularly at its 70 μm peak and around 20 μm. The imaging profiles also proved decent outer disk radii tracers, but mostly for the nearer sources. As such, the outer radius is still discussed but it is noted that this is not as constrained as the inner radius, which heavily affects the N-band emission in all cases.

5 Discussion

We have performed RT modelling of a sample of eight MYSOs invoking a star-disk-envelope configuration. In the following section, we discuss common features and traits identified in the sample. Detailed discussions of the individual sources can be found in Appendix A.

5.1 Comparing the MYSOs of the sample

Two factors are also defined to assist in the interpretation of the parameter space, resolving factor and flaring factor. The resolving factor was calculated to determine whether the resolving power of the MIDI observations affected the quantities derived. This was calculated by dividing the largest projected baseline (Vinković & Jurkić 2007) used to observe each source by the distance to the source. The flaring factor is simply the scale-height of the ‘small grains’ disk (disk 2) divided by the scale-height of the ‘large-grains’ disk (disk 1).

In order to make general comparisons between the sources in our sample, we calculate Pearson correlation coefficients. Correlations were inspected for basic parameters (bolometric luminosity, central object mass and distance) and the modelling parameters (cavity opening angle, disk flaring factor, envelope infall rate, cavity density, centrifugal radius and minimum and maximum disk radius). We do not calculate correlations for any parameters that were determined to be poorly constrained, namely disk mass and envelope radius.

Before considering the correlations they were checked for bias. To check whether the specific MIDI observations affected the correlations found, correlations between the resolving factor and the parameter space were calculated. Strong correlations were detected for two characteristics – the distance and stellar mass, which seem consistent, as a very distant source needs to be more massive or luminous to be detectable, but if a source is close and observed with small baseline lengths it may be comparably resolved as a source that is distant but observed with the largest baselines. We find that no significant correlations are found between the inclination and the rest of the parameter space and that therefore the results found for the sample are not a result of fortuitous orientations. We conclude that any strong correlations beyond these above are tangible results from the parameter space. The vast majority of the correlation coefficients are <0.5. This implies that while the sample all share the same base geometry, diversity still exists within the group. We discuss some correlations of particular interest in the next subsection.

Notable correlations

A strong correlation is detected between the scale-height of disk 1 and the cavity opening angle (0.72) implying that the more flared the disk, the larger the cavity opening angle (Fig. 2). If the disk scale-height is large, the disk is more flared and a larger emitting surface will be visible. This same affect occurs if the cavity opening angle increases as this in turn increases the available area for irradiation by the protostar and therefore generates a larger emitting region. The correlation between disk 2 scale-height and cavity opening angle is moderate (0.62) because one model does not possess this second disk component.

A correlation of 0.75 is found between the envelope infall rate and cavity opening angle, implying that as the envelope infall rate increases, the cavity opening angle will get larger (Fig. 3). The Ulrich geometry the envelope model uses is based on a rotating and infalling solution and results in a build of material at the equatorial plane. If the envelope infall rate is high, more material is moved from the outer envelope towards the midplane it becomes more concentrated while the regions closer to the poles of the envelope are less dense. The remaining envelope material would therefore be more vulnerable to the effects of any outflow activity from the central source, meaning that more material could be carved of the envelope, thereby widening the outflow cavity such that it has a larger opening angle, explaining why this trend is observed.

A strong correlation (0.92) exists between disk minimum radius and luminosity (Fig. 4). If considering straightforward dust sublimation this is to be expected, yet the minimum dust radius does not correspond to the actual dust sublimation for many sources. This is illustrated in Fig. 5. In Paper I causes of inner holes in disks were discussed atlength and it was determined that the most likely causes for inner holes in massive stars are (1) photoevaporation and (2) the presence of binary/multiple companions. Of these two causes, photoevaporation as a mechanism seems to be in agreement with this trend. If a source is more luminous, then it will have a greater capacity to generate photons of high enough energies to dissociate the inner disk material and create cleared inner regions. Therefore, the large luminosity leading to a large inner hole makes sense as the star has another mechanism through which to destroy the disk. This comes with the caveat, however, that the star must not also be creating UV photons to the point where an ionising ultra-compact HII region has begun to form around the source, as that is not observed for any of these sources. A strong correlation does not exist between minimum envelope radius and luminosity, as the minimum envelope radius of the final model for IRAS 17216-3801 vastly differs from its minimum disk radius. If IRAS 17216-3801 is removed, the correlation between minimum envelope radius and luminosity is becomes strong and the correlations between luminosity and both disk and envelope minimum radius become 0.93.

Fig. 2 Relationship between the scale-height of the large-grains disk component (D1 SH) and the cavity opening angle.

Fig. 3 Relationship between envelope infall rate and cavity opening angle.

Fig. 4 Relationship between the luminosity of the source and the minimum disk radius.

Fig. 5 Relation between luminosity and sublimation radius (orange dots). For comparison, the blue dots denote the disk minimum radii as in Fig. 4.

5.2 Substructure in MYSO disks

In our sample of eight MYSOs, selected merely on the basis of completeness of observations, our analysis indicates that five sources have evidence for substructure in their disks. A disk with substructure in our terminology is contrasted with an azimuthally symmetric disk extending inwards down to the dust sublimation radius. For all except one, this takes the form of a cleared inner hole. The selection of objects based on completeness of IR spatial and SED data impliesa number of things for our sample compared to massive stars in general. Firstly, the objects are relatively close by and only G305.20+0.21 lies at a distance greater than the typical 2–3 kpc distance expected for MYSOs. Secondly, they are statistically expected to be at the lower end of the massive stellar mass distribution. Additionally, they possess somewhat wide outflow cones, such that a significant portion of the IR radiation from the inner disk regions can escape. Thus, finding objects with evidence of ‘transition disks’ amongst this sample would not be totally unexpected. The presence of the inner holes tentatively suggests that within the idea of dissipation of the circumstellar environment, the equatorial disk structure is sensitive to advanced erosion prior to a complete disappearance of the overlying envelope. The likely causes of this inner clearing were discussed in Paper I, where it was suggested that the presence of a companion or photoevaporation are the most likely causes for the clearing in MYSOs. In order to confirm which of these mechanisms is responsible, high-resolution imaging with the ELT/METIS or a much denser filling of the uv-plane with VLTI/MATISSE could be employed to investigate dust sublimation and disk evaporation and how this relates to the physical presence of an inner disk rim.

The remaining source, NGC 2264 IRS1, showed substructure in the form of a gap-like spiral feature. The fact that these substructures are present has interesting implications for the massive star formation process. To date, substructures have rarely been reported around massive protostars. Sanna et al. (2019) detect a sub-Keplerian disk and jet system around the MYSO G023.01-00.41 and find that the CH₃OH emission at scales ~2000 au is significantly warped with respect to the plane. Maud et al. (2019) detect a ring-like enhancement in the disk of AFGL 2136, whose size scale matches the inner rim in our final model (see Appendix A for further details). Various other disks (e.g. Ilee et al. 2018)show irregular features as a result of fragmentation. The substructure that improves the fit of the data for NGC 2264 IRS1 manifests as a gap in the disk. Various physical processes have been suggested that result in a well-defined radial discontinuity within a disk with a significantly lower density or even devoid of material. The first is the presence of a forming planet. Pinte et al. (2019) present one of the first confirmations of a planet in a gas and dust gap in the protostellar disk of the 2.4 M_⊙ object HD 97048, but otherwise the presence of planets in these gaps has mostly been inferred. As discussed in the context of G305 in Paper I, planets are not observed around massive stars but whether they can never form or cease to exist as the protostar becomes more luminous with age has not been confirmed. An alternative method of gap creation is dust trapping. Dust traps (Whipple 1972) occur when a local gas pressure maximum in the disk attracts andtraps dust grains allowing clearing to occur in the disks and has also been proposed as a mechanism for observedgaps in low-mass disks (e.g. Gonzalez et al. 2017). Recently van der Marel et al. (2018) illustrate through physical-chemical modelling that the gaps observed in ALMA data of the disk of the Herbig Ae star HD 163296 can be reproduced from dust trapping alone, without the presence of a planetary body. Dust growth presents another avenue of gap formation whereby the amalgamation or pile-up of dust grains clears space in the disk (e.g. Zhang et al. 2015; Okuzumi et al. 2016). Instabilities, both large-scale (Lorén-Aguilar & Bate 2016) and gravitational (Takahashi & Inutsuka 2016), have also been shown through modelling to be able to cause gaps. Among these processes, the ones that would operate in a disk around a massive star are probably grain growth and instabilities, provided that the young massive star accretes from a high density (and massive) disk. Toroids, which have been shown to create instabilities in works such as Lorén-Aguilar & Bate (2016), have been observed around massive forming stars (such as those studied in Beltrán et al. 2005), presenting an avenue for this instability and therefore gap to occur. Gravitational instability is a more likely occurrence in massive systems where the central protostellar mass is higher and when multiple systems are present. Price et al. (2018) have shown through hydrodynamical simulations that spirals, shadows and asymmetrical structures can all be produced by a companion star on an inclined and eccentric orbit approaching periastron. Both the effects of binaries and dust growth take 10 s of thousands of years to occur, and are therefore also better suited to an MYSO in a later phase of formation.

6 Conclusions

In this paper we have presented a multi-scale and multi-wavelength analysis of a sample of eight MYSOs with the aim to make direct comparisons of their characteristics. We derive physical parameters from the observational data by means of radiative transfer modelling, fitting simultaneously spatially resolved and SED data. The RT model consists of a disk-outflow-envelope morphology and we find that such a geometry provides reasonable fits to the data. Envelopes were best described by Ulrich-type density distributions, suggesting that they could be dynamically unstable and implying that they are actively infalling environments.

Agreement exists between the characteristics found through this work and those found through previous studies. From the study of this sample it appears that the ~10 μm interferometric data (probing ≾100 au scales) are most sensitive to the inner disk geometry of the sources. This could be determined thanks to the combination of the interferometric data with the ~20 μm single-dish data which provides additional constraint on cavity material. Variation is seen at the small scales (≾100 au), in particular within the disk geometry. Substructure is detected for 75% of the sources, that is the disks deviate from an azimuthally symmetric disk extending inwards down to the dust sublimation radius. The most common form of this substructure is an inner hole, implying some dispersal mechanism is present within the inner environments of the MYSOs in addition to dust sublimation. A strong trend observed between the minimum dust radius and the luminosity of the central source implies that photoevaporation may be a likely candidate for the observed clearing. The final source with substructure, NGC 2264 IRS1, displays a gap-like structure. Given that increasing the fraction of large dust grains in the disk improved the fit to the N-band data for this source, and the high multiplicity fraction for MYSOs (Pomohaci et al. 2019), we suggest that dust growth and the presence of a binary/multiple system are the most likely candidates for any substructure in this disk. All the inferred substructures have been previously observed in evolved low-mass protostellar disks. This implies that despite their shorter evolutionary timescales, massive disk evolution follows a similar path as low-mass disk evolution. Future work to detect multiplicity and further improve the view of these inner regions through spectrointerferometric, imaging and image reconstruction methods, will assist in confirming these substructures and distinguishing their origins.

Acknowledgements

We thank the anonymous referee for their helpful pointers and Prof. Melvin Hoare and Prof. Steven Longmore for their discussions, all of which enabled the improvement of this work. We also thank the STFC for funding this PhD project. Themodelling presented in this work was undertaken using ARC3, part of the High Performance Computing facilities at the University of Leeds, UK.

Appendix A In-depth presentations of each source

This section is dedicated to each object. For each object we (1) provide an introduction on the object itself (2) discuss any interesting effects on each of the simulated scales when experimenting with the parameter space of the RT code, and (3) present the results and final best fitting model for each object.

Fig. A.1 Observed visibilities for W33A for each configuration (black) with the simulated visibilities (coloured) for different models to show the effects different components on the simulated MIDI data.

G305.20+0.21 is also included in the final sample but the results of its fitting have been previously published in Paper I. During the consideration of the full sample, we experimented with including the second disk component to see if improvements could be made to the fit of G305. Including the small grain disk component worsened the fit, and as such the final results of Paper I are carried through to the comparative analysis made in Sect. 5. The final geometry for G305 is characterised by a disk with an inner radius of 60 au, over three times the sublimation radius of the central protostar.

A.1 W33A

de Wit et al. (2010) investigated the MYSO W33A by fitting an RT model to the SED and MIDI visibilities. They found that the source’s MIDI data and SED could be reproduced by a model consisting of an envelope geometry with two bipolar outflow cavities and concluded that the majority of the mid-IR emission originates from the outflow cavity walls, in agreement with other works such as De Buizer (2006). Ilee et al. (2013) attribute the CO bandhead emission of W33A to a 1–2 au disk, while Davies et al. (2010) find a Brγ jet perpendicular to a disk found in CO emission and absorption using NIR adaptive optics integral field spectroscopy. At longer wavelengths, two large continuum sources are observed within W3, the brighter of which is MM1 (Galván-Madrid et al. 2010) or W33A. Maud et al. (2017) studied W33A with ALMA and SMA data to investigate whether a Keplerian disk may be present. They detected a change in position of the blue-shifted and red-shifted emission close to the peak of the continuum emission, which could be indicative of a rotating disk. They probed complex large-scale structure, with hints of a disk at smaller scales. The observations used in our study of the source consist of the MIDI observations from de Wit et al. (2010), the COMICS image presented in Wheelwright et al. (2012) and an SED compiled from the RMS database and de Wit et al. (2010).

The luminosity of the source is listed as 3.4 × 10⁴ L_⊙ by the RMS survey (Mottram et al. 2011). The kinematic distance to W33A was thought to be 3.8 kpc (Faúndez et al. 2004) but water maser measurements by Immer et al. (2013) recalculate this as 2.4 ± 0.2 kpc. Our work therefore updates the mid-IR view of W33A substantially, updating it for the correct distance, adding to previous methodologies by including imaging data as an aspect of the fitting process and by using a more sophisticated radiative transfer code and due to the handling of the interferometric field of view (as discussed in Paper I). A test was done to see if the fit could be reproduced by simply reducing the luminosity of the central source according to the change in the distance. When done, the simulated visibilities for all the MIDI configurations increased between 8.5–10.5 μm and the visibility at 7.5 μm for configuration i), taking them outside of the observed errorbars and worsening the fits. This may be because the sublimation radius of the source will have become smaller, presenting a more compact component for emission resulting in increased visibilities. The 24.5 μm density profile was also too broad and the silicate absorption feature of the SED became deeper and in violation of the observed MIDI spectra.

After determining that changing the distance alone could not provide a suitable fit to the dataset, recent literature, notably the existence of a potential disk (Maud et al. 2017), was taken into account when improving the fit. Maud et al. (2017) studied W33A with ALMA, tracing its wider environment. In their work they note the turbulent and messy nature of the source and image a stream of material that is spiral-like in nature on scales of ~5000 au. Although Keplerian signatures at the centre of the source are weak, they also state that while they do not detect a disk on scales of the spiral, a disk component could be viable at scales ≤500 au.

Including a disk within our models was experimented with, and improvements were found to the fits of the MIDI data. A disk can be included in the model with a maximum radius of 500 au, in agreement with the suspected scales of Maud et al. (2017). The disk can be made smaller than this without violating the fits but given the size of the protostellar environment and the suspected mass of the source, 500 au is used as the final maximum disk value. This implies that the disk’s emission from its inner regions is more important than the emission from its outer regions, which is consistent as the outer regions of a dense feature like a disk will be too cool to be traced by our high-resolution IR data. The disk is flat, with the large-grains disk component scale-height set to half that of the smaller grains disk. The model envelope has outflow cavities carved out at a 20° opening angle, twice that of de Wit et al. (2010), but with a similar cavity density. The envelope infall rates are identical between the two works. The inclination from de Wit et al. (2010) for the 3.8 kpc distance well satisfies all observables for the closer distance. The visibilities are displayed in Fig. A.1, the model 24.5 μm image and the radial profiles are displayed in Fig. A.2 and the SED is displayed in Appendix A.3.

To summarise, W33A has demonstrated the importance of having the correct distance to MYSOs when deriving geometrical traits with the final model for the new distance characterised by the presence of a flat disk, which was not a part of the previous 3.8 kpc model of de Wit et al. (2010).

A.2 NGC 2264 IRS1

NGC 2264 IRS1 presented an interesting case for the fitting process due to its close distance of ~740 pc. During the fitting of other objects at more typical MYSO distances (2–4 kpc) the COMICS profile was largely unaffected by the variation of disk parameters. In the case of NGC 2264 IRS1 however, the disk had profound effects on the profile and added another level of complexity to the fitting. Through the fitting it was determined that an Ulrich-type density distribution, as opposed to the power-law one of de Wit et al. (2009), results in a model that produces a satisfactory fit to the datasets. Including the envelope is necessary to obtain a good fit to the longer wavelength data in the SED, which was not reproduced by previous work by Grellmann et al. (2011). The final model consists of adisk-outflow system, surrounded by this Ulrich envelope. The central protostar has a luminosity of 4.17 × 10³ L_⊙ which is in good agreement with those found in the literature (Harvey et al. 1977, Schwartz et al. 1985). The stellar mass in the model is 9 M_⊙, similar to the mass postulated by Thompson et al. (1998) based on the source’s B0-B2 spectral type. An inclination of 15° was able to successfully fit all the observables, which is consistent with the existence of an outflow in the line of sight found by Schreyer et al. (1997).

Fig. A.2 Convolved model image (top, contours represent 5, 10, 25 and 75% of the peak flux) and the radial profile (bottom) of the best fitting model for the 2.4 kpc distance. The original COMICS image of W33A can be found in Wheelwright et al. (2012).

Fig. A.3 SED of the preferred final model for W33A. Data from the RMS survey are shown as diamonds, the MIDI spectrum is shown as red crosses and data points from de Wit et al. (2010) are shown as triangles.

Figure A.4 shows the MIDI data and two models which are discussed further later in this section. The twoN-band configurations of NGC 2264 IRS1 trace significantly different areas of the protostellar environment, with large differences in both PA and baseline (an illustration is provided in Fig. 1). Both the MIDI configurations predominately trace disk material. The 40.2 m baseline was characterised by its upward trend in visibility which proved challenging to fit. Often the issue was, again, that the model visibilities were too high between 7.5–9 μm especially, indicating that the emitting regions at larger scales within the disk were also too compact and needed to be minimised. A near-perfect fit for Configuration (ii) could be found when a disk was included with a 500 au maximum radius consisting of 20% large-grains (Model 1 in Wood et al. 2002) and 80% small grains (Galactic ISM grains as defined by Kim et al. 1994) but the fit for the 7–9 μm region of the 40.2 m baseline fit was poor in comparison. By adding substructure to the disk by using the 3D, spiral capabilities of the code, the fit for configuration (i) could be improved (as described in Sect. 4). Including this ‘spiral’ dramatically improves the fit of Configuration (i) but worsens the fit of configuration (ii) as can be seen in Fig. A.4. The overall change in the simulated visibilities for configuration (ii) (V ² ~ 0.02) is smaller than the improvement in visibility for configuration (i) (V ² ~ 0.08) so the model with substructure is referred to as the preferred fit.

Fig. A.4 Observed visibilities of NGC 2264 IRS1 for each MIDI configuration (black) with the simulated visibilities for each model image (coloured).

Given there is some evidence for PAH emission in the literature (e.g. Gibb et al. 2004), the inclusion of PAH emission in the code was experimented with. Adding PAHs/VSGs to all areas of the model protostellar environment increases the visibilities overall, vastly worsening the fit of configuration (i). However, including the PAHs in just the cavity and envelope results in no change to the MIDI visibilities. This consecrates the fact that the MIDI data are dominated by the disk for this source. Ultimately a compromise is made by adding PAHs to the cavity and envelope only. This can be explained physically as PAHs constitute some of the smallest components of protostellar material and behave like very small dust grains. As a result, they could be vulnerable to the effects of photoevaporation. This means that they could be pushed from the inner environment into the envelope and cavity material at greater distances from the central protostar. Additionally, the PAHs in the further reaches of the envelope material would feel a weaker force from the stellar photons and would not be so easily removed from the environment. Therefore, the idea that PAHs may be present in outer regions of the protostellar environment but not the inner regions, may be consistent (Siebenmorgen & Krügel 2010). The improvements noted by the inclusion of the spiral exist with or without the presence of PAHs.

Figure A.5 shows the COMICS image, model 24.5 μm image and the resulting radial profiles. The image is symmetrical and appears to show details of substructure, but de Wit et al. (2009) state that most of these are due to diffraction effects. Given its close distance, the 24.5 μm emission of NGC 2264 IRS1 is heavily influenced by the disk included in the model. The close distance, inclination and the need for substructure within the model meant that the disk mass for NGC 2264 IRS1 could be more accurately constrained than other sources to an order of ~0.1 M_⊙. Changing the cavity density exponent, and therefore how evenly the cavity material was distributed, could assist in the fitting of the source. Increasing the exponent from the default value of 0 (constant density throughout the cavity) to 0.25 (meaning the mass was slightly weighted towards the bottom of the cavity, close to the protostar), reduced the amount of 24.5 μm emission creating a better fit. This is consistent as more emitting material is at smaller radii within the protostellar environment, making the source appear less extended. A combination of the close-to-pole-on inclination and a foreground A_v of 25 allowed sufficient fitting of the silicate absorption feature in the SED and is close to A_v values quoted in the literature (Thompson et al. 1998; Thompson & Tokunaga 1978). Changing the cavity exponent to satisfy the COMICS profile flattened the peak of the SED slightly, improving the fit. The SED fit is shown in Fig. A.6.

A.3 S255 IRS3

S255 IRS3 has become of particular interest since an accretion-ejection event was observed around the source by Caratti o Garatti et al. (2017). A particularly bright 6.7 GHz methanol flare was discussed in Fujisawa et al. (2015) and a sub-millimetre flare was also found in Zinchenko et al. (2017). This flare was later tied to NIR and radio flares by Caratti o Garatti et al. (2017) and these were determined to be the result of an accretion burst of order 10⁻³ M_⊙ yr⁻¹. This was the first directly observed accretion burst around an MYSO, which has important implications as episodic accretion is a key component of the low-mass star formation process (Stamatellos et al. 2011).

All the interferometric data for S255 IRS3 came from Boley et al. (2013), with six configurations of data available for fitting. All the configurations, save one, trace very similar position angles of around 70°. The remaining configuration traced a position angle of 32.6° and was low resolution due to its 17.6 m baseline.

The COMICS image from de Wit et al. (2009) is used in the fitting process. Both IRS1 and IRS3 are visible in the COMICS image. In order to exclude scales that may be influenced by the emission of IRS1 we fit the inner 3″ of the source as opposed to the whole profile. All of the fitted data were taken at a time that appears to be an accretion lull for S255 IRS3, as thefluxes from the SED (which were taken at similar times to the MIDI and COMICS data) are in much better agreement with the non-burst fluxes from Caratti o Garatti et al. (2017) than the burst fluxes. Fluxes from the RMS surveyand Caratti o Garatti et al. (2017) were incorporated (blue diamonds in Fig. A.9), alongside additional fluxes from de Wit et al. (2009) (green crosses in Fig. A.9) and the MIDI flux spectrum.

Fig. A.5 COMICS 24.5μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.

Fig. A.6 Model SED of the best-fitting model (black). Multi-wavelength flux measurements from the RMS are represented as blue diamonds, the yellow diamond represents the COMICS flux density and the fluxes corresponding to the MIDI visibilities are also shown in red.

The plots for the final model are shown in Figs. A.7–A.9. The final luminosity for the central source in the model is 2.15 × 10⁴ L_⊙, which is in excellent agreement with that observed by Wang et al. (2011) of ~2.4 × 10⁴ L_⊙. The overall model is a disk and outflow system, in agreement with works such as Zinchenko et al. (2017), who studied the source with SMA and IRAM-30 m continuum and spectral lines and find signatures of rotation and attribute this to a disk-like structure. The presence of the disk also agrees with the detection of the accretion event; as for low-mass stars, such events occur as a result of accretion from a disk surrounding the protostar. The fact that the envelope is rotating and infalling could be linked to the fact that the source has displayed an accretion event. Infall from the wider environment is directly linked to disk accretion in the turbulent core model (Tan et al. 2014). With a centrifugal radius of 500 au, as opposed to 2000 au, more material is fed closer to the central star, providing a larger reservoir of easily accessible material for accretion. This need not be a constant flow of material, meaning that the envelope infall rate and disk accretion rate can fluctuate throughout the life of the stellar system, allowing for low accretion rates. A flatdisk with no inner clearing allows a good fit to be obtained for the MIDI visibilities of S255 IRS3. Two disk components were included in the successful model, a small grains disk and a large-grains disk, both with the same scale-height. Configuration (vi) is the worst fit of the set. When consulting Fig. 1 of Caratti o Garatti et al. (2017), we can see that this configuration is still likely to be tracing cavity material, as illustrated by Fig. 1. The simulated visibilities are lower than those observed implying that the smooth cavity surface presents toolarge an area for emission to match the data. This could be explained by the fact that the cavity material is clumpy and irregular (Caratti o Garatti et al. 2017). The cavity density of S255 IRS3 is of order 10⁻¹⁹g cm⁻³, which is comparable to that of W33A.

Fig. A.7 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured).

An inclination of 120° fit all the observables, which is consistent with the close-to-edge-on orientation suspected by Boley et al. (2013). The foreground extinction towards the source islarge (A_v ~70) similar to the findings of Caratti o Garatti et al. (2017). Given the close to edge-on inclination, a disk mass from 0.1–10 M_⊙ could be present and satisfy the observables which is poorly constrained due to the limitations of the IR as a disk mass tracer. A compromise was found for the longer wavelength regime of the SED between the data of Caratti o Garatti et al. (2017) and the data from de Wit et al. (2009) (that they do not list as upper limits).

Varying the disk parameter space has negligible effects on the COMICS profile, save for the disk outer radius and the centrifugal radius of the envelope, which is attributed to the sources inclination and distance. It was found that the 2000 au maximum disk radius and centrifugal radius used for G305 resulted in too much 24.5 μm emission to successfully fit the radial profile, implying that the disk must be less extended. Reducing these values to 500 au allowed a good fit for the COMICS profile.

The wider environment of S255 IRS3 differs compared to the other MYSOs in this sample, with an outer envelope radius that is five times smaller than the other sources and an envelope infall rate that is significantly lower. Raising the infall rate to a value comparable to those found for the other sources (~10⁻⁴ M_⊙ yr⁻¹) increases the peak of the SED significantly. While this appears to improve the fit to the SED, the fluxes that are better fit are almost certainly detecting flux from both IRS1 and IRS3 and as a result this envelope infall rate is likely to be an overestimation of what is actually present. A fit with the lower envelope infall rate that fits the smaller aperture ~70 μm point better is preferred. High-resolution, small aperture observations at ~100 μm could allow separate SEDs for the two MYSOs to be obtained and for this to be better investigated. While obtaining and analysing these longer wavelength data are outside the scope of this work, increasing the envelope infall rate also lowers all the simulated MIDI visibilities that only trace IRS3, not IRS1. Using an envelope infall rate of the same order of magnitude as W33A and NGC 2264 IRS1 vastly worsens the fits for all configurations, further showing the need for a lower envelope infall rate.

A.4 IRAS 17216-3801

Kraus et al. (2017) resolved IRAS 17216-3801 into a close binary using spectro-interferometric data from AMBER, CRIRES and GRAVITY. They find that the two stars are separated by ~170 au and that not only is there a disk surrounding the pair of stars but disks surrounding each individual object as well. HO-CHUNK3D does not provide the capabilities for modelling multiple systems but does present the opportunity to investigate the circumbinary material of IRAS 17216-3801. A single central object is included in the model, whose luminosity was obtained through the fitting of the SED as per the other sources in this work.

Fig. A.8 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.

Fig. A.9 Model SED of the best-fitting model (black) for S255 IRS3. The different symbols correspond to different datasets, purple crosses are data from Caratti o Garatti et al. (2017), blue diamonds are fluxes from the RMS database, squares are data from de Wit et al. (2009) (with the COMICS flux in yellow) and the fluxes corresponding to the MIDI visibilities are also shown in red. Open symbols correspond to upper limits.

Surrounding this central source is a disk-outflow-envelope system. The disk scale-heights are very similar to the defaults suggested by Whitney et al. (2013), with the small grains disk being slightly flatter. The inner disk radius is 100 au, which is comparable to the estimate of Kraus et al. (2017) for their circumbinary disk. The visibilities are shown in Figs. A.10 and A.11. IRAS 17216-3801 has the largest MIDI dataset of our sample with ten configurations. The position angles of this dataset range from approximately 14–175° but all are fairly low resolution with baselines of 35 m or shorter. This means that most of the MIDI configurations are sensitive only to the circumbinary material. As illustrated in Fig. 1, configurations (viii), (ix) and (x) are the only three with PAs suitable for tracing the circumprimary disk and only configuration (viii) has the potential to trace the circumsecondary disk. The baselines of configurations (viii), (ix) and (x) are some of the largest for this source; configuration viii) can probe scales of ~20–35 mas, configuration (ix) traces scales between ~20 and 50 mas and configuration (x) traces scales between 19 and 33 mas. As such these configurations not only have the PA to trace the circumstellar disks, but also the required resolution. Configurations (vi) and (vii) have comparable resolving power to configurations (viii), (ix) and (x) but their PAs mean that they are more likely to be tracing any material within the cavity carved out by the binary. The baselines of the remaining configurations are short, meaning that they would all be tracing the suspected circumbinary disk material. The average visibility across each dataset varies considerably, and no trend is observed between the average visibility of the configurations that may be tracing circumstellar disk material and those that will not. It was found that changing the minimum envelope radius to values smaller than 100 au improved the fits of configurations (viii), (ix) and (x) by lowering the simulated visibilities across each configuration. As previously mentioned each of these configurations has the potential to be tracing circumstellar disk material, so the fact that including additional dust within the circumbinary disk improves the fits to these configurations appears consistent. If the configurations are not tracing circumstellar disk material, having a smaller minimum envelope radius than the minimum disk radius implies thatthe binary stars cannot have destroyed the dust in the inner hole nor accreted it into their binary disks.

Fig. A.10 Observed visibilities of IRAS 17216-3801 for each configuration (black) with the simulated visibilities for each model image (coloured) (1).

In order to satisfy the VISIR data (Fig. A.12) a low inclination was required (15°) and a cavitydensity exponent of 1. This is because the VISIR profile is barely resolved. A density exponent of 1 as opposed to 0 (the default value) means that the density in the cavity decreases linearly with radius instead of remaining at a constant value. The innermost part of the cavity will then be brighter than its outer regions, reducing the extended emission and lowering the ~20 μm flux at larger radii, improving the fit. Similarly, close to pole-on inclinations reduce the extended emission as only one of the cavity lobes is directly visible in the line of sight.

The SED is shown in Fig. A.13. The ATLASGAL 870 μm point proved difficult to fit. The aperture of this measurement was 19.2″ but the region is fairly uncrowded with no obvious sources contaminating such a measurement. In an attempt to improve the fit to this data point, disks of bigger radius were trialled. The larger the extent of the disk, the more cold emission at sub-millimetre and millimetre wavelengths will be present. A 4000 au maximum radius disk showed small improvements in the millimetre of the regions but more substantial changes in the 20 μm region of the SED, making the dip-like feature between 10–20 μm flatten. This increase in radius caused all the simulated MIDI visibilities to decrease, worsening some fits and but improving the fit of configuration (ix). The fit of the VISIR profile was also violated as the model radial profile displayed more extended emission than the observed profile. This occurred for any disks larger than 2000 au so this was the final value chosen for the disk to preserve the quality of the VISIR fit. Making the envelope radius larger improved the fit to the longer wavelength end of the SED by shifting the SED to longer wavelengths, but this worsened the fit of the 20 μm region of the SED.

The mass of the central source was also found to affect the shape of 10–20 μm region of the SED. For example, if a 20 solar mass central source is used this 10–20 μm region of the SED is flatter. This is consistent as the density distribution of the Ulrich envelope is proportional to $M_{*}^{1/2}$ (Terebey et al. 1984), meaning that as the stellar mass increases, the density will decrease within the infall radius and reduce the amount of 20 μm flux. Configuration (v) proved the most difficult to fit. This occurred despite the fact that it is one of the lowest-resolution datasets and is likely to be tracing material within the circumbinary disk. A large (nearly eight-fold) upturn is present in the visibilities between 7 and 9 μm, which is often a hallmark feature of a compact component like a disk. To test whether the very large increase in visibility could be due to the presence of the inner rim specifically, the inner disk radius was changed to the disk to 134 au, which is the scale traced by this configuration at 7.5 μm. A minor upturn in the visibilities was observed but the fits for the other, higher-resolution configurations worsened. When the inner radius of the disk was varied between 100 and 150 au the visibilities changed marginally but the overall goodness of fit changed negligibly.

Fig. A.11 Observed visibilities of IRAS 17216-3801 for each configuration (black) with the simulated visibilities for each model image (coloured) (2).

The upturn in configuration (v) could be due to the binary system interacting with the circumbinary disk. Small-scale features, such as a puffed-up inner rim and disk-spiral, were tested but also present marginal changes, mostly likely due to the object’s distance. A very specific local substructure must be present at the scales traced by this baseline alone and not the others of similar PA/baseline to induce the observed increase in visibility. Various theoretical works have shown that binary or multiple stellar systems can result in the formation of such a small-scale substructure in protostellar disks. Nelson & Marzari (2016) used the SPH code VINE to model the circumbinary disk of the low-mass source GG Tau A. They found that the binary system ‘stirs’ the materialin the disk, leading to the formation of spiral structures, which extend towards the centre within the inner radius of the disk as streams in the earliest stages. Within the disk itself, they see significant areas of over-density. Thun et al. (2017) find that the inner rims of the circumbinary disks are often unstable, clumpy and recessing from their hydrodynamical modelling. Desai et al. (2019) focus on the study of unequal mass binaries in particular using a 3D hydrodynamics code CHYMERA. Here they find that gravitational instabilities are induced by the binary interaction, with the smaller of the companions creating a density wave at the earliest stages and at later stages also see the generation of spirals. Observationally, Tobin et al. (2016) detect spiral structures surrounding a three-star system. Although low mass, two of the stars have a separation comparable to that of IRAS 17216-3801 of ~180 au. These authors postulate that this distant star is the cause of disk fragmentation in the system and that the disk could have recently undergone a gravitational instability, inducing the formation of one or two of the companion stars. While all these cases focus on low-mass stars, the potential a binary system has to create gravitational instabilities and affect the stability and geometry of its circumbinary disk will be increased if those stars have higher mass. As such, the likelihood of over-densities that could create the compact feature required to fit this configuration is high and could therefore explain the poor fitting of these data.

A.5 Mon R2 IRS2

The Mon R2 region is home to five different sources and located at a distance of ~840 pc (Herbst & Racine 1976). A large extended HII region, illuminated mostly by scattered light, covers 30″ of the region and is visible in both 870 μm maps (Giannakopoulou et al. 1997) and the 24.5 μm image of de Wit et al. (2009). Mon R2s IRS2 lies within the centre-north region of this HII region and is its main source of illumination (Aspin & Walther 1990). A large CO bipolar outflow is present in the wider environment (Giannakopoulou et al. 1997) but it is unknown which MYSO in the region is responsible for it. Another massive protostellar system, Mon R2 IRS3, which has been resolved into two sources A and B, is the dominant source of the whole region. No free–free emission is associated with IRS2 and it is unresolved in UKIRT observations (Alvarez et al. 2004). Boley et al. (2013) find its mid-IR emission to be compact (~40 au), and attribute this to its lower luminosity than IRS3A and B. Henning et al. (1992) find a luminosity of the source of ~6000 L_⊙. Recent work by Jiménez-Serra et al. (2020) studies the gas of Mon R2 IRS2. They find that their ALMA images, studies of the H21α emission and 3D radiative transfer modelling show the presence of an ionised Keplerian disk and a neutral rotating disk. The disk is highly flared (as shown in Fig. 4 of the cited paper) and deviation of the H21α emission centroids from the disk plane imply warping in the disk, which they suggest could be due to the presence of a secondary object.

Fig. A.12 VISIR 19.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.

Fig. A.13 Model SED of the best-fitting model (black) of IRAS 17216-3801. Multi-wavelength flux measurements from the literature are represented as blue diamonds, the yellow diamond represents the VISIR flux density and the fluxes corresponding to the MIDI visibilities are also shown in red.

Our work finds that Mon R2 IRS2 is another source where the fits to the MIDI data are improved by the inclusion of an inner hole (20 au in radius). The dust sublimation radius of the source is small (~6 au) due to comparatively low source luminosity required to fit the other observables. The PAH emission lines for the source are weak (Gibb et al. 2004), and it was found that the MIDI fit could be improved by removing the PAHs from the disk. Removing the PAHs in theenvelope and cavity material had negligible effects on the MIDI fits, implying that the N-band emission of the source is disk-dominated. Removing PAHs from the entire environment greatly worsens the NIR slope of the SED fit. The final luminosity of 5520 L_⊙ is similar to that of Henning et al. (1992). The envelope of Mon R2 IRS2 is similar to that of the majority of the sample, with a comparable infall rate and the same outer radius. The cavity density is low (~10⁻²¹ g cm⁻³) and the cavity opening angle is 30°. The source is inclined at 130° so both cavity lobes are visible. Through the SED fitting it was difficult to simultaneously fit the silicate absorption feature and the 20 μm flux. Eventually the fitting of the silicate absorption feature was prioritised as this is the dataset with higher spatial resolution. The visibilities in Fig. A.14, the images in Fig. A.15 and the SED in Fig. A.16.

Fig. A.14 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for Mon R2 IRS2.

The position angles of the MIDI baselines for Mon R2 IRS2 mean that all five configurations are likely to be tracing disk material and configuration (i) is likely to be tracing cavity material. Three of the MIDI configurations have high visibilities and are tracing an unresolved component of the protostellar environment. The simulated visibilities of configuration (i) get lower if the outer radiusof the disk is made larger implying that a small, compact disk best satisfies this configuration. However, as R_c /$R_{\text{max}}^{\text{disk}}$ increased to values of 4000 au (in agreement with some of the largest observed disks around MYSOs e.g. Johnston et al. 2015) improvements were seen to the fit of the COMICS profile as the flux at larger radii increased, so a compromise ultimately had to be made. Configurations (iv) and (v) of the MIDI data proved the most difficult to fit and resulted in the inclusion of the hole. Without this the visibilities for both the configurations were too high. This hole is 40 au across, in agreement with the size of the compact emission from Boley et al. (2013). This implies that the mid-IR emission for this source, as traced by their work, is mostly dominated by emission from this inner rim. A compromise was found between thefitting of configurations (iv) and (v) as opposed to fitting one perfectly and leaving the other poorly fit.

A.6 M8EIR

This work finds that the 13.5 M_⊙ central object of Ilee et al. (2013) well reproduces the shape of the SED of M8EIR. A luminosity of ~12000 L_⊙ is sufficient to fit the multiple datasets, which is also in good agreement with the value of L obtained for the source in Ilee et al. (2013). In the preferred model a 5 × 10⁵au envelope, infalling at a rate of 1 × 10⁻³ M_⊙ yr⁻¹, surrounds this protostar. The presence of an envelope agrees with previous work by Simon et al. (1984) who interpreted M8EIR as two physically different components, a hot compact one and a broader cooler one, and stated that the broad component is not molecular cloud material, but an ‘intermediary phase’ of material between molecular cloud and hot disk material. A disk is present with a maximum radius of 2000 au and outflow cavities are carved out of the envelope with opening angles of 25° and densities of order 10⁻²⁰ g cm⁻³. Using a cavity density exponent of 0.25 was found to improve both the simulated 24.5 μm profile and the shape of the SED, similar to NGC 2264 IRS1 and IRAS 17216-3801. Using an inner hole (radius 30 au) again improved the fit of the simulated visibilities. Using the dust sublimation radius (~15 au) as the minimum radius led to the simulated visibilities for all configurations to be lower than those observed. In configuration (i) a sharp downturn is found at 8.5 μm. Such a quick visibility change can be caused by the turnover induced by the inner rim of the disk. When varying the inner disk radius around 30 au similar shapes were observed, but alongside further detriments to the other visibility datasets. As a result 30 au was chosen as the final inner radius. The simulated visibilities between 9 and 12 μm for configuration (ii) are ~0.02 higher than those observed, implying that the structure traced by this configuration is slightly larger than that modelled. Given the long baseline length of the configuration this could be explained if the inner rim is not completely smooth but clumpy in nature. If this inner radius were the dust destruction radius a sharp cut off would be expected, but given that an inner hole exists the idea that the medium could also not be perfectly smooth is not unlikely. A compromise was found between the fitting of configurations (iii) and (iv) as opposed to fitting one perfectly and the other very poorly. The visibilities are shown in Fig. A.17, the images in Fig. A.17 and the SED in Fig. A.19.

Fig. A.15 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux. The COMICS image is subject to a large amount of low level contamination, which can be seen in the 5% contours. In order to avoid this emission, only the inner regions of the source were considered in the fit.

Fig. A.16 Model SED of the best-fitting model of Mon R2 IRS2 (black). Multi-wavelength flux measurements from the RMS are represented as blue diamonds, the yellow diamond represents the VISIR flux density and the fluxes corresponding to the MIDIvisibilities are also shown in red. The unfilled diamonds represent the fluxes that were not considered in the fitting due to their suspected contamination and the fluxes from Henning et al. (1992) are included as squares.

A.7 AFGL 2136

AFGL 2136 IRS1 (hereafter AFGL 2136) was studied by de Wit et al. (2011) using a similar analysis to our own. The main difference was the amount of MIDI data available as only one configuration was fit in that work, as opposed to the four used here. Additional technical differences between our fitting and de Wit et al. (2011) are the multiplication of the model images with a Gaussian and the convolution of the COMICS image with a Gaussian PSF as opposed to an observed one (see Paper I). The COMICS image was presented in de Wit et al. (2011) but not included in the fitting process. The final model of de Wit et al. (2011) was used as a starting point for the modelling in our work. It was found that without multiplying by the Gaussian, a good fit could not be found to the MIDI configuration common to both this paper and de Wit et al. (2011). The model also provided poor fits to all the other MIDI configurations that were not included in de Wit et al. (2011).

Fig. A.17 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for M8EIR.

Maud et al. (2018) find that the dust emission implies a disk-like structure, so various disk-outflow-envelope geometries were trialled during the modelling. disk masses ranging between 10–30 solar masses are workable components of the model, in agreement with the kinematic measurements of Maud et al. (2018, 2019). The larger the mass, the less the 24.5 μm emission, improving the fit of the radial profile. For a 30 solar mass disk the millimetre slope of the SED fit is worse than for a 10 solar mass disk. The quality of the MIDI fits are comparable across the disk mass range. For the 30 solar disk mass model, the fit of configuration (ii) improves slightly as the simulated visibility at 8μm increases, but the visibilities for configuration (iv) also increase, which further worsens the fit to that configuration by a comparable amount. Given the similarity across the MIDI data a compromise was made between the SED and COMICS fits and 20 solar masses was selected as the final disk mass.

Similarity is found between the preferred best fitting model of our work and the model of de Wit et al. (2011) in that the minimum dust radius is larger than that of the dust sublimation radius. In this work the inner radius is four times larger than the sublimation radius, in de Wit et al. (2011) it was five. In de Wit et al. (2011) a gaseous disk and bloated star resided within this cleared inner hole. Including a gaseous disk is not part of the methodology of this work, but it could exist within the inner dust rim. As discussed in the examination of G305.20+0.21 in Paper I, the presence of a bloated source can also not be ruled out as our observables are sensitive to overall luminosity of the central source, not its specific radius or temperature. The final model includes an inclination for the system of 60°, and the central object is similar to that of de Wit et al. (2011); a bloated star of radius 25 R_⊙ with an overall luminosity of 1.51 × 10⁵ L_⊙. This luminosity matches the O-type classification ascribed by Maud et al. (2018, 2019) and in good agreement with the bolometric luminosity determined by Lumsden et al. (2013). Both a large and small grains disk were used in the final model, with matching scale-heights. The outflow cavity opening angle is tightly constrained at 22.5° to maintain the balance of fits of the three observables and the cavity density is of order 10⁻¹⁹ g cm⁻³. The centrifugal radius is 2000 au and the infall rate is of order 10⁻⁴ solar masses per year.

The visibilities of the final model are shown in Fig. A.20, the images are shown in Fig. A.20 and the SED is shown Fig. A.22. Configuration (ii) of the MIDI data proved nigh-impossible to well fit. Despite this, the preferred model recreated the shape of the visibilities well and was only V ² = 0.015 off in visibility. This implies that the protostellar environment traced by this configuration is larger than the model environment. The PA of this configuration was 14° meaning that, taking into account the source’s PA on the sky, it could be tracing envelope material close to the cavity wall or disk material. Both the continuum and CO contours of Maud et al. (2018) show the medium to be clumpy, despite the disk’s stability. Given that only smooth geometries were used in the modelling, it could be that the model is not reproducing this emission leading to this poor fitting. The other three MIDI configurations are matched well by the final geometry.

Fig. A.18 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.

Fig. A.19 Model SED of the best-fitting model (black) for M8EIR. All fluxes from the literature are shown as diamonds, except for the spectrum from Feldt et al. (2008). This was sampled electronically and the data are represented as purple crosses. Thefluxes corresponding to the MIDI visibilities are shown in red. Open symbols correspond to upper limits.

Fig. A.20 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for AFGL 2136.

Fig. A.21 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The model is more symmetric than the observed source meaning it has more emission in its eastern regions, explaining the larger amount of 24.5 μm flux visible for the model in the radial profiles. The contours in the images represent 5, 10, 25 and 75% of the peak flux.

Fig. A.22 Model SED of the best-fitting model (black). Multi-wavelength flux measurements from the RMS are represented as blue diamonds, the yellow diamond represents the COMICS flux density and the fluxes corresponding to the MIDI visibilities are also shown in red. The unfilled diamonds represent the fluxes that were not considered in the fitting due to their suspected contamination.

A further study at longer wavelengths, Maud et al. (2019), detect a ring-like structure at scales of 120 au. This is in very good agreement with the inner dust radius found in our model, which is four times the sublimation radius or ~125 au.

Appendix B List of MIDI observations

Appendix C List of SED fluxes

References

All Tables

All Figures

Fig. 1 Schematics of the suspected geometries of each of the sources with the position angles of each configuration of their MIDI data additionally shown. G305.20+0.21, the subject of Paper I, is also shown. The envelope is shown in blue, the cavity in yellow and the disk in red. Two different schematics are shown, one that applies to sources that have inclinations of less than 30 degrees. In this case, most of the cavity emission will be along the same region of the disk in the line of sight. The other schematic represents sources inclined to the point where cavity emission extends beyond the disk in the line of sight.

In the text

	Fig. 2 Relationship between the scale-height of the large-grains disk component (D1 SH) and the cavity opening angle.
In the text

	Fig. 3 Relationship between envelope infall rate and cavity opening angle.
In the text

	Fig. 4 Relationship between the luminosity of the source and the minimum disk radius.
In the text

	Fig. 5 Relation between luminosity and sublimation radius (orange dots). For comparison, the blue dots denote the disk minimum radii as in Fig. 4.
In the text

	Fig. A.1 Observed visibilities for W33A for each configuration (black) with the simulated visibilities (coloured) for different models to show the effects different components on the simulated MIDI data.
In the text

	Fig. A.2 Convolved model image (top, contours represent 5, 10, 25 and 75% of the peak flux) and the radial profile (bottom) of the best fitting model for the 2.4 kpc distance. The original COMICS image of W33A can be found in Wheelwright et al. (2012).
In the text

	Fig. A.3 SED of the preferred final model for W33A. Data from the RMS survey are shown as diamonds, the MIDI spectrum is shown as red crosses and data points from de Wit et al. (2010) are shown as triangles.
In the text

	Fig. A.4 Observed visibilities of NGC 2264 IRS1 for each MIDI configuration (black) with the simulated visibilities for each model image (coloured).
In the text

	Fig. A.5 COMICS 24.5μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.
In the text

	Fig. A.6 Model SED of the best-fitting model (black). Multi-wavelength flux measurements from the RMS are represented as blue diamonds, the yellow diamond represents the COMICS flux density and the fluxes corresponding to the MIDI visibilities are also shown in red.
In the text

	Fig. A.7 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured).
In the text

	Fig. A.8 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.
In the text

Fig. A.9 Model SED of the best-fitting model (black) for S255 IRS3. The different symbols correspond to different datasets, purple crosses are data from Caratti o Garatti et al. (2017), blue diamonds are fluxes from the RMS database, squares are data from de Wit et al. (2009) (with the COMICS flux in yellow) and the fluxes corresponding to the MIDI visibilities are also shown in red. Open symbols correspond to upper limits.

In the text

	Fig. A.10 Observed visibilities of IRAS 17216-3801 for each configuration (black) with the simulated visibilities for each model image (coloured) (1).
In the text

	Fig. A.11 Observed visibilities of IRAS 17216-3801 for each configuration (black) with the simulated visibilities for each model image (coloured) (2).
In the text

	Fig. A.12 VISIR 19.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.
In the text

	Fig. A.13 Model SED of the best-fitting model (black) of IRAS 17216-3801. Multi-wavelength flux measurements from the literature are represented as blue diamonds, the yellow diamond represents the VISIR flux density and the fluxes corresponding to the MIDI visibilities are also shown in red.
In the text

	Fig. A.14 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for Mon R2 IRS2.
In the text

Fig. A.15 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux. The COMICS image is subject to a large amount of low level contamination, which can be seen in the 5% contours. In order to avoid this emission, only the inner regions of the source were considered in the fit.

In the text

Fig. A.16 Model SED of the best-fitting model of Mon R2 IRS2 (black). Multi-wavelength flux measurements from the RMS are represented as blue diamonds, the yellow diamond represents the VISIR flux density and the fluxes corresponding to the MIDIvisibilities are also shown in red. The unfilled diamonds represent the fluxes that were not considered in the fitting due to their suspected contamination and the fluxes from Henning et al. (1992) are included as squares.

In the text

	Fig. A.17 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for M8EIR.
In the text

	Fig. A.18 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The contours in the images represent 5, 10, 25 and 75% of the peak flux.
In the text

	Fig. A.19 Model SED of the best-fitting model (black) for M8EIR. All fluxes from the literature are shown as diamonds, except for the spectrum from Feldt et al. (2008). This was sampled electronically and the data are represented as purple crosses. Thefluxes corresponding to the MIDI visibilities are shown in red. Open symbols correspond to upper limits.
In the text

	Fig. A.20 Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for AFGL 2136.
In the text

Fig. A.21 COMICS 24.5 μm image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image was convolved with the PSF of the observed object to accurately mimic the effects of the telescope specific to the observations. The model is more symmetric than the observed source meaning it has more emission in its eastern regions, explaining the larger amount of 24.5 μm flux visible for the model in the radial profiles. The contours in the images represent 5, 10, 25 and 75% of the peak flux.

In the text

	Fig. A.22 Model SED of the best-fitting model (black). Multi-wavelength flux measurements from the RMS are represented as blue diamonds, the yellow diamond represents the COMICS flux density and the fluxes corresponding to the MIDI visibilities are also shown in red. The unfilled diamonds represent the fluxes that were not considered in the fitting due to their suspected contamination.
In the text