Free Access
Issue
A&A
Volume 540, April 2012
Article Number A89
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201118689
Published online 02 April 2012

© ESO, 2012

1. Introduction

Massive stars, M ≳ 8   M, can play a dominant role in the evolution of their host galaxies via feedback phenomena such as ionizing radiation, strong stellar winds and supernova explosions. These processes exert a large influence on the outcome of subsequent star formation (SF). Massive SF can therefore be considered as an integral part of the overall SF process and a key ingredient for a complete description of it. However, our comprehension of many aspects of massive SF is still unsatisfactory. For example, it is still not known whether massive SF is intrinsically different to the generally much better understood process of SF at low masses or not (see e.g. Zinnecker & Yorke 2007).

The uncertainties regarding massive SF are in part the result of the short Kelvin-Helmholtz (KH) timescale associated with massive stars: they contract rapidly to sustain a high luminosity while still acquiring their mass. If the accretion proceeds spherically, the ensuing radiation pressure could limit the final stellar mass to approximately 40 M (Larson & Starrfield 1971; Kahn 1974; Wolfire & Cassinelli 1987). This led to the suggestion that massive SF is intrinsically different to the low mass case. However, recent radiation-hydrodynamical simulations clearly demonstrate that the radiation is actually channelled through optically thin polar cavities while accretion continues through an equatorial disc (Yorke & Sonnhalter 2002; Krumholz et al. 2009; Kuiper et al. 2010). Therefore, deviations from spherical symmetry, which are suggested by the spectral energy distributions (SEDs) of young massive stars (Harvey et al. 1977; Guertler et al. 1991), remove the limit imposed by spherical accretion. Nonetheless, confirming that massive stars build up their final mass via disc accretion by direct observations is challenging. Massive young stellar objects (MYSOs) are generally located at kpc distances and their AU-scale circumstellar environment remains unresolved with single dish telescopes at most wavelength regions (see e.g. Cesaroni et al. 2007).

Long-baseline infrared interferometry is one of the few techniques which offers the milli-arcsecond (mas) resolution required to probe for MYSO discs (see e.g. Follert et al. 2010; Vehoff et al. 2010; de Wit et al. 2011; Grellmann et al. 2011). In particular, an exemplary result is delivered by Kraus et al. (2010) using the Very Large Telescope Interferometer (VLTI) and the near-infrared beam-combiner AMBER. These authors reconstruct a K-band synthesis image of a 20 M MYSO with an angular resolution of 2.4 mas. The resultant image shows hot material in a disc-like geometry of approximately 20 AU in size. As a result, a consensus is beginning to emerge that the disc accretion scenario can account for the formation of stars up to at least  ~30   M. This implies that the MYSO circumstellar environment deviates significantly from spherical symmetry. In this scenario, rotation, accretion discs and outflow activity play a crucial role in providing the asymmetric environment required to facilitate accretion (see e.g. Beuther et al. 2002).

The combination of a high stellar luminosity with accretion and outflow activity results in the circumstellar environment of MYSOs being shaped by several forces. Determining the relative importance of these phenomena in sculpting the circumstellar environments of MYSOs requires spatially resolved observations. In particular, the mid-infrared (MIR) wavelength region provides a wealth of information on the geometry of the accretion environment. Furthermore, this information can be accessed by both interferometric techniques and diffraction limited single-dish imaging with 8 m class telescopes. MIR interferometric observations reveal that the N-band emission of MYSOs on scales of 50 to 100 AU is dominated by warm dust located in the envelope (de Wit et al. 2007; Linz et al. 2009; Vehoff et al. 2010). It is expected that such envelopes will contain outflow cavities as there are many detections of outflow activity associated with MYSOs (e.g. G35.20–0.74 in De Buizer 2006; IRAS 20126+4104 in De Buizer 2007; Cep A HW2 in de Wit et al. 2009). In de Wit et al. (2010), we argue that, specifically, the warm dust responsible for the N-band emission of W33A is located in the envelope close to the walls of the cavities evacuated by the outflow.

The aim of this paper is twofold. We aim to assess the source of the MIR emission of a sample of MYSOs and thus test the premise that their MIR emission is dominated by warm dust located in outflow cavity walls. This will then allow us to determine whether outflows play an important role in shaping the environments of MYSOs (as opposed to, for example, rotation). We address these goals by comparing spatially resolved images of a sample (20 objects) to appropriate models. Our approach is similar to the work presented previously in DW09: de Wit et al. (2009, hereafter DW09) with the addition of 2D modelling. We present MIR VLT images observed at 20   μm for a sample of MYSOs. We then compare the spatially resolved images to two-dimensional, axis-symmetric dust radiative transfer models that feature rotating, collapsing envelopes and outflow cavities (the models of Whitney et al. 2003a,b). The paper is structured as follows. We describe the sample selection, the observations and data reduction procedures in Sect. 2. The MIR images are described in Sect. 3 and the radiative transfer model and its use are detailed in Sect. 4. The results of the modelling are presented in Sect. 5 and discussed in Sect. 6. Finally, our conclusions are given in Sect. 7.

2. Selection, observations and data reduction

We aim to perform our analysis on a representative sample of MYSOs, in order to draw general conclusions. Selecting such a sample is not trivial. Initial attempts employed IRAS data (see e.g. Palla & Stahler 1991; Molinari et al. 1996; Sridharan et al. 2002). As a result, such catalogues suffered from source confusion due to the large beam of IRAS. To rectify this problem and generate a Galaxy-wide, unbiased, MYSO sample, we have conducted a survey to detect and characterise MYSOs: the RMS survey (Lumsden et al. 2002; Urquhart et al. 2008). The RMS is based on data of the MSX survey (Egan et al. 2003). MSX data offers a significant improvement over IRAS data in terms of resolution (e.g. arcsecond rather than arcminute resolution), and thus enables the selection of a representative sample of MIR bright MYSOs.

From the RMS database1 we selected objects with distance estimates within 3 kpc2 and with 21 μm MSX fluxes larger than 30 Jy. The flux limit was imposed in order to obtain a decent signal-to-noise ratio in the wings of the resolved profiles. In total, we observed 19 objects at 20 μm using the VLT Imager and Spectrometer for the mid InfraRed (VISIR, Lagage et al. 2004) mounted at the Cassegrain focus of UT3 of the VLT (see Table1). Observations were conducted using the imaging mode of the instrument and the Q3 filter which has a central wavelength of 19.5 μm and a half-band-width of 0.4 μm. During the observations, VISIR was equipped with a DRS 256  ×  256 Si:As detector with an angular pixel size of 0.127′′. This configuration enabled us to obtain oversampled, diffraction limited images, which with 8 m class telescopes results in an angular resolution of approximately 0.6′′.

In this paper, we also present the imaging observations and analysis of a key MYSO object W33A (G012.9090–00.2607). It was observed at the slightly longer wavelength of 24.5 μm with the COMICS instrument (COoled Mid Infrared Camera and Spectrometer, Kataza et al. 2000) mounted on the Cassegrain focus of the Subaru telescope. We employed the imaging facility of COMICS and used a filter centred at 24.5 μm (see DW09 for filter response functions). During the observations, COMICS was equipped with a 320  ×  240 Si:As IBC detector which provides oversampled diffraction limited images with a pixel size of 0.13′′. A log of all the observations is presented in Table 1.

For both the VISIR and COMICS observations, standard stars were observed to provide a reference point-spread-function (PSF). The standard stars selected are MIR bright,  ~100 Jy at 25 μm, isolated sources which are expected to be unresolved. Data reduction, consisting of shifting and adding nodded and chopped images, was conducted using routines written in idl. The resultant images were not astrometrically corrected. Therefore, the images are presented in terms of distance from the peak of the intensity distribution of the target source. The azimuthally averaged intensity profiles of the standards are displayed in Fig. 1. The radial profiles of the VISIR standard stars are generally consistent until a level of approximately 0.001 times the peak flux and appear to behave as point sources. The standard star of the COMICS data (HD 124897) is of lower signal-to-noise ratio but is similar out to 1 percent of the peak flux to the standards presented in DW09.

thumbnail Fig. 1

The azimuthally averaged intensity profiles of the PSF standard stars. HD 124897 is the PSF standard taken with Subaru/COMICS. The solid line with square points is the average VISIR PSF.

Open with DEXTER

Table 1

Log of the observations.

3. Observational results

We assessed the MIR images of our target sources by comparing their azimuthally averaged intensity profiles to the instrumental PSF determined from the standard stars. We find that the majority of the MYSOs observed, 14 of 20, are spatially resolved. The images of the target MYSOs are presented in Figs. 2 and 3, which contain the unresolved and resolved MYSOs respectively.

The resolved MYSOs generally appear as a single source within the VISIR field of view, 30″  ×  30″. Most display an approximately circular symmetric morphology at the resolution of the observations (~0.6′′), although some show clear signs of structure (e.g. G263.7759−00.4281, Fig. 3: panel A). Several objects are associated with additional, distinct sources of localised emission, which often have a cometary type morphology (e.g. G269.1586−01.1385 & G349.7215+00.1203, Fig. 3: panels D & M). This is suggestive of a Hii region (Hoare et al. 2007). Indeed, both the objects mentioned above are detected at radio wavelengths at coordinates consistent with the locations of the cometary 20 μm sources (see Urquhart et al. 2007). It is of interest to note that 24.5 μm images of regions containing both MYSOs and ultracompact Hii regions show a clear morphological difference between the two types of source, i.e. compact vs. extended (DW09).

thumbnail Fig. 3

Resolved sources. The images are scaled logarithmically. North is to the top of the page and east is to the left. The contours typically represent 2, 5, 10, 25 and 75 percent of the peak flux.

Open with DEXTER

Approximately one third of the sample are unresolved. Using the Kolmogorov-Smirnov test, we investigated whether the resolved and unresolved samples were drawn from intrinsically different distributions of MYSOs in terms of distance and luminosity. There is no significant difference between the distance and luminosity distributions of the two samples, as might be expected given that the objects are drawn from a single catalogue. However, the angular size of a centrally heated source of flux is dependent upon both the luminosity of central object and the distance to it (angular size  ∝ , see e.g. Vinković & Jurkić 2007). We find that the resolved objects generally have a higher value of than the unresolved objects. More specifically, the hypothesis that the unresolved sources have the same distribution as the resolved sources can be discarded at a level of 98 percent significance. Thus, we conclude that resolved sources appear larger than the unresolved sources as the resolved sources are generally more luminous and less distant than their unresolved counterparts.

4. Radiative transfer modelling

To investigate the nature of the circumstellar structure of the resolved sources, we compare the observed images to 2D axis-symmetric radiative transfer models appropriate for massive stellar objects surrounded by in-falling, rotating envelopes. By simultaneously comparing the resolved MIR images and SEDs to the radiative transfer calculations, we aim to understand the source of the 20 μm MIR emission of MYSOs. The models and the modelling procedure are described in detail in the following sections.

4.1. The radiation transfer code

We employ the 2D, axis-symmetric dust radiation transfer code of Whitney et al. (for details see Whitney et al. 2003a,b). The code calculates radiation transfer through a dusty structure which consists of a proto-stellar envelope surrounding a central star. A low density polar cavity can be inserted into the overall envelope structure. The code also allows for a dusty circumstellar disc. The density distribution of the envelope is described by the treatment of a collapsing, rotating envelope of Ulrich (1976) and Terebey et al. (1984). The key parameters determining the envelope density are the mass infall rate and the centrifugal radius. The dust sublimation radius is calculated self-consistently from the stellar luminosity and the dust sublimation temperature. The shape of the polar cavities can be described in two ways: using a polynomial function or following streamlines. The streamline is conical on large scales, while the apparent opening angle of the polynomial function can be specified (the opening angle at the stellar surface is 180°). If included, the dust disc is flared and follows the α prescription (see Shakura & Sunyaev 1973) to account for flux due to accretion. The reader is referred to Whitney et al. (2003a,b) for more details on the code and the possible geometries. We have used this code extensively and the reader is referred to de Wit et al. (2010, 2011) for more details on the application of the code to model MYSOs and their environments.

thumbnail Fig. 4

Panel a) logarithmically scaled 24.5 μm image of the standard W33A model. Panel b) the same image as panel a) but convolved with the COMICS instrumental PSF. Panel c) the COMICS image of W33A. Panel d) the azimuthally averaged radial intensity profile at 24.5 μm of W33A (points with error bars) alongside the model radial profile and that of a PSF standard (the lower line). The outflow axis of the model is aligned with the SE/NW orientation of W33A’s outflow.

Open with DEXTER

4.2. Modelling methodology

Our previous use of the Whitney et al. code consisted of constructing a series of dedicated models in order to reproduce the SEDs, interferometric and auxiliary spatial information of MYSOs. We will exploit these customised models appropriate for MYSO environments in our approach here, partially motivated by the lack of such models in the well-known SED grid (Robitaille et al. 2007). In particular, we will use a model geometry that provides a successful fit to the source W33A (de Wit et al. 2010). Its defining feature is that the MIR emission in the N-band originates in warm dust located in the walls of outflow cavities. The model proto-stellar envelope extends inwards down to the dust sublimation radius and it has low density outflow cavities as observed in several MIR observations of MYSOs (see e.g. De Buizer 2007). We note that the model does not explicitely include a dust disc (according to the prescription of the Whitney et al. code). This is motivated by the aforementioned MIR interferometric observations. Adding a disc generally results in high visibilities in the MIR whereas MYSOs typically exhibit low visibilities, indicative of a more extended dust distribution (see e.g. de Wit et al. 2010).

We illustrate the model 24.5 μm image for the W33A model in panel a) of Fig. 4. Also shown are the PSF convolved image (panel b), the COMICS image of W33A and the corresponding radial intensity profile, panels c) and d) respectively. The shape of the outflow cavities can clearly be delineated in the model image in panel a). This demonstrates that the dust located in the outflow cavity walls is warmed-up and emits a significant fraction of the MIR flux. W33A is know to drive an outflow in the SE/NW direction (see Davies et al. 2010, and references therein). The orientation of the outflow is traced by nebulosity in the NIR. While the 24.5 μm image of this object is resolved (panel c), it appears relatively symmetric. However, the image is slightly extended towards the SE, in the direction of the blue lobe of its outflow. In panel d) we compare the model’s radial intensity profile with the observed one and find they are essentially identical. We underline that we did not perform a model fit to the 24.5 μm image, we simply used the final W33A model presented in de Wit et al. (2010) to create the corresponding image.

Motivated by the good match of the W33A model and image and the knowledge that the SEDs of MYSOs are known to be similar, we proceeded to fit the observations of the remaining MYSOs in the following way. The SEDs of the MYSOs observed were constructed from the continuum data assembled on the RMS database (which makes use of the work of various authors) and the methodology of Mottram et al. (2011). The SED coverage differs from one source to the other. Generally one can distinguish between sources with continuum flux measurements up to 100 μm and those with sub-mm observations. We note that we do not fit explicitly the 9.7 μm silicate feature, contrary to the modelling efforts in DW09. This is simply because this information is not available for the majority of our sources.

To reproduce both the SEDs and 20 μm intensity distributions of the resolved sample, we used the W33A model geometry summarised above. We only altered parameters that can be expected to vary between objects: the inclination, the envelope infall rate (which scales the total dust mass in the envelope) and the opening angle of the outflow cavity (which may vary with time). We emphasise that we use one basic model to reproduce the observations. There may be differences between sources, for example the centrifugal radius and the properties of circumstellar discs – if any are present – may vary from source to source. The model will not account for such differences. However, it provides a simple test of the hypothesis that the bulk of the Q-band emission of MYSOs can be attributed to warm dust in outflow cavity walls. Here we outline the methodology followed in reproducing the observations.

For a given MYSO, the stellar luminosity was initially set to the value in Table 1 and the system inclination (the angle between the polar axis and the line of sight) was estimated from its MIR image. The stellar luminosity and infall rate were then varied until the model SED was consistent with that observed. This set the luminosity and provided first estimates of the other parameters. Once the SED was satisfactorily reproduced, the inclination and outflow opening angle were varied in an attempt to recreate the observed spatial intensity profile. This had an affect on the reproduction of the SED. Therefore, the infall rate was allowed to vary to enable a fit to both the SED and the radial profile. The model fitting was done by hand rather than using a χ2 minimisation technique as the computational time required to make a grid of model images for each object was prohibitively large.

We demonstrate the effect of varying the free parameters in Fig. 5. In general, the more inclined models appear more extended. This might be expected if much of the MIR emission originates in outflow cavity walls. At large inclinations, the projected area of the outflow cavity walls is larger and thus the resultant flux distribution is more extended. The cavity opening angle has a similar effect on the intensity distribution. This is because increasing the opening angle increases the distance at which the majority of the cavity wall is exposed to direct irradiation from the central object. This also increases the projected area of the outflow cavities and the extension seen in the azimuthally averaged intensity distribution. Since both parameters affect the observed extension, they are slightly degenerate. However, the two parameters have different effects on the gradient of the radial profile (see Fig. 5).

The final SEDs are presented in Fig. 6 and the associated azimuthally averaged intensity profiles are displayed in Fig. 7. The parameters of the selected models are presented in Table 2. Since the sample exhibits a varied morphology (see Fig. 3), we discuss the modelling results for each object individually in Appendix A. We discuss the general results in the following section.

thumbnail Fig. 5

An exploration of the parameter space immediately surrounding the model of G310.0135+00.3892. The plots show the effect of varying the free parameters (the inclination, opening angle and infall rate). For comparison, a PSF standard is also shown in the intensity distribution figures.

Open with DEXTER

Table 2

Key parameters of the individual models.

thumbnail Fig. 6

The SEDs of the resolved MYSOs compared to the final model SEDs.

Open with DEXTER

thumbnail Fig. 7

The azimuthally averaged intensity profiles of the resolved MYSOs (upper circles) alongside the model radial profiles (squares) and the associated PSF standard (lower circles). The lower limit to the y axis is set to the root-mean-square noise in the background of the MYSO images. The upper limit of the x is axis is set to the distance at which the MYSO profiles fall to the level of the background noise. The radial profile of G263.7759–00.4281 is also shown averaged over the upper and lower half of its bipolar morphology separately (short and long dashed lines respectively). The error bars represent the RMS within a given annuli and thus represent an upper limit on the uncertainty in the flux distribution.

Open with DEXTER

5. Modelling results

In general, the observed SEDs could be reproduced relatively well. There are several cases where the observed NIR fluxes could not be reproduced exactly. However, we do not weight the NIR fluxes heavily as these can be strongly dependent on dust opacities and the circumstellar geometry on small scales. In the majority of cases, barring objects associated with Hii regions, the spatial distribution of the 20 μm emission of the MYSOs observed can also be reproduced by our 2D axis-symmetric radiative transfer (RT) model. The individual model fitting results for the SEDs and the 20 μm intensity profiles are shown in Figs. 6 and 7. Clearly, because of the unequal SED coverage from source to source, the model fidelity differs. We do not attempt to provide a unique model for each single source but evaluate whether the W33A model is capable of providing satisfactory fits, by changing a few of the model’s most basic parameters. From this starting point, we find that in general, if the image exhibits a single source with an approximately symmetric morphology, both the object’s SED and the 20 μm intensity profile are well reproduced. By extension of the case for W33A, the majority of the 20 μm emission in these instances is very likely emission from warm dust in cavity walls.

When an object’s morphology appears either very complex or cometary, the model generally fails to reproduce the intensity profile. For example, this is the case for G343.5024−00.0145, which is a known ultracompact Hii region. For this object, a decent fit to the SED (panels I in Figs. 6 and 7) predicts a “MYSO” structure that is much more compact than the one observed, even with extreme values for the inclination and opening angle. A cometary appearance alone will often indicate an ionized nature for the region, and this is generally confirmed via radio continuum observations. That the model cannot reproduce the observations of MYSOs associated with Hii regions is perhaps to be expected as Hii regions represent a source of flux not considered by the model. Furthermore, the appearance of a Hii region likely signifies the end of the embedded accretion phase of a MYSO, and it is this phase the model is designed to represent.

We give an overview of the values for the varied parameters in Table 2. The opening angles of the best-fitting models are generally narrow (~10°). This implies that the outflows of the MYSOs observed are relatively well collimated. The envelope infall rates of the models, a parametrisation of the density distribution, are generally of the order 10-4   M   yr-1. The corresponding values of visual extinction, AV, are typically between 10 and 30, in agreement with the measured AVs of MYSOs (Porter et al. 1998), although several of the more inclined models have AVs in excess of 100. The total masses of the envelopes are generally of the order 10 000 M. We note that these are dependent on the envelope outer radii, which are set at 5    ×    105 AU. In the case of the W33A model, reducing the outer radius by a factor of 2 decreases the total mass by a factor of approximately 5. This is a result of adopting a low but constant density, ρ = 8 × 10-20   g   cm-3, within the outflow cavity. The envelope masses calculated are comparable to the clouds that the objects are located in (see Urquhart et al. 2011, and references therein). Even allowing for an uncertainty of a factor of a few, this indicates that these objects are still heavily embedded in their natal material and their environment still contains a large reservoir of material available to form further stars.

6. Discussion

6.1. On the MIR emission of MYSOs

Over the past decade, observational evidence indicating that the formation of massive stars, up to a mass of at least  ~30   M, is accompanied by phenomena characteristic of low mass SF has accumulated. Such phenomena include: circumstellar discs (see e.g. Kraus et al. 2010; Davies et al. 2010; Masqué et al. 2011; Goddi et al. 2011a), molecular outflows (see e.g. Beuther et al. 2002; Cyganowski et al. 2009) and jet activity (see e.g. Rodriguez & Bastian 1994; Guzmán et al. 2010).

Despite the many common characteristics of low and high mass SF, some observational phenomena seem to be exclusively associated with young, embedded high-mass stars. One example is the Class II 6.7-GHz methanol maser, whose activity is uniquely associated with massive star forming regions (e.g. Walsh et al. 1997). These sites evidently produce sufficient IR radiation in order to pump this transition, unlike sites of low mass SF. It was initially thought that this maser activity originated in circumstellar discs. As a result, a correlation between the extension of several MYSOs in the MIR and their maser emission led to the suggestion that the MIR emission of MYSOs traces discs (see e.g. De Buizer et al. 2000). However, observations with 8 m class telescopes were able to spatially resolve the MIR emission of a number of masing sources to reveal that their MIR emission is aligned with, rather than perpendicular to, their large CO molecular outflows (see e.g. De Buizer 2006, 2007). These findings associate the MIR emission of massive YSOs with outflow cavities, rather than circumstellar discs. Nevertheless, the kinematic diversity of methanol masers in a large sample of massive outflow sources does not support a uniquely identifiable source within the close MYSO environment (Cyganowski et al. 2009; Moscadelli et al. 2011). Therefore, the source of the MIR emission of MYSOs has yet to be conclusively established.

Recently, several MYSOs with compact MIR emission were investigated using long-baseline interferometry. In several cases, the spatially resolved N-band emission on scales of 100 AU was found to be consistent with thermal envelope emission (de Wit et al. 2007; Linz et al. 2009) and in particular with emission arising from the cavity walls (de Wit et al. 2010). However, in some other cases, the mas-scale MIR emission of MYSOs is suggested to be (at least partially) located in a circumstellar disc (see e.g. Follert et al. 2010; de Wit et al. 2011; Grellmann et al. 2011). In other words, there are hints that the N-band emission, resolved with interferometry, is not completely dominated by envelope emission. The analysis of the Q-band emission presented in this paper favours the cavity wall emission interpretation for a sample of MYSOs. Upon first inspection, this appears at odds with the scenario of a disc dominating the N-band emission. Here, we consider the possible contribution to the Q-band emission by a circumstellar disc in order to assess whether our finding agrees with theoretical expectations.

The recent study by Zhang & Tan (2011) discusses in detail the appearance of embedded young massive stars in the IR wavelength region. The authors use a dust radiative transfer code, as we do. However, their geometry differs from the one employed in this paper in that their code includes an accretion disc that extends from the stellar surface to the centrifugal radius. Therefore, their results allow us to assess the possible disc contribution to the Q-band emission of MYSOs. The RT calculations of Zhang & Tan (2011) reveal that outflow cavities have a large influence on the shape of the SED. Furthermore, in their calculation for a 60° inclination (see e.g. their Fig. 9), the outflow cavities are clearly significant features in the images, up to a wavelength of 70 μm. In the NIR, scattered light from the cavities makes a significant contribution to the near-IR emission. In the N-band, thermal emission from the disc and cental object are also evident. In the Q-band however, the emission from the base of the cavity clearly dominates the flux.

The relative contributions of the disc and envelope to the modelled SED depend critically on exactly those parameters with which we fit our observations, i.e. opening angle, mass infall rate and inclination. We note that the opening angle applied by Zhang & Tan (2011) is rather wide (2θ ≈ 90°), which would imply an object in a rather advanced phase of formation, following the expected widening of the opening angle as function of time (see e.g. Kuiper et al. 2010). Smaller opening angles increase the extinction towards the central regions and thus favour the dominance of cavity wall emission in the total MIR emission. A similar reasoning applies for the mass infall rate. It is therefore not surprising that the observed MYSOs, which are believed to be in a less evolved phase than the phase represented by the fiducial model of Zhang & Tan (2011), are dominated by the cavity wall emission. This argues against a large effect, if any, by a circumstellar disc on our VISIR images. This conclusion can also be inferred from the SED of an accretion disc, since it peaks at wavelengths much shorter than 20 μm. Furthermore, the envelope contains much more cool material than does any reasonable accretion disc, whose sizes for MYSOs are estimated to be less than 500 AU (see e.g. Patel et al. 2005; Hoare 2006; Reid et al. 2007; Jiménez-Serra et al. 2007; de Wit et al. 2011; Goddi et al. 2011b). We conclude that it is reasonable to expect that dust emission from the envelope (the cavity walls) dominates the continuum SED at wavelengths longward of the N-band. This is in agreement with our finding that the observed extension of the Q-band flux can be modelled as emission from outflow cavity walls.

6.2. The structure of MYSO envelopes: rotation versus outflows

We have shown that the morphology of the 20 μm emission of MYSOs is consistent with that predicted by models incorporating outflow cavities. We followed a similar methodology to DW09 and have improved on the modelling part of their analysis. DW09 compared their images to 1-D radiative transfer models. This was motivated by the rather spherical appearance of most sources, which suggests that their MIR emission originates from a relatively symmetrical dusty envelope. They found that, in general, their data were best fit by an envelope with a density gradient of ρ ∝ r-1. This density law is shallower than that obtained via a similar analysis performed in the sub-mm (see e.g. Mueller et al. 2002). It was therefore suggested that this flattening of the density law towards smaller spatial scales could be evidence for the onset of rotational support of the envelope at approximately 1000 AU. This suggestion is probably invalidated by the results presented in this paper based on more realistic RT models.

Thermal emission from an envelope with evacuated outflow cavities fits the SED and the N-band (de Wit et al. 2010) and Q-band (this paper) morphology of the luminous MYSO W33A. Additionally, in this particular case, the envelope model also fits the observed morphologies of the near-IR scattering nebula and the 350 μm emission. This provides strong evidence in favour of this model and our explanation of the 20 μm emission (warm dust in the cavity walls). Moreover, the W33A image (Fig. 4) is slightly extended towards the SE which is the direction of the blue lobe of the large scale outflow, reinforcing this interpretation. Importantly, we note that W33A is not the only case in our sample for which the extension at 20 μm is in the direction of the outflow position angle (see for example the sources G263.7759−00.4281 and G332.2941+02.2799). Therefore, our association of the observed MIR flux with outflow cavity walls can be substantiated.

The brightness of the cavity walls, as exemplified in panel a) of Fig. 4, is a consequence of the density contrast between the in-falling envelope and the outflow cavities. The radiation of the central source cannot penetrate far in the dense envelope, but can heat dust in the walls of the outflow cavity to large radii. Once the model images are convolved with an instrumental PSF, it is no longer obvious that the flux traces outflows, despite the source being resolved. The flux from the outflow cavity walls can appear more extended than symmetric envelope emission. A centrally concentrated density distribution (steep density law) in a spherical envelope produces a small thermal emission region in the MIR. The shallower the density law the larger the (normalised) emission region becomes. Therefore, the influence of the outflow cavities can mimic the effect of a shallow density gradient. We verified that the sizes of the emission regions differ strongly between a  −2 and  −1 density law, but the difference between a flat density law and a  −1 density law are relatively minor. Still the SEDs differ significantly. The steeper density law produces higher MIR fluxes. As an exercise, we fitted the standard 2D W33A model with a spherical model and find that the SED and intensity profiles are matched best by a density law between  −0 and  −1. This explains the results presented in DW09. Given the more realistic nature of the 2D models we have employed, we suggest that the preference for shallow density laws found by DW09 is the result of much of the MIR emission of MYSO emanating from warm dust in the cavity walls. An explanation in terms of rotational support on scales of  ~1000 AU cannot be substantiated.

6.3. SED fitting

Fitting SEDs has long been used as a tool to provide important constraints on the morphology of circumstellar environments (see e.g. Guertler et al. 1991). Still, the inclusion of spatially resolved observations at different wavelength regimes is indispensable for a correct interpretation of the observed emission (see e.g. van der Tak et al. 2000). This is because SED fitting can result in degenerate solutions. To assess the added value provided by the VISIR images presented in this paper, we use the results presented in Mottram et al. (2011). These authors used the grid of models provided by Robitaille et al. (2007) to fit the SEDs of a large percentage of MYSOs from the RMS survey in order to determine their bolometric luminosity. The grid of Robitaille et al. (2007) was created using the same RT code of Whitney et al. which we use in this paper. Therefore, we used the code to create images for the best fitting models found in Mottram et al. (2011) (in the filter used) and assessed whether the model images are consistent with our observations.

To illustrate the results of this exercise, we focus on the object G310.0135+00.3892. This is an object of particular interest as it is the object studied by Kraus et al. (2010) via a high angular resolution VLTI aperture synthesis image in the K-band. We generate 20 μm images of the set of models that best fit the object’s SED as found by Mottram et al. (2011). The result is shown in Fig. 8. The SED fitting tool preferentially returns models with large cavity opening angles (~30–50°), but these all fail to reproduce the observed intensity profile. The best fitting model of Kraus et al. (2010), which is also based on SED fitting with the model grid of Robitaille et al. (2007) and the code of Whitney et al., provides a better match. However, it still over-predicts the flux in the central 1′′ from the source (see Fig. 8). This is likely due to the large opening angle of the Kraus et al. model (40°) generating a MIR emission region which is too extended. The VISIR observations clearly demonstrate that, for the assumed intermediate inclination, a narrow opening angle is required. We conclude that the spatially resolved VLT-VISIR images at 20 μm provide valuable constraints on the nature of the envelopes surrounding young massive stars.

thumbnail Fig. 8

The azimuthally averaged intensity profile of G310.0135 (upper circles) compared to a PSF standard (lower circles) and the models that fit the SED of this object from Mottram et al. (2011, dashed lines). The solid line is the radial profile of the model of Kraus et al. (2010). The error bars represent the rms within a given annuli and thus represent an upper limit on the uncertainty in the flux distribution.

Open with DEXTER

7. Conclusions

In this paper we present diffraction limited, MIR imaging of a sample of MYSOs drawn from the RMS survey. By comparing the spatially resolved images of the MYSOs to 2D radiative transfer models, we constrain the structure of the envelopes that surround them. Here we list the salient results.

  • We spatially resolve the MIR (19.5 or 24.5 μm) emission of 14 MYSOs. The remainder of the MYSOs observed (6) are unresolved at the current resolution (~0.6′′).

  • The model of the MYSO W33A that we have developed previously can, in most cases, reproduce the images and SEDs of the MYSOs. This is not the case for objects associated with a Hii region. That the model can reproduce the infrared emission of a sample of MYSOs suggests that the circumstellar environments of MYSOs are relatively uniform. The large envelope masses of the models are consistent with the identification of MYSOs as objects in their accretion phase.

  • It is found that the extent of the MIR emission of the spatially resolved MYSOs can generally be attributed to warm dust located in the walls of outflow cavities. We suggest that the relatively shallow intensity profile gradients found by previous MIR diffraction limited imaging of MYSOs is indicative of outflow cavities (rather than alternative explanations such as rotational flattening).

  • We show that the varied morphology observed can be attributed to outflow cavities seen at a variety of inclinations. This emphasises that outflows are an ubiquitous feature of massive star formation.


2

Some of the distance estimates have since been revised to larger values.

Acknowledgments

H.E.W. acknowledges the financial support of the MPIfR. W.J.d.W. is grateful for the hospitality of G. Weigelt and the IR interferometry group at the MPIfR where this paper was finalised. J. Mottram is thanked for providing details of previous fits to the SEDs of many of the sources in this paper. This paper made use of information from the Red MSX Source survey database at www.ast.leeds.ac.uk/RMS which was constructed with support from the Science and Technology Facilities Council of the UK.

References

Online material

thumbnail Fig. 2

Unresolved sources. The images are scaled logarithmically. North is to the top of the page and east is to the left. The contours typically represent 1, 5, 25 and 75 percent of the peak flux.

Open with DEXTER

Appendix A: Notes on individual objects

Here we discuss the images and modelling results for each resolved MYSO.

A: G263.7759−00.4281

The 20 μm morphology of G263.775986−0100.4281 (IRAS 08448-4343) is reminiscent of the cavities of a bipolar outflow seen close to edge on. The direction of the elongation seen in the MIR image is consistent with that seen in the NIR via 2MASS images. The radio luminosity of the object is consistent with a jet rather than a Hii region (Hoare et al. 2007). This supports the notion that the bulk of the MIR emission arises in the walls of cavities evacuated by an outflow. Furthermore, the extension in the NW/SE direction is aligned with the H2 jet detected by Giannini et al. (2005), confirming that this object drives an outflow and the MIR morphology is associated with outflow activity.

The SED of this object can be reproduced with the selected model with an inclination of i = 87, which is consistent with the notion that this object is seen close to edge on. The intensity profile of the model is not entirely consistent with the data. However, the axis-symmetric model used cannot account for the slight asymmetry of the image. Indeed, the model can approximately recreate the extension of the southern part of the observed bipolar structure but not the northern half (see Fig. 7, panel A). The model cannot reproduce asymmetry at these wavelengths in an edge on configuration. Therefore, we surmise that the model reproduces the data as well as can be expected.

B: G265.1438+01.4548

This object exhibits a fairly symmetric morphology in the MIR. If the generic model featuring outflows is correct, this would suggest that this object is seen at a relatively low inclination. Modelling of this object’s CO bandhead emission is consistent with this scenario (Bik & Thi 2004; Wheelwright et al. 2010). The SED and intensity profile are relatively well reproduced by the generic model with an inclination of i = 32, which is consistent with the hypothesis that this object is viewed at a low inclination. The model slightly under-predicts the flux at distances greater than  ~2′′   from the centre of the intensity distribution. However, examination of the image reveals that this may not be related to the source as the core of the emission is concentrated within 2′′. Therefore, this slight discrepancy is neglected.

C: G268.3957−00.4842

G268.3957−00.4842 exhibits a relatively symmetric morphology in the MIR. As before, it is suggested that this is the result of a relatively low inclination. Indeed, the object’s SED and intensity profile are reproduced by a model with an inclination of i = 30, consistent with the hypothesis that this object is observed at a low inclination.

D: G269.1586−01.1383

The 20 μm image of G269.1586−01.1383 incorporates an additional, localised source of flux approximately 7′′   to the NE. Both sources exhibit an extended, slightly cometary morphology reminiscent of a Hii region. However, radio observations indicate that while the northern source is a Hii region (Urquhart et al. 2007), the southern source is likely to be a YSO. This is confirmed via low resolution NIR spectra. The spectrum of the southern source exhibits H2 line emission, suggesting that it drives an outflow.

The complex appearance of this object is difficult to recreate with the 2D, axis-symmetric codes of Whitney et al. In principle, a bipolar cometary morphology could be recreated by warm dust emission from the walls of cavities carved by a bipolar outflow with a large opening angle. However, the image contains no hint of a bipolar structure. Therefore, it might be expected that the model cannot reproduce the intensity profile of this source, as is found to be the case. We note that this object exhibits a notably different morphology to the rest of the sample, which may indicate it represents a different phase of massive star formation than the other MYSOs. This is partially substantiated by the fact that this object is associated with a Hii region.

E: G310.0135+00.3892

G310.0135+00.3892 (IRAS 13481-6124) was observed with the VLTI and AMBER by Kraus et al. (2010). These authors reconstructed images of this object in the K-band and detected a disc seen under moderate inclination. In addition, Kraus et al. (2010) discovered signs of an outflow orientated in the NE/SW direction, perpendicularly to the disc. The 20 μm morphology of this object is slightly extended in the NE/SW direction, i.e. along the direction of the outflow. This is consistent with the notion that the majority of the MIR emission traces warm dust in the walls of outflow cavities.

The observed SED and intensity profile are reproduced relatively well by a model with an inclination of i = 32, close to the intermediate inclination derived by Kraus et al. (2010). The model under-predicts the flux at distances greater than  ~1.5′′   from the peak of the intensity distribution. However, examination of the image reveals that at this distance the morphology is not symmetric. The axis-symmetric model cannot reproduce this morphology. Therefore, it is surmised that the model reproduces the data as well as possible and this slight discrepancy is neglected.

F: G318.0489+00.0854

The 20 μm image of G318.0489+00.0854 exhibits a slightly cometary morphology. This object is associated with two Hii regions detected via radio observations, both within several arcseconds (Urquhart et al. 2007). The near/far ambiguity over the distance to this source has yet to be resolved. As a result, it is impossible to unambiguously model its intensity profile and SED. Therefore, we do not attempt to reproduce the observations of this object.

G: G332.2941+02.2799

G332.2941+02.2799 (IRAS 16019-4903) exhibits a significantly extended morphology in the NE/SW direction. This is also seen in the NIR in 2MASS images and at 10 μm (see Mottram et al. 2007). This object is known to be associated with several bipolar outflows and the direction of the most prominent outflow is aligned with the extension seen in the infrared images (Henning et al. 2000). As a result, the observed morphology is consistent with the notion that the MIR flux largely arises from warm dust in outflow cavity walls.

Since these observations were taken, it has been determined that the luminosity of G332.2941+02.2799 is 844 L (see Mottram et al. 2011). This places the spectral type of this object at approximately B2/B3. Therefore, this object is an intermediate mass YSO rather than a massive YSO and should be excluded from the sample. Nonetheless, we attempt to recreate the observations with the generic model developed for the sample of MYSOs. The model cannot fully recreate the extension of the intensity profile observed, even with a high inclination and stellar luminosity. This may imply that the model developed for MYSOs is not directly applicable to intermediate mass YSOs. Alternatively, the strong outflow may cause shock heating in the cavities, which is not included in the model.

H: G332.986886−0100.4871

G332.986886−0100.4871 exhibits a fairly symmetric morphology. As before, it is suggested that this is the result of a relatively low inclination. The object’s SED is well reproduced by a model with an inclination of i = 15°, as is the object’s intensity profile. The model intensity profile deviates from that observed at a distance of  ~1.75′′ from the centre. However, this is where the object’s flux distribution falls to the level of the background and thus this discrepancy is not considered significant.

I: G339.622186−0100.1209

This object exhibits a slightly cometary morphology. As discussed in the case of G269.1586−01.1383, this is difficult to recreate at 20 μm with an axis-symmetric code. We note that this object is the most luminous in the sample with L = 5.2 × 105   L. Mottram et al. (2011) find that there is not a significant population of radio quiet MYSOs with such a luminosity. Sources with luminosities greater than L ~ 105   L are almost inevitably associated with a Hii region. Therefore, it might be expected that this object is associated with a Hii region. While this object is a radio non-detection (Urquhart et al. 2007), it is possible that a weak Hii region has escaped detection due to the large distance to this object (13 kpc). This could explain the cometary morphology of this source.

A relatively inclined model (i = 60°) reproduces the object’s SED and intensity profile relatively well. The mm flux is under predicted, but at the large distance to this object, several sources could have been included in the mm beam, resulting in an over-estimation of the flux. We note that, although the model recreates the observed intensity profile, the model image exhibits bi-polar structure while the observed morphology is appears more mono-polar. As discussed above, this source may harbour a Hii region, which could explain this discrepancy.

J: G343.502486−0100.0145

The image of G343.502486−0100.0145 exhibits a complex, notably extended morphology. Extended nebulosity towards the SE is also visible in 2MASS and Glimpse IRAC images. The extension of the 20 μm image is broadly consistent with that seen at

shorter wavelengths, although these data exhibit a more complex morphology. This source was detected at 4.8 GHz (Urquhart et al. 2007), and is thus classified as a Hii in the RMS catalogue. Therefore, it might be expected that the model developed to recreate the observations of radio quiet MYSOs cannot reproduce the complex morphology exhibited by this object. Indeed, the model that recreates the majority of the object’s SED fails to match the spatial extension of the 20 μm flux. We suggest that this is likely due to some of the observed flux originating in dust surrounding a Hii region, rather than the central source.

K: G343.521386−0100.5171

The image of G343.521386−0100.5171 reveals three localised sources. The central source clearly dominates the flux. However, one of the secondary sources is within 2′′ of the primary source. It is apparent that secondary flux will contaminate the intensity profile of the primary. Indeed, the image reveals that this source is less centrally concentrated than its isolated counterparts. A limited range of angles was used in constructing the azimuthally averaged intensity profile of this object in an attempt to limit contamination, but it is not clear how effective this was. The model fails to simultaneously reproduce the broad central maximum of the image and the observed SED. We suggest that this is the result of contamination of the intensity profile by flux from the neighbouring secondary source.

L: G345.0061+01.7944

G345.0061+01.7944 exhibits an extended morphology similar to that of G332.2941+02.2799. This suggests that the MIR image traces outflow structure seen close to edge-on. Indeed, the NIR spectrum of this object exhibits H2 emission indicative of outflow activity. Unlike the case of G332.2941+02.2799, a satisfactory fit to both the SED and the intensity profile was found. This is consistent with the hypothesis that the difficulties in reproducing the observations of G332.2941+02.2799 were the result of applying a model developed for high luminosity MYSOs to the case of an intermediate mass YSO.

M: G349.7215+00.1203

G349.7215+00.1203 is the second most luminous object in the sample with L = 3.1 × 105   L. The image of this object reveals a cometary shaped source of flux to the NW of the primary object. Radio observations confirm that the additional source of flux is a Hii region (Urquhart et al. 2007). The model reproduces most of the features of the SED and intensity profile of this object. The model under-predicts the flux at  ~100 μm. However, this could be due to the measured flux including a contribution from the Hii region. The model intensity profile exhibits a slightly different slope than that observed, although it is essentially within the uncertainty. We suggest that this could also be due to the presence of the Hii region.

All Tables

Table 1

Log of the observations.

Table 2

Key parameters of the individual models.

All Figures

thumbnail Fig. 1

The azimuthally averaged intensity profiles of the PSF standard stars. HD 124897 is the PSF standard taken with Subaru/COMICS. The solid line with square points is the average VISIR PSF.

Open with DEXTER
In the text
thumbnail Fig. 3

Resolved sources. The images are scaled logarithmically. North is to the top of the page and east is to the left. The contours typically represent 2, 5, 10, 25 and 75 percent of the peak flux.

Open with DEXTER
In the text
thumbnail Fig. 4

Panel a) logarithmically scaled 24.5 μm image of the standard W33A model. Panel b) the same image as panel a) but convolved with the COMICS instrumental PSF. Panel c) the COMICS image of W33A. Panel d) the azimuthally averaged radial intensity profile at 24.5 μm of W33A (points with error bars) alongside the model radial profile and that of a PSF standard (the lower line). The outflow axis of the model is aligned with the SE/NW orientation of W33A’s outflow.

Open with DEXTER
In the text
thumbnail Fig. 5

An exploration of the parameter space immediately surrounding the model of G310.0135+00.3892. The plots show the effect of varying the free parameters (the inclination, opening angle and infall rate). For comparison, a PSF standard is also shown in the intensity distribution figures.

Open with DEXTER
In the text
thumbnail Fig. 6

The SEDs of the resolved MYSOs compared to the final model SEDs.

Open with DEXTER
In the text
thumbnail Fig. 7

The azimuthally averaged intensity profiles of the resolved MYSOs (upper circles) alongside the model radial profiles (squares) and the associated PSF standard (lower circles). The lower limit to the y axis is set to the root-mean-square noise in the background of the MYSO images. The upper limit of the x is axis is set to the distance at which the MYSO profiles fall to the level of the background noise. The radial profile of G263.7759–00.4281 is also shown averaged over the upper and lower half of its bipolar morphology separately (short and long dashed lines respectively). The error bars represent the RMS within a given annuli and thus represent an upper limit on the uncertainty in the flux distribution.

Open with DEXTER
In the text
thumbnail Fig. 8

The azimuthally averaged intensity profile of G310.0135 (upper circles) compared to a PSF standard (lower circles) and the models that fit the SED of this object from Mottram et al. (2011, dashed lines). The solid line is the radial profile of the model of Kraus et al. (2010). The error bars represent the rms within a given annuli and thus represent an upper limit on the uncertainty in the flux distribution.

Open with DEXTER
In the text
thumbnail Fig. 2

Unresolved sources. The images are scaled logarithmically. North is to the top of the page and east is to the left. The contours typically represent 1, 5, 25 and 75 percent of the peak flux.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.