Unveiling the traits of massive young stellar objects through a multi-scale survey
AAstronomy & Astrophysics manuscript no. first © ESO 2021February 11, 2021
Unveiling the traits of massive young stellar objects through amulti-scale survey * A. J. Frost , R. D. Oudmaijer , W. J. de Wit and S. L. Lumsden Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, Leuven, 3001, Belgiume-mail: [email protected] School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago, ChileReceived 23rd October 2020 / Accepted 8th February 2021
ABSTRACT
Context.
The rarity and deeply embedded nature of young massive stars has limited the understanding of the formation of stars withmasses larger than 8M (cid:12) . Previous work has shown that complementing spectral energy distributions with interferometric and imagingdata can probe the circumstellar environments of massive young stellar objects (MYSOs) well. However, complex studies of singleobjects often use di ff erent approaches in their analysis. Therefore the results of these studies cannot be directly compared. Aims.
This work aims to obtain the physical characteristics of a sample of MYSOs at ∼ ∼ Methods.
We apply the same multi-scale method and analysis to a sample of MYSOs. High-resolution interferometric data(MIDI / VLTI), near-di ff raction-limited imaging data (VISIR / VLT, COMICS / Subaru), and a multi-wavelength spectral energy distri-bution are combined. By fitting simulated observables derived from 2.5D radiative transfer models of disk-outflow-envelope systemsto our observations, the properties of the MYSOs are constrained.
Results.
We find that the observables of all the MYSOs can be reproduced by models with disk-outflow-envelope geometries, anal-ogous to the Class I geometry associated with low-mass protostars. The characteristics of the envelopes and the cavities within themare very similar across our sample. On the other hand, the disks seem to di ff er between the objects, in particular with regards to whatwe interpret as evidence of complex structures and inner holes. Conclusions.
The MYSOs of this sample have similar large-scale geometries, but variance is observed among their disk properties.This is comparable to the morphologies observed for low-mass young stellar objects. A strong correlation is found between theluminosity of the central MYSO and the size of the transition disk-like inner hole for the MYSOs, implying that photoevaporation orthe presence of binary companions may be the cause.
Key words. stars: formation – stars:imaging – stars:early-type – – techniques:interferometric – IR:stars
1. Introduction
Massive stars (those with masses greater than or equal to 8M (cid:12) )shape the Universe at all scales, with their winds, outflows,and supernovae (SNe) a ff ecting phenomena such as interstellarchemistry, the nature of molecular clouds, and the morphologyof galactic superwinds (Leitherer 1994). However, despite theimportance of massive stars, their formation is still not com-pletely understood. While forming low-mass stars lose their na-tal material (envelope) faster than their disk, revealing a clearstar-disk system, young massive stars remain deeply embed-ded. This, combined with the fact that they are numericallyrare, makes them more di ffi cult to investigate. One crucial stageof massive star formation is the massive young stellar object(MYSO). Here, the forming star is deeply embedded in its en-velope, rendering it optically invisible. Energy does escape anMYSO at micron wavelengths and longer with the most lumi- * Based on observations made with ESO telescopes at the ParanalObservatory under programme IDs 097.C-0320(A), 75.C-0755, 74.C-0389, 082.C-899(A), 076.C-0725(B), 84.C-0183(D), 84.C-0183(C),84.C-0183(B), 84.C-0183(B), 90.C-0717(A), 89.C-0968(A), 82.C-0899(A), 84.C-1072(A), 84.C-1072(B), 60.A-9224(A), 273.C-5044(A),74.C-0389(B), 75.C-0755(A), 75.C-0755(B), 75.C-0755(A), 75.C-0755(B), 77.C-0440(A) and 381.C-0607(A). nosity emitted at around 100 µ m, meaning they are IR bright.Despite the high bolometric luminosity of MYSOs, there is noionising radiation. This implies that the star is still accreting ata high rate, making the MYSO a key stage of the massive starformation process. With current instrumentation available forthe IR and (sub-)millimetre wavelengths, similar spatial reso-lutions ( ∼ − M (cid:12) yr − . This implies that the detecteddisks and outflows observed around MYSOs are engaged in theaccretion process, as they are in low-mass star formation. Thelow-mass star formation process has been assigned four distinctevolutionary stages based on the evolution of the envelope anddisk (Shu & Adams 1987). However, this has not been achievedfor massive star formation as massive protostellar disk evolu-tion in particular is yet to be constrained. Now that they have Article number, page 1 of 34 a r X i v : . [ a s t r o - ph . S R ] F e b & A proofs: manuscript no. first been confirmed to exist, determining the nature of massive pro-tostellar disks and their role in the accretion process for massiveprotostars 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 & Interferometer(VLTI). For example, de Wit et al. (2010) determined the im-portance of the bipolar outflow cavities as a source of mid-IRemission through the study of the MYSO W33A. Alternatively,Wheelwright et al. (2012) and de Wit et al. (2009) used Q-bandimaging to investigate the envelopes, allowing them to probecooler regions of the protostellar environments. These studiesmade use of radiative transfer modelling in order to comparesynthetic images to observed ones. These studies also consid-ered spectral energy distributions (SEDs), which provide a multi-wavelength view of MYSOs but only provide, at best, indirectspatial information. As such, they must be combined with spa-tially resolved observations to provide geometrical informationon the MYSOs.Another potential source of IR photons are the warmcircumstellar disks and the dust at sublimation temperature.In Frost et al. 2019 (hereafter Paper I), we combined theabove techniques to provide a multi-scale view of the MYSOG305.20 + ∼ µ m, Q-band direct imaging dataat ∼ µ m, and an SED with a 2.5D radiative transfer model,the characteristics of the MYSO were obtained. We found thatthe broad ∼ µ emission traces not only the cavity emission (deWit et al. 2010) but is also sensitive to the emission originatingin the disk. This is in agreement with previous work by Boleyet al. (2013), who found that their data from the MID-infraredInterferometer (MIDI, VLTI) can trace both disks and outflows.Furthermore, the inner dust radius, corresponding the inner dustdisk boundary, had to be over three times the sublimation radiusto satisfy the observables. The SED is dominated by emissioncoming from the outflow cavities and envelope, constituting re-processed emission from the star and hot inner disk regions. Thefact that the disk is likely to have an inner hole and has ratherlow envelope density suggests that this source is more evolved inits 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 addi-tional seven MYSOs. In total, this constitutes the largest MYSOsample to date to be subject to a multi-wavelength detailed anal-ysis of this sort. In doing so, we obtain a detailed set of character-istics for each MYSO that can be consistently compared to eachof the others, which has been lacking in the literature. In Section2, we describe the observations. In Section 3, the methodologyis briefly reiterated. In Section 4, the main results of the fittingprocess are described (an in-depth presentation of each source isprovided in Appendix A). In Section 5, we discuss the results ofthe sample. Finally, we summarise and conclude in Section 6.
2. Observations
The sample of MYSOs consists of eight sources located at a va-riety of distances, listed in Table 1. The sources were chosenbased on the availability of extensive observational data cover-age. The types of dataset used between this paper and Paper I arethe same, although the exact instrument and source of the datavaries slightly, the details of which we describe here.
An important source for the mid-IR interferometric data of thissample of MYSOs is Boley et al. (2013). These data were allcollected with the MIDI instrument. MIDI is a two-telescopebeam combiner that operates at N-band and delivers spectrallydispersed interferometric observables. It was fed by the VLT in-terferometer until its decommissioning in 2015. For all excepttwo sources, the MIDI data used were reduced and previouslypresented in Boley et al. (2013). In addition to the data from Bo-ley et al. (2013), one configuration of MIDI data from Grellmannet al. (2011) was used for one source (NGC 2264 IRS1) andthe MIDI data from de Wit et al. (2010) was used for W33A.The data of W33A were observed using a similar sequence tothe Boley et al. (2013) data, using the HIGHSENS mode (wherecorrelated and total flux measurements are taken separately andthe photometric observations immediately follow the interfero-metric measurements to produce the final visibilities) and thesame prism. They used calibrators from Cohen et al. (1999) aswell as a di ff erent version (1.6) of the data reduction softwareMIA + EWS, but the observing procedure followed was very sim-ilar to that of Boley et al. (2013). For the Grellmann et al. (2011)data an estimate of the sky background was made and subtractedover an area corresponding to the size of the mask over whichthe correlated spectra were extracted. The data from Grellmannet al. (2011) were taken in the SCIPHOT mode using a grismwith a spectral resolution of λ /∆ λ ≈ ff erence being the specific u-v coverageof each set of observations. This depends on the specific observ-ing configuration, and the areas probed for each source are illus-trated in Figure 1. A full list of the MIDI observations used inour work can be found in Appendix B. This paper uses single-dish Q-band data ( ∼ µ n) obtained withthe VISIR instrument at the VLT (Lagage et al. 2004) andthe COMICS instrument at the Subaru telescope. VISIR is amid-IR imager and spectrograph (mounted on Unit Telescope3 (UT3) of the VLT) and was used in service mode to observethe sources G305 and IRAS 17216-3801 as a part of ESO run097.C-0320(A) (PI: A. J. Frost). The observations were takenwith the Q-band filter which has a central wavelength of 19.5 µ mand a half-band-width of 0.4 µ m, tracing cooler material than thatof MIDI and therefore providing an alternative view of the pro-tostellar environments. The data reduction and further details onthis VISIR data were referenced in Paper I. COMICS observa-tions were used at a central wavelength of 24. µ m for the remain-ing sources, somewhat longer than the VISIR observations butstill sensitive to cavity wall and envelope emission. Images weretaken with a pixel size of 0.13 × and a field of viewof approximately 40 ×
30 arcsec . The COMICS data presentedhere were reduced as part of de Wit et al. (2009) and processedwith IRAF (Tody 1986).All Q-band observations were accompanied by PSF stan-dards and all target objects are spatially resolved. HD 123139was observed as PSF standard on the night of the VISIR obser-vations with a measured FWHM of 0.48" and α Tau was used asthe PSF during the post-processing of the COMICS observationsand has a FWHM of 0.34". We performed aperture photometryon the Q-band images in a standard way comparing the flux be-haviour of flux standard throughout the night. In our modelling,
Article number, page 2 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. 1.
Schematics of the suspected geometries of each of the sources with the position angles of each configuration of their MIDI data additionallyshown. G305.20 + ff erentschematics are shown, one that applies to sources that have inclinations of less than 30 degrees. In this case, most of the cavity emission will bealong the same region of the disk in the line of sight. The other schematic represents sources inclined to the point where cavity emission extendsbeyond the disk in the line of sight. we reduce the 2D images to azimuthal averaged radial profiles,which are then compared to synthetic data. The fluxes for the SEDs were compiled on an object-by-objectbasis. The Red MSX Source (RMS) Survey (Lumsden et al.2013) database was used in addition to the literature to create anobservational SED for each object. The literature was consultedto ascertain whether, as was the case for G305, other sources could be contaminating various fluxes, and if this was the casesuch measurements were not considered in the fitting. All theincluded SED fluxes are listed in Appendix C.
3. Methodology
The methods of this work closely follow those of Paper I. Inshort, the multi-wavelength dataset described in Section 2 isused to probe the protostellar environments of MYSOs at 10masand 100mas scales and as a whole. By fitting one RT model toall these observations simultaneously, the characteristics of the
Article number, page 3 of 34 & A proofs: manuscript no. first
Table 1.
List of the MYSOs included in this study with a column stating whether VISIR or COMICS data were used in the fitting process.
References: (1) Lumsden et al. (2013); (2) Mottram et al. (2011a); (3) (Mottram et al. 2011b); (4) Immer et al. (2013); (5) Grellmann et al. (2011);(6) Kamezaki et al. (2014); (7) (Wang et al. 2011); (8) Burns et al. (2016); (9) Henning et al. (1992); (10) Herbst & Racine (1976); (11) Krauset al. (2017); (12) Boley et al. (2013); (13) Linz et al. (2009); (14) Damiani et al. (2019); (15) Urquhart et al. (2012); (16) Urquhart et al. (2014).
Name RA (J2000) Dec. (J2000) log( LL (cid:12) ) Distance VISIR (V) or COMICS (C)?(h:m:s) (d:m:s) (kpc)G305.20 + VW33A 18:14:39.0 -17:52:03 4.5 CNGC 2264 IRS1 06:41:10.15 + CS255 IRS3 06:12:54.02 17:59:23.60 4.7 CMon R2 IRS2 06:07:45.8 -06:22:53.2 3.8 CIRAS 17216-3801 17:25:06.51 -38:04:00.4 4.8 VM8EIR 18:04:53.18 -24:26:41.4 3.8 CAFGL 2136 18:22:26.38 -13:30:12.0 5 CMYSO environment at these multiple scales are probed, provid-ing a multi-scale view of the MYSOs.HOCHUNK-3D (Whitney et al. 2013) is the code used togenerate the RT models from which the required simulated ob-servables are extracted. The model consists of two main geo-metrical components; an envelope with bipolar outflow cavitiesand a disk, which surround a central object. Polynomial-shapedor streamlined cavities present the two main geometrical optionsfor the bipolar outflow cavities, and two envelope density dis-tributions - power-law and Ulrich (Ulrich 1976) (a rotating andinfalling envelope) - are explored. The disk is described by twoseparate density components. One is a large-grain disk compo-nent and the other a small grains disk component, which are as-signed di ff erent dust types. The base stellar parameters sampledare the stellar radius, mass, and temperature. We find that it isthe stellar luminosity, as opposed to the individual radius andtemperature, which a ff ect our observables, so this is discussedin later sections. The disk is modelled via two separate diskcomponents, a large-grains disk (disk 1 - described by Model1 from Wood et al. 2002) and a small grains disk (disk 2 - de-scribed by the interstellar medium (ISM) grain model of Kimet al. 1994). For each component the sample parameters werethe inner radius, outer radius, and scale-height. In Paper I, onlythe large-grain disk component was included; in this work, bothdisk components are used. This allows for di ff erent geometricalarrangements of grains in the disk to be experimented with. Thedisk scale-height and the inner and outer radii are free param-eters for each of the disk components. Additionally, we fit forthe total disk mass and the relative contribution by each compo-nent. The scale-height exponent and the radial density exponentare fixed at values 2.2 and 1.2 respectively. These values are rec-ommended by Whitney et al. (2013) based on the hydrostaticcalculations of D’Alessio et al. (1998). The overall disk parame-ters were the total mass of the disk and the ratio of this disk massbetween the two disk components. For one source, the spiral fea-tures of the code were also required, where the free parameterswere the pitch angle of the spiral, the summation parameter (thatdetermines the width of the arms), the fraction of mass entrainedin the arms and the radius where the spiral arms begin. The cav-ity parameters investigated were the cavity opening angles, thecavity density exponent, and the coe ffi cient for cavity densitydistribution. The envelope free parameters were the centrifugalradius, the minimum and maximum envelope radii, and the en-velope infall rate.The dust types within the models are not varied through-out the fitting process. The dust type for the large-grains diskfollows a dust distribution based on models from Wood et al. (2002), with grain sizes ranging between 5 µ m and 1mm. Thecavity and ‘small grains’ disk share the same dust type, basedon that of Kim et al. (1994) for the di ff use ISM (with maxi-mum grain sizes of 0.01 µ m). For the envelope grains, Whitneyet al. (2003) produced a size distribution of grains in agreementwith R v ∼
4, which is typical of dense regions of molecularclouds (Whittet et al. 2001). To include the e ff ects of ices inthe cooler outer regions of the envelope a layer of water icewas included as the outer 5% of the grains’ surface. This finalpercentage was decided upon after variation testing and com-parison to polarisation observations, as described in Whitneyet al. (2003). Polyaromatic-hydrocarbons (PAHs) and very smallgrains (VSGs) can also be accounted for in the code, resultingin emission that is visible in the SED and also has e ff ects on theother observables. If PAH emission is recorded for a source inthe literature, the dust files which are identical to the above filesare used, with the exception that they do not include the grainpopulation with sizes smaller than 200Å. These are instead cal-culated for using the dust distribution of Draine & Li (2007).PAHs were not detected for the source studied in Paper I, so theresults of including these PAHs on our simulated observables arepublished here for the first time.The interferometric observables are simulated in Fourierspace. HOCHUNK-3D model images are generated at 1 µ m in-tervals within the MIDI N-band range, the fast-Fourier trans-form is taken of the images, and the specific u-v point or visibil-ity according to the position angle and baseline of the observa-tions is extracted, corresponding to the data point obtained withMIDI. By applying this to each image, a visibility profile withwavelength is obtained. In order to simulate the Q-band images,model images are generated at the required wavelength and con-volved to the telescope’s resolution. Azimuthally averaged radi-ally profiles are then used to compare the model, source and PSFimages. The SEDs were generated with IDL post-processingscripts (Whitney et al. 2013), which could then be compared withphotometric data points compiled from the literature. The fitting process is the same as that of Paper I, with the SEDconstituting the starting point for the fitting, before the high-resolution datasets were considered. The Ulrich-type envelopeand cavities successfully used in de Wit et al. (2010) and Wheel-wright et al. (2012) were the starting geometry for the envelopeand cavity specifications in this work. The fitting process startswith the simplest geometry, which consists of only of an enve-
Article number, page 4 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs lope. If a fit could not be obtained, then the model geometry be-comes progressively more complex by changing outflow cavityparameters and adding disks and eventually disk substructure.When combining a variety of observations, the balance be-tween the fits of those observations needs to be considered.When the SED alone is utilised, a common method of uncer-tainty quantification has been the chi-square minimisation. TheSED Fitter tool of Robitaille et al. (2007), which utilises a largegrid of models and has been used to analyse the data of manyprotostars, employs this minimisation. Once a method is ex-panded past SED considerations, which as previously discussedis important to better constrain protostellar environments, the un-certainties in the fit parameters must be and have been consid-ered di ff erently. This is because di ff erent observations, at di ff er-ent wavelengths and using di ff erent techniques will all trace dif-ferent scales of MYSO environments in contrasting ways. Fewpapers dedicated to MYSOs in the literature fit RT models tomore than one set of observations, but the papers that are mostsimilar to our approach have tackled this is di ff erent ways.Linz et al. (2009) fit MIDI interferometric data and an SEDfor one MYSO. They used the aforementioned SED Fitter toolto constrain the fits on their SED and compared the model visi-bilities associated with a variety of acceptable models from thatprocess in order to decide upon their best fitting model. A similarapproach was taken in de Wit et al. (2010), where the SED Fit-ter was again used to provide some quantitative analysis of theSED fits, but the overall fit between the MIDI and SED was de-termined by eye. In de Wit et al. (2011), where one MIDI datasetand an SED were used, the overall fit was determined by eye.Johnston et al. (2013) modelled the SED and 2MASS imageprofiles of the MYSO AFGL 2591 using a Monte Carlo radia-tive transfer code. The reduced combined chi-square was calcu-lated for the line profiles and their SEDs, giving equal weightto each observation within the final chi-square. They studied theconstraints on their parameter space by systematically varying afew key parameters and investigating how the chi-squared variedwith each parameter.Boley et al. (2013) determined their best fits among 2DGaussian models fit to their MIDI data by performing a gridsearch over their parameter space and the chi-square minimisa-tion. The uncertainties on the fitting parameters were found us-ing a Monte Carlo approach. This involved generating syntheticobservations, which were normally distributed about the meanvisibility of their observations, and repeating this to optimise.This was a viable approach given that only one observable (theirMIDI data) was being used in their work. In Boley et al. (2016), avariety of di ff erent models were applied to a combined dataset ofMIDI observations, T-ReCS aperture masking data and K-bandinterferometric data for MYSO. Gaussian modelling was againused to fit their visibilities, and the chi-square minimisation (asin Boley et al. (2013)) was used to determine good fits; howevera di ff erent fit found for the di ff erent wavelength datasets and nooverall chi-square values were computed. Temperature-gradientdisk models were also used to fit the K-band visibilities and theN-band correlated fluxes and a reduced chi-squared was calcu-lated across the di ff erent observations; however, the weightingwas not explicitly explained. Finally, radiative transfer modelswere employed via the use of the 2003 version of HO-CHUNK(Whitney et al. 2003). They did not attempt to reproduce theSED throughout this modelling; only the K-band visibilities andN-band fluxes were modelled. No errors are formalised on theirfinal parameters, and the chi-square was not utilised during thisportion 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. Followingthe method of De Buizer et al. (2017), they performed the chi-square minimisation to find five best fitting models for each ofthe MYSOs. Five parameters alone were varied within the fittingprocess in order to investigate how the SED is a ff ected by modelscorresponding to di ff erent evolutionary phases that occur withinthe turbulent core model (Tan et al. 2014).The aforementioned works constitute those most similar tothe work on this sample, but more observations are used hereand all are fit simultaneously. Hardly any of these works cal-culated a global chi-square for di ff erent datasets, nor did theyderive formal errorbars on their parameters when RT codes wereused. This appears to be the inherent compromise of combiningmany datasets and using a complex modelling approach. Includ-ing multiple observations and a large model parameter space pro-vides a much better picture of a source, but formalising the com-bined error using these methods becomes arduous. Calculatingan overall chi-square and giving each observation equal weight,as has been done in some of the works mentioned in the previoussubsection, was not determined to be an appropriate approach forthe work on our sample. The MIDI spectrum and VISIR profileshave dozens of points covering a very small wavelength regime,while the SED covers a wavelength regime from one to thou-sands of micron with a (most frequently) small number of datapoints. This means that there will be an obvious bias within anychi-squared value calculated for the combined datasets, in favourof the 7-13 µ m and ∼ µ m datasets. Another reason to not com-bine the chi-squared is that each observation traces very di ff erentscales so, similarly, the innermost and smallest geometrical fea-tures of the source, which are traced by the data-rich MIDI spec-trum, will be unfairly focused on through the use of a combinedchi-square.Overall, it can be concluded that while the chi-squared maynot be an appropriate way to provide a blanket estimate of a fitwhen multiple observations that trace di ff erent scales are used, itcan assist in the determination of which observables are key fortracing di ff erent physical characteristics of MYSOs. In order tocomment on the uncertainties involved in the fitting of the wideparameter space provided by the RT code, the chi-squared wascalculated for each observation. The chi-squared was calculatedseparately for each configuration of the MIDI observations andthen averaged. The image profile chi-squares were calculatedfrom the radial profiles up to the point where the backgroundnoise exceeded the rms. The SED chi-squared was calculated byconvolving the model SED with filter fluxes corresponding to theavailable data points. While the SED was the starting point forthe fitting (as in de Wit et al. (2010); Linz et al. (2009)) the SEDfitter tool is not used. As such, the fits of these models avoidthe inherent biases of the Robitaille grid of models, which arediscussed at length in Robitaille (2017).Following the above considerations and investigations intothe individual chi-squares, we find that the di ff erent parameterscan be traced within certain ranges. The disk inner radius couldbe tightly constrained within 5 ◦ , as could the cavity opening an-gles. Generally, significant changes were found in the MIDI fitsfor variations of the order 0.1 in scale-height in the two disk com-ponents. Similarly the cavity density exponent was constrainedwithin ± Article number, page 5 of 34 & A proofs: manuscript no. first mass could be constrained to an order of magnitude, and the en-velope outer radius could be constrained within a factor of five.For sources that were particularly close or orientated such thatthe disk emission dominated or had a mass-dependent substruc-ture such as a spiral, the disk mass could be better constrained,and we discuss these ranges on a source by source basis in theAppendix.
4. Results
In this section, we present the results derived when applying theabove methods to the sample of MYSOs. Table 2 describes sum-marises the key features of each MYSO, while Table 3 lists theirspecific parameters. An in-depth discussion of the fitting and fea-tures of each source is provided in Appendix A.G305.20 + + ff ects on the infall properties, andtherefore the density distribution, of the envelopes.The envelope properties of the MYSOs were traced predom-inately by the longer wavelength emission of the SED and the20 µ m emission, which is sensitive to the centrifugal radius (forcloser sources) and the envelope infall rate. The larger the en-velope infall rate, the larger the density within the infall radius(according to the Ulrich density distribution Ulrich 1976) andthe larger the amount of material available for re-emission of the20 µ m flux. The infall rate reaches values of ∼ − M (cid:12) yr − for oursample. This is consistent with infall rates reported by McKee& Tan (2003) for the turbulent core model. Accretion rates ofthis magnitude are reported in numerical work (e.g. Hosokawaet al. 2010, Kuiper et al. 2011), which is consistent assumingthat envelope infall rates translate to dynamical mass accretionrates onto the star. The envelope outer radii are all similar acrossthe sample, but our high-resolution data are not sensitive to suchlarge scales. All the MYSOs studied within our work have ob-servables that can be fit with an Ulrich-type envelope, imply-ing that all of the MYSOs are actively infalling and contain thesubsequent density enhancements in the midplane that a ff ect the20 µ m emission. Follow-up kinematic studies could confirm thisto be the case following methods such as those detailed in Ohashi(2004) (e.g. searching for inverse P Cygni line profiles in thesources’ millimetre spectra). S255 IRS3 appears to be the ex-ception in terms of envelope parameters, given its low envelopeinfall rate and smaller outer radius (by a factor of five). S255 isthe only source confirmed to be episodically accreting so the factthat this source’s envelope infall rate is lower is particularly in-teresting, as it could be related to the accretion activity and thefact that the observations in this paper correspond to an accretionlull. All the MYSOs in this sample have disks as a component oftheir final models. These disks are not toroids, which are 100s ofsolar masses (Beltrán & de Wit 2016), but are objects that (whilethey could exist at larger radii) are present in the innermost( ∼ ff ects on the MIDI visibilities. For the majority of the modeldisks, the small grains disk component has a scale-height twicethat of the large-grains disk. Such a ratio is in agreement with thehydrostatic calculations of D’Alessio et al. (1998) and models ofobserved images (Cotera et al. (2001) and Wood et al. (2002)).IRAS 17216-3801 has a flaring factor very close to the valuemost commonly observed throughout the sample of 2 at 1.75,but the fact that a binary system is likely to be having dynamicale ff ects on this disk could explain this. Two of the remaining disksare flat (scale-heights of the disks were 0.1) and one of themrequired no small grains to be present in the disk to successfullyfit the datasets.Whitney et al. (2013) recommended that the distribution ofthe dust be split such that 20% of the grains are in disk 1 (thelarge-grains disk component) and the remainder are in disk 2 (thesmall grain component). This was based on the expected dustsettling timescales observed for low-mass source. Even thoughthe exact nature of dust amalgamation within massive protostel-lar environments is yet to be accurately constrained, we find thatthis approximation suited all of the observations well, save forthose of G305. The fact that no disk 2 was required to fit theMIDI visibilities for this source was discussed as a significant el-ement that characterises this object in a possible transition phasein its disk evolution towards full disk dispersion, as discussed inPaper I.Inner clearing, as found for the disk of G305.20 + ff erent form of substructure.The spiral parameters of the code were utilised, and we find thatincluding a spiral with a loose pitch angle of 30 ◦ , which con-tains 90% of the mass of the disk and begins and ends at 30 and60au, respectively, assisted in fitting the observables. These pa-rameters essentially create substructure that manifests as a diskwith a gap rather than an obvious spiral. NGC 2264 IRS1 is byfar the closest MYSO of the sample and, as such, was traced atthe smallest scales, so the fact that evidence is found for disksubstructure only in the nearest of all sources in our sample isdeemed consistent.Of the disk parameters, one weaker constrained parameter isthe disk mass. The only source where the disk mass was con-strained to less than a factor of ten was NGC 2264 IRS1. This islikely due to the fact that it is near, that it has a close to face-oninclination, and that it has a mass-dependent substructure. Theouter disk radius is only loosely constrained by the observables.A larger disk creates more millimetre emission, which can af-fect the longer wavelength regime of the SED and the maximumdisk radius is also tied to the centrifugal radius, which a ff ectsthe shape of the SED particularly at its 70 µ m peak and around20 µ m. The imaging profiles also proved decent outer disk radiitracers, but mostly for the nearer sources. As such, the outer ra-dius is still discussed but it is noted that this is not as constrainedas the inner radius, which heavily a ff ects the N-band emission inall cases. Article number, page 6 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Table 2.
Short descriptions of the geometry inferred from the parameters of the final models of each MYSO.
Name GeometryG305.20 + R sub ( R min = R envmin < R diskmin )M8EIR Flared disk with a cleared inner hole 20au in radiusAFGL 2136 125au radius inner hole
5. Discussion
We have performed RT modelling of a sample of eight MYSOsinvoking a star-disk-envelope configuration. In the followingsection, we discuss common features and traits identified in thesample. Detailed discussions of the individual sources can befound in Appendix A.
Two factors are also defined to assist in the interpretation ofthe parameter space, resolving factor and flaring factor. The re-solving factor was calculated to determine whether the resolvingpower of the MIDI observations a ff ected the quantities derived.This was calculated by dividing the largest projected baseline(Vinkovi´c & Jurki´c 2007) used to observe each source by the dis-tance to the source. The flaring factor is simply the scale-heightof the ‘small grains’ disk (disk 2) divided by the scale-height ofthe ‘large-grains’ disk (disk 1).In order to make general comparisons between the sources inour sample, we calculate Pearson correlation coe ffi cients. Corre-lations were inspected for basic parameters (bolometric luminos-ity, central object mass and distance) and the modelling param-eters (cavity opening angle, disk flaring factor, envelope infallrate, cavity density, centrifugal radius and minimum and maxi-mum disk radius). We do not calculate correlations for any pa-rameters that were determined to be poorly constrained, namelydisk mass and envelope radius.Before considering the correlations they were checked forbias. To check whether the specific MIDI observations a ff ectedthe correlations found, correlations between the resolving fac-tor and the parameter space were calculated. Strong correlationswere detected for two characteristics - the distance and stellarmass, which seem consistent, as a very distant source needs tobe more massive or luminous to be detectable, but if a source isclose and observed with small baseline lengths it may be com-parably resolved as a source that is distant but observed withthe largest baselines. We find that no significant correlations arefound between the inclination and the rest of the parameter spaceand that therefore the results found for the sample are not a resultof fortuitous orientations. We conclude that any strong correla-tions beyond these above are tangible results from the parameterspace. The vast majority of the correlation coe ffi cients are < Fig. 2.
Relationship between the scale-height of the large-grains diskcomponent (D1 SH) and the cavity opening angle.
A strong correlation is detected between the scale-height of disk1 and the cavity opening angle (0.72) implying that the moreflared the disk, the larger the cavity opening angle (Figure 2). Ifthe disk scale-height is large, the disk is more flared and a largeremitting surface will be visible. This same a ff ect occurs if thecavity opening angle increases as this in turn increases the avail-able area for irradiation by the protostar and therefore generatesa larger emitting region. The correlation between disk 2 scale-height and cavity opening angle is moderate (0.62) because onemodel does not possess this second disk component.A correlation of 0.75 is found between the envelope infallrate and cavity opening angle, implying that as the envelope in-fall rate increases, the cavity opening angle will get larger (Fig-ure 3). The Ulrich geometry the envelope model uses is based ona rotating and infalling solution and results in a build of materialat the equatorial plane. If the envelope infall rate is high, morematerial is moved from the outer envelope towards the midplaneit becomes more concentrated while the regions closer to thepoles of the envelope are less dense. The remaining envelopematerial would therefore be more vulnerable to the e ff ects ofany outflow activity from the central source, meaning that more Article number, page 7 of 34 & A proofs: manuscript no. first T a b l e . P a r a m e t e r s o f t h e p r e f e rr e d m od e l s f o r t h e f u ll s a m p l e o f s ou r ce s . P A i ss ho r t h a nd f o r po s iti on a ng l ea nd i s d e fi n e d a s ea s t fr o m no r t h . E x f . s t a nd s f o r ca v it yd e n s it y e xpon e n t ( W h it n e y e t a l . ) , D S H i ss ho r t f o r‘ d i s k1 s ca l e - h e i gh t ’ a nd D S H i ss ho r t f o r‘ d i s k2 s ca l e - h e i gh t ’ . S ou r ce M (cid:63) L (cid:63) i d R m i n e nv R m a x e nv R c ˙ M i n f a ll θ c a v n c a v E x f M d i s k R m i nd i s k R m a xd i s k D u s t F r ac ti on D S HD S H A f o r v P A ( M (cid:12) )( L (cid:12) )( ◦ )( kp c )( a u )( a u )( a u )( M (cid:12) y r − )( ◦ )( g c m − )( M (cid:12) )( a u )( a u )( ◦ ) W A . ( R s ub ) × . × − × − ( R s ub ) . . . G . + . × . × − × − . - NG C I R S . ( R s ub ) × × − × − . . ( R s ub ) . . . - . S I R S . ( R s ub ) × . × − × − ( R s ub ) . . . I R A S - . × × − × − . . . M on R I R S . × . × − × − . . . . M E I R . . × × − × − . . . . . A F G L . × × − . × − . . . - Fig. 3.
Relationship between envelope infall rate and cavity openingangle.
Fig. 4.
Relationship between the luminosity of the source and the mini-mum disk radius. material could be carved of the envelope, thereby widening theoutflow cavity such that it has a larger opening angle, explainingwhy this trend is observed.A strong correlation (0.92) exists between disk minimum ra-dius and luminosity (Figure 4). If considering straightforwarddust sublimation this is to be expected, yet the minimum dust ra-dius does not correspond to the actual dust sublimation for manysources. This is illustrated in Figure 5. In Paper I causes of innerholes in disks were discussed at length and it was determinedthat the most likely causes for inner holes in massive stars are
Article number, page 8 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. 5.
Relation between luminosity and sublimation radius (orangedots). For comparison, the blue dots denote the disk minimum radii asin Fig. 4.
1) photoevaporation and 2) the presence of binary / multiple com-panions. Of these two causes, photoevaporation as a mechanismseems to be in agreement with this trend. If a source is moreluminous, then it will have a greater capacity to generate pho-tons of high enough energies to dissociate the inner disk materialand create cleared inner regions. Therefore, the large luminosityleading to a large inner hole makes sense as the star has anothermechanism through which to destroy the disk. This comes withthe caveat, however, that the star must not also be creating UVphotons to the point where an ionising ultra-compact HII regionhas begun to form around the source, as that is not observed forany of these sources. A strong correlation does not exist betweenminimum envelope radius and luminosity, as the minimum enve-lope radius of the final model for IRAS 17216-3801 vastly di ff ersfrom its minimum disk radius. If IRAS 17216-3801 is removed,the correlation between minimum envelope radius and luminos-ity is becomes strong and the correlations between luminosityand both disk and envelope minimum radius become 0.93. In our sample of eight MYSOs, selected merely on the basisof completeness of observations, our analysis indicates that fivesources have evidence for substructure in their disks. A diskwith substructure in our terminology is contrasted with an az-imuthally symmetric disk extending inwards down to the dustsublimation radius. For all except one, this takes the form of acleared inner hole. The selection of objects based on complete-ness of IR spatial and SED data implies a number of things forour sample compared to massive stars in general. Firstly, theobjects are relatively close by and only G305.20 + / METIS or a much denser filling of theuv-plane with VLTI / MATISSE could be employed to investigatedust sublimation and disk evaporation and how this relates to thephysical presence of an inner disk rim.The remaining source, NGC 2264 IRS1, showed substruc-ture in the form of a gap-like spiral feature. The fact thatthese substructures are present has interesting implications forthe massive star formation process. To date, substructures haverarely been reported around massive protostars. Sanna et al.(2019) detect a sub-Keplerian disk and jet system around theMYSO G023.01-00.41 and find that the CH OH emission atscales ∼ (cid:12) object HD 97048, but otherwise the presence of planets inthese gaps has mostly been inferred. As discussed in the contextof G305 in Paper I, planets are not observed around massive starsbut whether they can never form or cease to exist as the protostarbecomes more luminous with age has not been confirmed. Analternative method of gap creation is dust trapping. Dust traps(Whipple 1972) occur when a local gas pressure maximum inthe disk attracts and traps dust grains allowing clearing to oc-cur in the disks and has also been proposed as a mechanism forobserved gaps 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 thedisk of the Herbig Ae star HD 163296 can be reproduced fromdust trapping alone, without the presence of a planetary body.Dust growth presents another avenue of gap formation wherebythe amalgamation or pile-up of dust grains clears space in thedisk (e.g. Zhang et al. (2015), Okuzumi et al. (2016)). Instabil-ities, both large-scale (Lorén-Aguilar & Bate 2016) and grav-itational (Takahashi & Inutsuka 2016), have also been shownthrough modelling to be able to cause gaps. Among these pro-cesses, the ones that would operate in a disk around a massivestar are probably grain growth and instabilities, provided thatthe young massive star accretes from a high density (and mas-sive) disk. Toroids, which have been shown to create instabili-ties in works such as Lorén-Aguilar & Bate (2016), have beenobserved around massive forming stars (such as those studiedin Beltrán et al. (2005)), presenting an avenue for this instabil-ity and therefore gap to occur. Gravitational instability is a morelikely occurrence in massive systems where the central proto-stellar mass is higher and when multiple systems are present.Price et al. (2018) have shown through hydrodynamical simula-tions that spirals, shadows and asymmetrical structures can all beproduced by a companion star on an inclined and eccentric or-bit approaching periastron. Both the e ff ects of binaries and dustgrowth take 10s of thousands of years to occur, and are thereforealso better suited to an MYSO in a later phase of formation. Article number, page 9 of 34 & A proofs: manuscript no. first
6. Conclusions
In this paper we have presented a multi-scale and multi-wavelength analysis of a sample of eight MYSOs with the aimto make direct comparisons of their characteristics. We derivephysical parameters from the observational data by means ofradiative transfer modelling, fitting simultaneously spatially re-solved and SED data. The RT model consists of a disk-outflow-envelope morphology and we find that such a geometry pro-vides reasonable fits to the data. Envelopes were best describedby Ulrich-type density distributions, suggesting that they couldbe dynamically unstable and implying that they are actively in-falling environments.Agreement exists between the characteristics found throughthis work and those found through previous studies. From thestudy of this sample it appears that the the ∼ µ m interfero-metric data (probing (cid:47) ∼ µ m single-dish data which provides additional constraint oncavity material. Variation is seen at the small scales ( (cid:47) / multiple system are the most likely candidates forany substructure in this disk. All the inferred substructures havebeen 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 diskevolution. Future work to detect multiplicity and further improvethe view of these inner regions through spectrointerferometric,imaging and image reconstruction methods, will assist in con-firming these substructures and distinguishing their origins. Acknowledgements.
We thank the anonymous referee for their helpful pointersand Prof. Melvin Hoare and Prof. Steven Longmore for their discussions, allof which enabled the improvement of this work. We also thank the STFC forfunding this PhD project. The modelling presented in this work was undertakenusing ARC3, part of the High Performance Computing facilities at the Universityof Leeds, UK.
References
Alvarez, C., Hoare, M., Glindemann, A., & Richichi, A. 2004, A&A, 427, 505Aspin, C. & Walther, D. M. 1990, A&A, 235, 387Beichman, C. A. 1988, Astrophysical Letters and Communications, 27, 67Beltrán, M. T., Brand, J., Cesaroni, R., et al. 2006, A&A, 447, 221Beltrán, M. T., Cesaroni, R., Neri, R., et al. 2005, A&A, 435, 901Beltrán, M. T. & de Wit, W. J. 2016, A&A Rev., 24, 6Boley, P. A., Kraus, S., de Wit, W.-J., et al. 2016, A&A, 586, A78Boley, P. A., Linz, H., van Boekel, R., et al. 2013, A&A, 558, A24Burns, R. A., Handa, T., Nagayama, T., Sunada, K., & Omodaka, T. 2016, MN-RAS, 460, 283Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nature Physics,13, 276Chini, R., Kreysa, E., Mezger, P. G., & Gemuend, H. P. 1986a, A&A, 154, L8 Chini, R., Kruegel, E., & Kreysa, E. 1986b, A&A, 167, 315Cohen, M., Walker, R. G., Carter, B., et al. 1999, AJ, 117, 1864Cotera, A. S., Whitney, B. A., Young, E., et al. 2001, ApJ, 556, 958Cutri, R. M. & et al. 2012, VizieR Online Data Catalog, II / / IRAC Star For-mation Reference SurveyFeldt, M., Pascucci, I., Chesneau, O., et al. 2008, in The Power of Optical / IR In-terferometry: Recent Scientific Results and 2nd Generation, ed. A. Richichi,F. Delplancke, F. Paresce, & A. Chelli, 263Frost, A. J., Oudmaijer, R. D., de Wit, W. J., & Lumsden, S. L. 2019, A&A, 625,A44Fujisawa, K., Yonekura, Y., Sugiyama, K., et al. 2015, The Astronomer’s Tele-gram, 8286Galván-Madrid, R., Zhang, Q., Keto, E., et al. 2010, ApJ, 725, 17Giannakopoulou, J., Mitchell, G. F., Hasegawa, T. I., Matthews, H. E., & Mail-lard, J.-P. 1997, ApJ, 487, 346Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004,ApJS, 151, 35Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984Grellmann, R., Ratzka, T., Kraus, S., et al. 2011, A&A, 532, A109Guertler, J., Henning, T., Kruegel, E., & Chini, R. 1991, A&A, 252, 801Harvey, P. M., Campbell, M. F., & Ho ff mann, W. F. 1977a, ApJ, 215, 151Harvey, P. M., Campbell, M. F., & Ho ff mann, W. F. 1977b, ApJ, 215, 151Henning, T., Chini, R., & Pfau, W. 1992, A&A, 263, 285Herbst, W. & Racine, R. 1976, AJ, 81, 840Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478Ilee, J. D., Cyganowski, C. J., Brogan, C. L., et al. 2018, ApJ, 869, L24Ilee, J. D., Wheelwright, H. E., Oudmaijer, R. D., et al. 2013, MNRAS, 429,2960Immer, K., Reid, M. J., Menten, K. M., Brunthaler, A., & Dame, T. M. 2013,A&A, 553, A117Ishihara, D., Onaka, T., Kataza, H., et al. 2010, VizieR Online Data Catalog,2297Itoh, Y., Tamura, M., Suto, H., et al. 2001, PASJ, 53, 495Jiménez-Serra, I., Báez-Rubio, A., Martín-Pintado, J., Zhang, Q., & Rivilla,V. M. 2020, ApJ, 897, L33Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19Johnston, K. G., Shepherd, D. S., Robitaille, T. P., & Wood, K. 2013, A&A, 551,A43Kamezaki, T., Imura, K., Omodaka, T., et al. 2014, ApJS, 211, 18Kastner, J. H., Weintraub, D. A., & Aspin, C. 1992, ApJ, 389, 357Kim, S.-H., Martin, P. G., & Hendry, P. D. 1994, ApJ, 422, 164Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339Kraus, S., Kluska, J., Kreplin, A., et al. 2017, ApJ, 835, L5Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20Lagage, P. O., Pel, J. W., Authier, M., et al. 2004, The Messenger, 117, 12Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599Leitherer, C. 1994, in Reviews in Modern Astronomy, Vol. 7, Reviews in ModernAstronomy, ed. G. Klare, 73–102Linz, H., Henning, T., Feldt, M., et al. 2009, A&A, 505, 655Liu, M., Tan, J. C., De Buizer, J. M., et al. 2019, arXiv e-prints, 901.01958, ApJ,submittedLongmore, S. N., Burton, M. G., Minier, V., & Walsh, A. J. 2006, MNRAS, 369,1196Lorén-Aguilar, P. & Bate, M. R. 2016, MNRAS, 457, L54Lumsden, S. L., Hoare, M. G., Urquhart, J. S., et al. 2013, ApJS, 208, 11Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2019, A&A, 627, L6Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2018, A&A, 620, A31Maud, L. T., Hoare, M. G., Galván-Madrid, R., et al. 2017, MNRAS, 467, L120McKee, C. F. & Tan, J. C. 2003, ApJ, 585, 850Merello, M., Evans, Neal J., I., Shirley, Y. L., et al. 2015, ApJS, 218, 1Mezger, P. G., Chini, R., Kreysa, E., Wink, J. E., & Salter, C. J. 1988, A&A, 191,44 Article number, page 10 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Mottram, J. C., Hoare, M. G., Urquhart, J. S., et al. 2011a, A&A, 525, A149Mottram, J. C., Hoare, M. G., Urquhart, J. S., et al. 2011b, A&A, 525, A149Mueller, K. E., Shirley, Y. L., Evans, Neal J., I., & Jacobson, H. R. 2002, ApJS,143, 469Nelson, A. F. & Marzari, F. 2016, ApJ, 827, 93Ohashi, N. 2004, in IAU Symposium, Vol. 221, Star Formation at High AngularResolution, ed. M. G. Burton, R. Jayawardhana, & T. L. Bourke, 75Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016,ApJ, 821, 82Pinte, C., van der Plas, G., Menard, F., et al. 2019, arXiv e-prints,arXiv:1907.02538Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2Pomohaci, R., Oudmaijer, R. D., & Goodwin, S. P. 2019, MNRAS, 484, 226Price, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270Richardson, K. J., White, G. J., Gee, G., et al. 1985, MNRAS, 216, 713Robitaille, T. P. 2017, A&A, 600, A11Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169,328Sanna, A., Kölligan, A., Moscadelli, L., et al. 2019, A&A, 623, A77Schreyer, K., Helmich, F. P., van Dishoeck, E. F., & Henning, T. 1997, A&A,326, 347Schwartz, P. R., Thronson, H. A., J., Odenwald, S. F., et al. 1985, ApJ, 292, 231Shu, F. H. & Adams, F. C. 1987, in IAU Symposium, Vol. 122, CircumstellarMatter, ed. I. Appenzeller & C. Jordan, 7–22Siebenmorgen, R. & Krügel, E. 2010, A&A, 511, A6Simon, M., Peterson, D., Longmore, A., Storey, J., & Tokunaga, A. 1984, inBAAS, Vol. 16, 938Simpson, J. P., Burton, M. G., Colgan, S. W. J., et al. 2009, ApJ, 700, 1488Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2011, ApJ, 730, 32Takahashi, S. Z. & Inutsuka, S.-i. 2016, AJ, 152, 184Tan, J. C., Beltrán, M. T., Caselli, P., et al. 2014, in Protostars and Planets VI, ed.H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 149Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529Thompson, R. I., Corbin, M. R., Young, E., & Schneider, G. 1998, ApJ, 492,L177Thompson, R. I. & Tokunaga, A. T. 1978, ApJ, 226, 119Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016, Nature, 538, 483Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE)Conference Series, Vol. 627, Proc. SPIE, ed. D. L. Crawford, 733Ulrich, R. K. 1976, ApJ, 210, 377Urquhart, J. S., Hoare, M. G., Lumsden, S. L., et al. 2012, MNRAS, 420, 1656Urquhart, J. S., Moore, T. J. T., Csengeri, T., et al. 2014, MNRAS, 443, 1555van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJ, 867, L14Vinkovi´c, D. & Jurki´c, T. 2007, ApJ, 658, 462Wang, Y., Beuther, H., Bik, A., et al. 2011, A&A, 527, A32Wheelwright, H. E., de Wit, W. J., Oudmaijer, R. D., et al. 2012, A&A, 540, A89Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211Whitney, B. A., Robitaille, T. P., Bjorkman, J. E., et al. 2013, ApJS, 207, 30Whitney, B. A., Wood, K., Bjorkman, J. E., & Wol ff , M. J. 2003, ApJ, 591, 1049Whittet, D. C. B., Gerakines, P. A., Hough, J. H., & Shenoy, S. S. 2001, ApJ,547, 872Wood, K., Wol ff , M. J., Bjorkman, J. E., & Whitney, B. 2002, ApJ, 564, 887Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7Zhang, Y. & Tan, J. C. 2018, ApJ, 853, 18Zinchenko, I., Liu, S. Y., Su, Y. N., et al. 2012, ApJ, 755, 177Zinchenko, I., Liu, S. Y., Su, Y. N., et al. 2015, ApJ, 810, 10Zinchenko, I., Liu, S. Y., Su, Y. N., & Sobolev, A. M. 2017, A&A, 606, L6 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 inter-esting e ff ects on each of the simulated scales when experiment-ing with the parameter space of the RT code, and 3) present theresults and final best fitting model for each object.G305.20 + Appendix A.1: W33A de Wit et al. (2010) investigated the MYSO W33A by fitting anRT model to the SED and MIDI visibilities. They found that thesource’s MIDI data and SED could be reproduced by a modelconsisting of an envelope geometry with two bipolar outflowcavities and concluded that the majority of the mid-IR emissionoriginates from the outflow cavity walls, in agreement with otherworks such as De Buizer (2006). Ilee et al. (2013) attribute theCO bandhead emission of W33A to a 1-2au disk, while Davieset al. (2010) find a Br γ jet perpendicular to a disk found inCO emission and absorption using NIR adaptive optics integralfield spectroscopy. At longer wavelengths, two large continuumsources are observed within W3, the brighter of which is MM1(Galván-Madrid et al. 2010) or W33A. Maud et al. (2017) stud-ied W33A with ALMA and SMA data to investigate whether aKeplerian disk may be present. They detected a change in po-sition of the blue-shifted and red-shifted emission close to thepeak of the continuum emission, which could be indicative of arotating disk. They probed complex large-scale structure, withhints of a disk at smaller scales. The observations used in ourstudy of the source consist of the MIDI observations from deWit et al. (2010), the COMICS image presented in Wheelwrightet al. (2012) and an SED compiled from the RMS database andde Wit et al. (2010).The luminosity of the source is listed as 3.4 × L (cid:12) by theRMS survey (Mottram et al. 2011b). The kinematic distance toW33A was thought to be 3.8kpc (Faúndez et al. 2004) but wa-ter maser measurements by Immer et al. (2013) recalculate thisas 2.4 ± µ m and the visibility at 7.5 µ mfor configuration i), taking them outside of the observed error-bars and worsening the fits. This may be because the sublima-tion radius of the source will have become smaller, presentinga more compact component for emission resulting in increasedvisibilities. The 24.5 µ m density profile was also too broad andthe silicate absorption feature of the SED became deeper and inviolation of the observed MIDI spectra.After determining that changing the distance alone could notprovide a suitable fit to the dataset, recent literature, notablythe existence of a potential disk (Maud et al. 2017), was takeninto account when improving the fit. Maud et al. (2017) stud-ied W33A with ALMA, tracing its wider environment. In theirwork they note the turbulent and messy nature of the source andimage a stream of material that is spiral-like in nature on scalesof ∼ ≤ Article number, page 11 of 34 & A proofs: manuscript no. first
Fig. A.1.
Observed visibilities for W33A for each configuration (black) with the simulated visibilities (coloured) for di ff erent models to show thee ff ects di ff erent components on the simulated MIDI data. disk value. This implies that the disk’s emission from its innerregions is more important than the emission from its outer re-gions, which is consistent as the outer regions of a dense featurelike a disk will be too cool to be traced by our high-resolutionIR data. The disk is flat, with the large-grains disk componentscale-height set to half that of the smaller grains disk. The modelenvelope has outflow cavities carved out at a 20 ◦ opening an-gle, twice that of de Wit et al. (2010), but with a similar cavitydensity. The envelope infall rates are identical between the twoworks. The inclination from de Wit et al. (2010) for the 3.8kpcdistance well satisfies all observables for the closer distance. Thevisibilities are displayed in Figure A.1, the model 24.5 µ m imageand the radial profiles are displayed in Figure A.2 and the SEDis displayed in A.3.To summarise, W33A has demonstrated the importance ofhaving the correct distance to MYSOs when deriving geometri-cal traits with the final model for the new distance characterisedby the presence of a flat disk, which was not a part of the previ-ous 3.8kpc model of de Wit et al. (2010). Appendix A.2: NGC 2264 IRS1
NGC 2264 IRS1 presented an interesting case for the fittingprocess due to its close distance of ∼ ff ected by the variation of diskparameters. In the case of NGC 2264 IRS1 however, the disk hadprofound e ff ects on the profile and added another level of com-plexity to the fitting. Through the fitting it was determined thatan Ulrich-type density distribution, as opposed to the power-lawone of de Wit et al. (2009), results in a model that produces a sat-isfactory fit to the datasets. Including the envelope is necessary to obtain a good fit to the longer wavelength data in the SED, whichwas not reproduced by previous work by Grellmann et al. (2011).The final model consists of a disk-outflow system, surroundedby this Ulrich envelope. The central protostar has a luminosityof 4.17 × L (cid:12) which is in good agreement with those found inthe literature (Harvey et al. (1977a), Schwartz et al. (1985)). Thestellar mass in the model is 9M (cid:12) , similar to the mass postulatedby Thompson et al. (1998) based on the source’s B0-B2 spectraltype. An inclination of 15 ◦ was able to successfully fit all theobservables, which is consistent with the existence of an outflowin the line of sight found by Schreyer et al. (1997).Figure A.4 shows the MIDI data and two models which arediscussed further later in this section. The two N-band config-urations of NGC 2264 IRS1 trace significantly di ff erent areasof the protostellar environment, with large di ff erences in bothPA and baseline (an illustration is provided in Figure 1). Boththe MIDI configurations predominately trace disk material. The40.2m baseline was characterised by its upward trend in visibil-ity which proved challenging to fit. Often the issue was, again,that the model visibilities were too high between 7.5-9 µ m espe-cially, indicating that the emitting regions at larger scales withinthe disk were also too compact and needed to be minimised. Anear-perfect fit for Configuration ii) could be found when a diskwas included with a 500au maximum radius consisting of 20%large-grains (Model 1 in Wood et al. (2002)) and 80% smallgrains (Galactic ISM grains as defined by Kim et al. (1994))but the fit for the 7-9 µ m region of the 40.2m baseline fit waspoor in comparison. By adding substructure to the disk by us-ing the 3D, spiral capabilities of the code, the fit for configura-tion i) could be improved (as described in Section 4). Includingthis ‘spiral’ dramatically improves the fit of Configuration i) butworsens the fit of configuration ii) as can be seen in Figure A.4. Article number, page 12 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.2.
Convolved model image (top, contours represent 5, 10, 25 and75% of the peak flux) and the radial profile (bottom) of the best fittingmodel for the 2.4kpc distance. The original COMICS image of W33Acan be found in Wheelwright et al. (2012).
Fig. A.3.
SED of the preferred final model for W33A. Data from theRMS survey are shown as diamonds, the MIDI spectrum is shown as redcrosses and data points from de Wit et al. (2010) are shown as triangles.
Fig. A.4.
Observed visibilities of NGC 2264 IRS1 for each MIDI con-figuration (black) with the simulated visibilities for each model image(coloured).
The overall change in the simulated visibilities for configurationii) ( V ∼ V ∼ / VSGs to all ar-eas of the model protostellar environment increases the visibili-ties overall, vastly worsening the fit of configuration i). However,including the PAHs in just the cavity and envelope results in nochange to the MIDI visibilities. This consecrates the fact that theMIDI data are dominated by the disk for this source. Ultimately acompromise is made by adding PAHs to the cavity and envelopeonly. This can be explained physically as PAHs constitute someof the smallest components of protostellar material and behavelike very small dust grains. As a result, they could be vulnerableto the e ff ects of photoevaporation. This means that they could bepushed from the inner environment into the envelope and cavitymaterial at greater distances from the central protostar. Addi-tionally, the PAHs in the further reaches of the envelope materialwould feel a weaker force from the stellar photons and would notbe so easily removed from the environment. Therefore, the ideathat PAHs may be present in outer regions of the protostellar en-vironment but not the inner regions, may be consistent (Sieben- Article number, page 13 of 34 & A proofs: manuscript no. first
Fig. A.5.
COMICS 24.5 µ m image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image wasconvolved with the PSF of the observed object to accurately mimic the e ff ects of the telescope specific to the observations. The contours in theimages 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 bluediamonds, the yellow diamond represents the COMICS flux density andthe fluxes corresponding to the MIDI visibilities are also shown in red. morgen & Krügel 2010). The improvements noted by the inclu-sion of the spiral exist with or without the presence of PAHs. Figure A.5 shows the COMICS image, model 24.5 µ m imageand the resulting radial profiles. The image is symmetrical andappears to show details of substructure, but de Wit et al. (2009)state that most of these are due to di ff raction e ff ects. Given itsclose distance, the 24.5 µ m emission of NGC 2264 IRS1 is heav-ily influenced by the disk included in the model. The close dis-tance, inclination and the need for substructure within the modelmeant that the disk mass for NGC 2264 IRS1 could be more ac-curately constrained than other sources to an order of ∼ (cid:12) .Changing the cavity density exponent, and therefore how evenlythe cavity material was distributed, could assist in the fitting ofthe source. Increasing the exponent from the default value of 0(constant density throughout the cavity) to 0.25 (meaning themass was slightly weighted towards the bottom of the cavity,close to the protostar), reduced the amount of 24.5 µ m emissioncreating a better fit. This is consistent as more emitting materialis at smaller radii within the protostellar environment, makingthe source appear less extended. A combination of the close-to-pole-on inclination and a foreground A v of 25 allowed su ffi cientfitting of the silicate absorption feature in the SED and is closeto A v values quoted in the literature (Thompson et al. (1998),Thompson & Tokunaga (1978)). Changing the cavity exponentto satisfy the COMICS profile flattened the peak of the SEDslightly, improving the fit. The SED fit is shown in Figure A.6. Article number, page 14 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.7.
Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured).
Appendix A.3: S255 IRS3
S255 IRS3 has become of particular interest since an accretion-ejection event was observed around the source by Caratti oGaratti et al. (2017). A particularly bright 6.7GHz methanol flarewas discussed in Fujisawa et al. (2015) and a sub-millimetre flarewas also found in Zinchenko et al. (2017). This flare was latertied to NIR and radio flares by Caratti o Garatti et al. (2017)and these were determined to be the result of an accretion burstof order 10 − M (cid:12) yr − . This was the first directly observed accre-tion burst around an MYSO, which has important implicationsas episodic accretion is a key component of the low-mass starformation process (Stamatellos et al. 2011).All the interferometric data for S255 IRS3 came from Bo-ley et al. (2013), with six configurations of data available forfitting. All the configurations, save one, trace very similar posi-tion angles of around 70 ◦ . The remaining configuration traced aposition angle of 32.6 ◦ and was low resolution due to its 17.6mbaseline.The COMICS image from de Wit et al. (2009) is used in thefitting process. Both IRS1 and IRS3 are visible in the COMICSimage. In order to exclude scales that may be influenced by the emission of IRS1 we fit the inner 3" of the source as opposed tothe whole profile. All of the fitted data were taken at a time thatappears to be an accretion lull for S255 IRS3, as the fluxes fromthe SED (which were taken at similar times to the MIDI andCOMICS data) are in much better agreement with the non-burstfluxes from Caratti o Garatti et al. (2017) than the burst fluxes.Fluxes from the RMS survey and Caratti o Garatti et al. (2017)were incorporated (blue diamonds in Figure A.9), alongside ad-ditional fluxes from de Wit et al. (2009) (green crosses in FigureA.9) and the MIDI flux spectrum.The plots for the final model are shown in Figures A.7 - A.9.The final luminosity for the central source in the model is 2.15 × L (cid:12) , which is in excellent agreement with that observed byWang et al. (2011) of ∼ × L (cid:12) . The overall model is a diskand outflow system, in agreement with works such as Zinchenkoet al. (2017), who studied the source with SMA and IRAM-30mcontinuum and spectral lines and find signatures of rotation andattribute this to a disk-like structure. The presence of the diskalso agrees with the detection of the accretion event; as for low-mass stars, such events occur as a result of accretion from a disksurrounding the protostar. The fact that the envelope is rotatingand infalling could be linked to the fact that the source has dis- Article number, page 15 of 34 & A proofs: manuscript no. first
Fig. A.8.
COMICS 24.5 µ m image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image wasconvolved with the PSF of the observed object to accurately mimic the e ff ects of the telescope specific to the observations. The contours in theimages represent 5, 10, 25 and 75% of the peak flux. played an accretion event. Infall from the wider environment isdirectly linked to disk accretion in the turbulent core model (Tanet al. 2014). With a centrifugal radius of 500au, as opposed to Fig. A.9.
Model SED of the best-fitting model (black) for S255 IRS3.The di ff erent symbols correspond to di ff erent datasets, purple crossesare data from Caratti o Garatti et al. (2017), blue diamonds are fluxesfrom the RMS database, squares are data from de Wit et al. (2009) (withthe COMICS flux in yellow) and the fluxes corresponding to the MIDIvisibilities are also shown in red. Open symbols correspond to upperlimits. − gcm − , which is comparable to that of W33A.An inclination of 120 ◦ fit all the observables, which is con-sistent with the close-to-edge-on orientation suspected by Boleyet 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 from0.1-10M (cid:12) could be present and satisfy the observables whichis poorly constrained due to the limitations of the IR as a diskmass tracer. A compromise was found for the longer wavelength
Article number, page 16 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.10.
Observed visibilities of IRAS 17216-3801 for each configuration (black) with the simulated visibilities for each model image (coloured)(1). 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 notlist as upper limits).Varying the disk parameter space has negligible e ff ects onthe COMICS profile, save for the disk outer radius and the cen-trifugal radius of the envelope, which is attributed to the sourcesinclination and distance. It was found that the 2000au maximumdisk radius and centrifugal radius used for G305 resulted in toomuch 24.5 µ m emission to successfully fit the radial profile, im-plying that the disk must be less extended. Reducing these valuesto 500au allowed a good fit for the COMICS profile.The wider environment of S255 IRS3 di ff ers compared to theother MYSOs in this sample, with an outer envelope radius thatis five times smaller than the other sources and an envelope infallrate that is significantly lower. Raising the infall rate to a valuecomparable to those found for the other sources ( ∼ − M (cid:12) yr − )increases the peak of the SED significantly. While this appears toimprove the fit to the SED, the fluxes that are better fit are almostcertainly detecting flux from both IRS1 and IRS3 and as a result this envelope infall rate is likely to be an overestimation of whatis actually present. A fit with the lower envelope infall rate thatfits the smaller aperture ∼ µ m point better is preferred. High-resolution, small aperture observations at ∼ µ m could allowseparate SEDs for the two MYSOs to be obtained and for thisto be better investigated. While obtaining and analysing theselonger wavelength data are outside the scope of this work, in-creasing the envelope infall rate also lowers all the simulatedMIDI visibilities that only trace IRS3, not IRS1. Using an en-velope infall rate of the same order of magnitude as W33A andNGC 2264 IRS1 vastly worsens the fits for all configurations,further showing the need for a lower envelope infall rate. Appendix A.4: IRAS 17216-3801
Kraus et al. (2017) resolved IRAS 17216-3801 into a close bi-nary using spectro-interferometric data from AMBER, CRIRESand GRAVITY. They find that the two stars are separated by ∼ Article number, page 17 of 34 & A proofs: manuscript no. first
Fig. A.11.
Observed visibilities of IRAS 17216-3801 for each config-uration (black) with the simulated visibilities for each model image(coloured) (2). of stars but disks surrounding each individual object as well.HO-CHUNK3D does not provide the capabilities for modellingmultiple systems but does present the opportunity to investigatethe circumbinary material of IRAS 17216-3801. A single centralobject is included in the model, whose luminosity was obtainedthrough the fitting of the SED as per the other sources in thiswork.Surrounding this central source is a disk-outflow-envelopesystem. The disk scale-heights are very similar to the defaultssuggested by Whitney et al. (2013), with the small grains diskbeing slightly flatter. The inner disk radius is 100au, which iscomparable to the estimate of Kraus et al. (2017) for their cir-cumbinary disk. The visibilities are shown in Figures A.10 andA.11. IRAS 17216-3801 has the largest MIDI dataset of our sam-ple with ten configurations. The position angles of this datasetrange from approximately 14-175 ◦ but all are fairly low resolu-tion with baselines of 35m or shorter. This means that most of theMIDI configurations are sensitive only to the circumbinary mate-rial. As illustrated in Figure 1, configurations viii), ix) and x) arethe only three with PAs suitable for tracing the circumprimarydisk and only configuration viii) has the potential to trace thecircumsecondary disk. The baselines of configurations viii), ix)and x) are some of the largest for this source; configuration viii)can probe scales of ∼ ∼ ◦ ) and a cavity density exponent of1. This is because the VISIR profile is barely resolved. A den-sity exponent of 1 as opposed to 0 (the default value) meansthat the density in the cavity decreases linearly with radius in-stead of remaining at a constant value. The innermost part of thecavity will then be brighter than its outer regions, reducing theextended emission and lowering the ∼ µ m flux at larger radii,improving the fit. Similarly, close to pole-on inclinations reducethe extended emission as only one of the cavity lobes is directlyvisible in the line of sight.The SED is shown in Figure A.13. The ATLASGAL 870 µ mpoint proved di ffi cult to fit. The aperture of this measurementwas 19.2" but the region is fairly uncrowded with no obvioussources contaminating such a measurement. In an attempt to im-prove the fit to this data point, disks of bigger radius were tri-alled. The larger the extent of the disk, the more cold emission Article number, page 18 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.12.
VISIR 19.5 µ m image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image wasconvolved with the PSF of the observed object to accurately mimic the e ff ects of the telescope specific to the observations. The contours in theimages represent 5, 10, 25 and 75% of the peak flux. at sub-millimetre and millimetre wavelengths will be present.A 4000au maximum radius disk showed small improvements inthe millimetre of the regions but more substantial changes in the20 µ m region of the SED, making the dip-like feature between10-20 µ m flatten. This increase in radius caused all the simulatedMIDI visibilities to decrease, worsening some fits and but im-proving the fit of configuration ix). The fit of the VISIR profilewas also violated as the model radial profile displayed more ex-tended emission than the observed profile. This occurred for anydisks larger than 2000au so this was the final value chosen forthe disk to preserve the quality of the VISIR fit. Making the en-velope radius larger improved the fit to the longer wavelengthend of the SED by shifting the SED to longer wavelengths, butthis worsened the fit of the 20 µ m region of the SED.The mass of the central source was also found to a ff ect theshape of 10-20 µ m region of the SED. For example, if a 20 solarmass central source is used this 10-20 µ m region of the SED isflatter. This is consistent as the density distribution of the Ulrichenvelope is proportional to M / ∗ (Terebey et al. 1984), mean-ing that as the stellar mass increases, the density will decreasewithin the infall radius and reduce the amount of 20 µ m flux.Configuration v) proved the most di ffi cult to fit. This occurreddespite the fact that it is one of the lowest-resolution datasetsand is likely to be tracing material within the circumbinary disk.A large (nearly eight-fold) upturn is present in the visibilities be- Fig. A.13.
Model SED of the best-fitting model (black) of IRAS 17216-3801. Multi-wavelength flux measurements from the literature are rep-resented as blue diamonds, the yellow diamond represents the VISIRflux density and the fluxes corresponding to the MIDI visibilities arealso shown in red. tween 7 and 9 µ m, which is often a hallmark feature of a compactcomponent like a disk. To test whether the very large increase invisibility could be due to the presence of the inner rim specif- Article number, page 19 of 34 & A proofs: manuscript no. first ically, the inner disk radius was changed to the disk to 134au,which is the scale traced by this configuration at 7.5 µ m. A minorupturn in the visibilities was observed but the fits for the other,higher-resolution configurations worsened. When the inner ra-dius of the disk was varied between 100-150au the visibilitieschanged marginally but the overall goodness of fit changed neg-ligibly.The upturn in configuration v) could be due to the binarysystem interacting with the circumbinary disk. Small-scale fea-tures, such as a pu ff ed-up inner rim and disk-spiral, were testedbut also present marginal changes, mostly likely due to the ob-ject’s distance. A very specific local substructure must be presentat the scales traced by this baseline alone and not the others ofsimilar PA / baseline to induce the observed increase in visibil-ity. Various theoretical works have shown that binary or multiplestellar systems can result in the formation of such a small-scalesubstructure in protostellar disks. Nelson & Marzari (2016) usedthe SPH code VINE to model the circumbinary disk of the low-mass source GG Tau A. They found that the binary system ‘stirs’the material in the disk, leading to the formation of spiral struc-tures, which extend towards the centre within the inner radiusof the disk as streams in the earliest stages. Within the disk it-self, they see significant areas of over-density. Thun et al. (2017)find that the inner rims of the circumbinary disks are often unsta-ble, clumpy and recessing from their hydrodynamical modelling.Desai et al. (2019) focus on the study of unequal mass binariesin particular using a 3D hydrodynamics code CHYMERA. Herethey find that gravitational instabilities are induced by the binaryinteraction, with the smaller of the companions creating a den-sity wave at the earliest stages and at later stages also see thegeneration of spirals. Observationally, Tobin et al. (2016) detectspiral structures surrounding a three-star system. Although lowmass, two of the stars have a separation comparable to that ofIRAS 17216-3801 of ∼ ff ect the stability and geometry of its circumbinary disk will beincreased if those stars have higher mass. As such, the likelihoodof over-densities that could create the compact feature requiredto fit this configuration is high and could therefore explain thepoor fitting of these data. Appendix A.5: Mon R2 IRS2
The Mon R2 region is home to five di ff erent sources and lo-cated at a distance of ∼ µ m maps (Gian-nakopoulou et al. 1997) and the 24.5 µ m image of de Wit et al.(2009). Mon R2s IRS2 lies within the centre-north region ofthis H ii region and is its main source of illumination (Aspin &Walther 1990). A large CO bipolar outflow is present in the widerenvironment (Giannakopoulou et al. 1997) but it is unknownwhich MYSO in the region is responsible for it. Another mas-sive protostellar system, Mon R2 IRS3, which has been resolvedinto two sources A and B, is the dominant source of the wholeregion. No free-free emission is associated with IRS2 and it isunresolved in UKIRT observations (Alvarez et al. 2004). Boleyet al. (2013) find its mid-IR emission to be compact ( ∼ ∼ (cid:12) . Recent work by Jiménez-Serra et al. (2020) studies the gas of Mon R2IRS2. They find that their ALMA images, studies of the H21 α emission and 3D radiative transfer modelling show the presenceof an ionised Keplerian disk and a neutral rotating disk. The diskis highly flared (as shown in Figure 4 of the cited paper) anddeviation of the H21 α emission centroids from the disk planeimply warping in the disk, which they suggest could be due tothe presence of a secondary object.Our work finds that Mon R2 IRS2 is another source wherethe fits to the MIDI data are improved by the inclusion of aninner hole (20au in radius). The dust sublimation radius of thesource is small ( ∼ ff ects on the MIDI fits, implying that the N-bandemission of the source is disk-dominated. Removing PAHs fromthe entire environment greatly worsens the NIR slope of the SEDfit. The final luminosity of 5520L (cid:12) is similar to that of Henninget al. (1992). The envelope of Mon R2 IRS2 is similar to that ofthe majority of the sample, with a comparable infall rate and thesame outer radius. The cavity density is low ( ∼ − gcm − ) andthe cavity opening angle is 30 ◦ . The source is inclined at 130 ◦ so both cavity lobes are visible. Through the SED fitting it wasdi ffi cult to simultaneously fit the silicate absorption feature andthe 20 µ m flux. Eventually the fitting of the silicate absorptionfeature was prioritised as this is the dataset with higher spatialresolution. The visibilities in Figure A.14, the images in FigureA.15 and the SED in Figure A.16.The position angles of the MIDI baselines for Mon R2 IRS2mean that all five configurations are likely to be tracing diskmaterial and configuration i) is likely to be tracing cavity ma-terial. Three of the MIDI configurations have high visibilitiesand are tracing an unresolved component of the protostellar envi-ronment. The simulated visibilities of configuration i) get lowerif the outer radius of the disk is made larger implying that asmall, compact disk best satisfies this configuration. However, as R c / R diskmax increased to values of 4000au (in agreement with someof the largest observed disks around MYSOs e.g. Johnston et al.(2015)) improvements were seen to the fit of the COMICS pro-file as the flux at larger radii increased, so a compromise ulti-mately had to be made. Configurations iv) and v) of the MIDIdata proved the most di ffi cult to fit and resulted in the inclusionof the hole. Without this the visibilities for both the configura-tions were too high. This hole is 40au across, in agreement withthe size of the compact emission from Boley et al. (2013). Thisimplies that the mid-IR emission for this source, as traced bytheir work, is mostly dominated by emission from this inner rim.A compromise was found between the fitting of configurationsiv) and v) as opposed to fitting one perfectly and leaving theother poorly fit. Appendix A.6: M8EIR
This work finds that the 13.5M (cid:12) central object of Ilee et al.(2013) well reproduces the shape of the SED of M8EIR. A lu-minosity of ∼ (cid:12) is su ffi cient to fit the multiple datasets,which is also in good agreement with the value of L obtainedfor the source in Ilee et al. (2013). In the preferred model a5 × au envelope, infalling at a rate of 1 × − M (cid:12) yr − , sur-rounds this protostar. The presence of an envelope agrees withprevious work by Simon et al. (1984) who interpreted M8EIRas two physically di ff erent components,a hot compact one and a Article number, page 20 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.14.
Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for Mon R2 IRS2. broader cooler one, and stated that the broad component is notmolecular cloud material, but an ‘intermediary phase’ of mate-rial between molecular cloud and hot disk material. A disk ispresent with a maximum radius of 2000au and outflow cavitiesare carved out of the envelope with opening angles of 25 ◦ anddensities of order 10 − gcm − . Using a cavity density exponentof 0.25 was found to improve both the simulated 24.5 µ m profileand the shape of the SED, similar to NGC 2264 IRS1 and IRAS17216-3801. Using an inner hole (radius 30au) again improvedthe fit of the simulated visibilities. Using the dust sublimationradius ( ∼ µ m. Such aquick visibility change can be caused by the turnover inducedby the inner rim of the disk. When varying the inner disk radiusaround 30au similar shapes were observed, but alongside furtherdetriments to the other visibility datasets. As a result 30au waschosen as the final inner radius. The simulated visibilities be-tween 9-12 µ m for configuration ii) are ∼ ff wouldbe expected, but given that an inner hole exists the idea that themedium could also not be perfectly smooth is not unlikely. Acompromise was found between the fitting of configurations iii)and iv) as opposed to fitting one perfectly and the other verypoorly. The visibilities are shown in Figure A.17, the images inFigure A.17 and the SED in Figure A.19. Appendix A.7: AFGL 2136
AFGL 2136 IRS1 (hereafter AFGL 2136) was studied by de Witet al. (2011) using a similar analysis to our own. The main di ff er-ence was the amount of MIDI data available as only one configu-ration was fit in that work, as opposed to the four used here. Ad-ditional technical di ff erences between our fitting and de Wit et al.(2011) are the multiplication of the model images with a Gaus-sian and the convolution of the COMICS image with a Gaussian Article number, page 21 of 34 & A proofs: manuscript no. first
Fig. A.15.
COMICS 24.5 µ m image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image wasconvolved with the PSF of the observed object to accurately mimic the e ff ects of the telescope specific to the observations. The contours in theimages 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 beseen 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 rep-resented as blue diamonds, the yellow diamond represents the VISIRflux density and the fluxes corresponding to the MIDI visibilities arealso shown in red. The unfilled diamonds represent the fluxes that werenot considered in the fitting due to their suspected contamination andthe fluxes from Henning et al. (1992) are included as squares.
PSF as opposed to an observed one (see Paper I). The COMICSimage was presented in de Wit et al. (2011) but not includedin 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 wasfound that without multiplying by the Gaussian, a good fit couldnot be found to the MIDI configuration common to both this pa-per and de Wit et al. (2011). The model also provided poor fitsto all the other MIDI configurations that were not included in deWit et al. (2011).Maud et al. (2018) find that the dust emission implies a disk-like structure, so various disk-outflow-envelope geometries weretrialled during the modelling. disk masses ranging between 10-30 solar masses are workable components of the model, in agree-ment with the kinematic measurements of Maud et al. (2018)and Maud et al. (2019). The larger the mass, the less the 24.5 µ memission, improving the fit of the radial profile. For a 30 solarmass disk the millimetre slope of the SED fit is worse than for a10 solar mass disk. The quality of the MIDI fits are comparableacross the disk mass range. For the 30 solar disk mass model, thefit of configuration ii) improves slightly as the simulated visibil-ity at 8 µ m increases, but the visibilities for configuration iv) alsoincrease, which further worsens the fit to that configuration by acomparable amount. Given the similarity across the MIDI data a Article number, page 22 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.17.
Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for M8EIR. compromise was made between the SED and COMICS fits and20 solar masses was selected as the final disk mass.Similarity is found between the preferred best fitting modelof our work and the model of de Wit et al. (2011) in that theminimum dust radius is larger than that of the dust sublimationradius. In this work the inner radius is four times larger thanthe sublimation radius, in de Wit et al. (2011) it was five. In deWit et al. (2011) a gaseous disk and bloated star resided withinthis cleared inner hole. Including a gaseous disk is not part ofthe methodology of this work, but it could exist within the innerdust rim. As discussed in the examination of G305.20 + ◦ , and the cen-tral object is similar to that of de Wit et al. (2011); a bloated starof radius 25R (cid:12) with an overall luminosity of 1.51 × L (cid:12) . Thisluminosity matches the O-type classification ascribed by Maudet al. (2018) and Maud et al. (2019) and in good agreement withthe 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 istightly constrained at 22.5 ◦ to maintain the balance of fits of thethree observables and the cavity density is of order 10 − gcm − .The centrifugal radius is 2000au and the infall rate is of order10 − solar masses per year.The visibilities of the final model are shown in Figure A.20,the images are shown in Figure A.20 and the SED is shownFigure A.22. Configuration ii) of the MIDI data proved nigh-impossible to well fit. Despite this, the preferred model recreatedthe shape of the visibilities well and was only V = ff invisibility. This implies that the protostellar environment tracedby this configuration is larger than the model environment. ThePA of this configuration was 14 ◦ meaning that, taking into ac-count the source’s PA on the sky, it could be tracing envelopematerial close to the cavity wall or disk material. Both the con-tinuum and CO contours of Maud et al. (2018) show the mediumto be clumpy, despite the disk’s stability. Given that only smoothgeometries were used in the modelling, it could be that the modelis not reproducing this emission leading to this poor fitting. Theother three MIDI configurations are matched well by the finalgeometry. Article number, page 23 of 34 & A proofs: manuscript no. first
Fig. A.18.
COMICS 24.5 µ m image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image wasconvolved with the PSF of the observed object to accurately mimic the e ff ects of the telescope specific to the observations. The contours in theimages 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 thespectrum from Feldt et al. (2008). This was sampled electronically andthe data are represented as purple crosses. The fluxes corresponding tothe MIDI visibilities are shown in red. Open symbols correspond toupper limits.
A further study at longer wavelengths, Maud et al. (2019),detect a ring-like structure at scales of 120au. This is in verygood agreement with the inner dust radius found in our model,which is four times the sublimation radius or ∼ Appendix B: List of MIDI observationsAppendix C: List of SED fluxes
Article number, page 24 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Fig. A.20.
Observed visibilities for each configuration (black) with the simulated visibilities for each model image (coloured) for AFGL 2136.Article number, page 25 of 34 & A proofs: manuscript no. first
Fig. A.21.
COMICS 24.5 µ m image (top left), convolved model image (top right) and subsequent radial profiles (bottom). The model image wasconvolved with the PSF of the observed object to accurately mimic the e ff ects of the telescope specific to the observations. The model is moresymmetric than the observed source meaning it has more emission in its eastern regions, explaining the larger amount of 24.5 µ m flux visible forthe model in the radial profiles. The contours in the images represent 5, 10, 25 and 75% of the peak flux.Article number, page 26 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs Table B.1.
List of MIDI observations of for the sample. The starred dataset is data from Grellmann et al. (2011). The data for W33A are from deWit et al. (2010). All other data are from Boley et al. (2013).
Object Configuration Date Telescopes Projected baseline Position angle ESO Run ID(UTC) (m) ( ◦ )G305.20 + / / / / Article number, page 27 of 34 & A proofs: manuscript no. first
Table C.1.
Fluxes available for use in the SED fitting of W33A.
Origin Wavelength ( µ m) Flux (Jy)2MASS 1.235 (1.21 ± × − ± × − ± ± ±
312 33.3 ± ±
412 23.2 ±
514 50.0 ±
621 144 ± ± ± ± ± ± Table C.2.
Fluxes available for use in the SED fitting of NGC 2264 IRS1.
Origin Wavelength ( µ m) Flux (Jy) Reference2MASS 1.24 0.040 ± ± ± ± ± ± ± ±
614 148 ±
921 13.7 ± ±
33 de Wit et al. (2009)KAO 53 980 ± ± ± ± Article number, page 28 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Table C.3.
Fluxes available for use in the SED fitting of S255 IRS3.
Origin Wavelength ( µ m) Flux (Jy) ReferenceNICMOS 2 0.00283 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± / IRC 8.6 103 ± ± ±
10 Caratti o Garatti et al. (2017)160 760 ± ±
150 Poglitsch et al. (2010)160 849 ± ± ± ± ± ± × − Skrutskie et al. (2006)1.66 (6.46 ± × − ± × − WISE 3.35 0.689 ± ± ± ± ± ±
514 154 ±
921 213 ± ±
20 de Wit et al. (2009)UKIRT 350 230 ± ± ± ± ± ± ± Article number, page 29 of 34 & A proofs: manuscript no. first
Table C.4.
Fluxes available for use in the SED fitting of IRAS 17216-3801.
Origin Wavelength ( µ m) Flux (Jy) Reference2MASS J-Band 1.24 0.0134 ± ± ± ± ± ± ± ± ± ± ± ± ± ±
10 This workATLASGAL 870 35.1 ± ± ± ± ± Table C.5.
Fluxes available for use in the SED fitting of Mon R2 IRS2.
Origin Wavelength ( µ m) Flux (Jy) ReferenceUKIDSS H-band 1.662 0.002 (2 ± × − Lawrence et al. (2007)K-band 2.159 0.0771 ± ± ± ± ± ±
122 Henning et al. (1992)1300 2640 ± Table C.6.
Fluxes available for use in the SED fitting of M8EIR (1).
Origin Wavelength ( µ m) Flux (Jy) Reference2MASS J-Band 1.24 (3.05 ± × − Skrutskie et al. (2006)H-Band 1.66 2.62 ± ± ± ± ± ± ± ± ±
20 de Wit et al. (2009)IRTF 1.24 (6.7 ± × − Simon et al. (1984)1.66 0.83 ± ± ± ± ± ± ±
42 Kastner et al. (1992)450 42.4 ± ± ± ± Article number, page 30 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Table C.7.
Fluxes available for use in the SED fitting of M8EIR (2).
Origin Wavelength ( µ m) Flux (Jy) ReferenceISO SWS 16.4 232 ±
70 Feldt et al. (2008)ISO LWS 2.34 11.0 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± Article number, page 31 of 34 & A proofs: manuscript no. first
Table C.8.
Fluxes available for use in the SED fitting of M8EIR (3).
Origin Wavelength ( µ m) Flux (Jy) ReferenceISO LWS 20.1 227 ±
45 Feldt et al. (2008)21.5 260 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± Article number, page 32 of 34. J. Frost, R. D. Oudmaijer, W. J. de Wit and S. L. Lumsden: A multi-scale survey of MYSOs
Table C.9.
Fluxes available for use in the SED fitting of AFGL 2136.
Origin Wavelength ( µ m) Flux (Jy) Reference2MASS J-Band 1.24 1.59 × − Skrutskie et al. (2006)H-Band 1.66 (8.52 ± × − Ks-Band 2.16 0.123 ± ± ± ± ± ±
814 270 ± ± ±
10 de Wit et al. (2009)IRAS 60 4810 ±
900 Beichman (1988)100 5700 ± ±
50 Mueller et al. (2002)450 72.7 ± ± ± Article number, page 33 of 34 & A proofs: manuscript no. first