VLTI/MIDI 10 micron interferometry of the forming massive star W33A
aa r X i v : . [ a s t r o - ph ] N ov Draft version November 5, 2018
Preprint typeset using L A TEX style emulateapj v. 03/07/07
VLTI/MIDI 10 MICRON INTERFEROMETRY OF THE FORMING MASSIVE STAR W33A
W.J. de Wit , M.G. Hoare , R.D. Oudmaijer , J.C. Mottram Draft version November 5, 2018
ABSTRACTWe report on resolved interferometric observations with VLTI/MIDI of the massive young stellarobject (MYSO) W33A. The MIDI observations deliver spectrally dispersed visibilities with valuesbetween 0.03 and 0.06, for a baseline of 45m over the wavelength range 8-13 µ m. The visibilitiesindicate that W33A has a FWHM size of approximately 120 AU (0.030”) at 8 µ m which increases to240 AU at 13 µ m, scales previously unexplored among MYSOs. This observed trend is consistent withthe temperature falling off with distance. 1D dust radiative transfer models are simultaneously fit tothe visibility spectrum, the strong silicate feature and the shape of the mid infrared spectral energydistribution (SED). For any powerlaw density distribution, we find that the sizes (as implied by thevisibilities) and the stellar luminosity are incompatible. A reduction to a third of W33A’s previouslyadopted luminosity is required to match the visibilities; such a reduction is consistent with new highresolution 70 µ m data from Spitzer’s MIPSGAL survey. We obtain best fits for models with shallowdust density distributions of r − . and r − . and for increased optical depth in the silicate featureproduced by decreasing the ISM ratio of graphite to silicates and using optical grain properties byOssenkopf et al. (1992). Subject headings: stars: formation — stars: early type — techniques: interferometric INTRODUCTION
Massive young stellar objects are luminous, embed-ded infrared (IR) sources that show many signs thatthey are still actively accreting mass. Their luminos-ity (
L > L ⊙ ) is such that they are expected to beionizing their surroundings to produce an H II region,yet they only have weak radio emission due to ionizedwinds or jets (Hoare 2002). Most appear to be drivingbipolar molecular flows and (sub-)millimetre interferom-etry is beginning to reveal evidence of rotating, disc-likestructures on scales of hundreds of AU (Patel et al. 2005;Beltr´an et al. 2006; Torrelles et al. 2007). There is greatinterest in knowing the distribution of the infalling andoutflowing material on smaller scales to provide cluesto which physical processes are controlling the dynam-ics and setting the final mass of the star (Beuther et al.2007; Hoare et al. 2007).The mid-IR (8-13 µ m) emission from MYSOs is thoughtto arise in the warm ( ∼
300 K) dust in the envelope heatedby the young star. Previous modelling (e.g. Churchwellet al. 1990) has indicated that the size of the mid-IRemission region should be unresolved at the typical dis-tances of MYSOs by single-dishes and this is borne outby observations (Kraemer et al. 2001; Mottram et al.2007). Exceptions occur when MYSOs are viewed closeto edge-on when a dense torus can completely obscurethe bright central region and only the warm dust in theoutflow cavities is seen (de Buizer 2006).Here we present new results from mid-IR interferomet-ric observations of the MYSO W33A using the VLTIMIDI instrument. W33A has a kinematic distance of3.8 kpc and a luminosity as derived from IRAS fluxesof 1 × L ⊙ (Fa´undez et al. 2004). It has weak, com-pact, optically thick radio continuum emission (Rengara- School of Physics and Astronomy, University of Leeds,LS2 9JT, UK; [email protected], [email protected],[email protected], [email protected] jan & Ho 1996; van der Tak & Menten 2005) and broad( ∼
100 km s − ), single-peaked H I recombination emis-sion lines (Bunn et al. 1995) consistent with an ionizedstellar wind origin. IR images from 2MASS and Spitzer’sGLIMPSE survey clearly show a large scale monopolarnebula emerging to the SE, which is most likely to de-marcate the blue-shifted lobe of a bipolar outflow.The cold dust envelope has been studied at a few arc-second resolution at millimetre wavelengths (van der Taket al. 2000), probing the imprint of the star formationprocess on the circumstellar dust distribution. Many dif-ferent scenarios have been put forward predicting thisdistribution. In case of W33A, G¨urtler et al. (1991) mod-elled the SED with a constant density spherical envelopeand found an inner cavity radius of 135 AU (35 mas at3.8 kpc) with the dust close to the sublimation tempera-ture. MIDI provides data on comparable size scales of 50mas, and since the mid-IR size will be somewhat largerthan the sublimation zone we expect the emission to bewell-resolved. OBSERVATIONS AND DATA REDUCTION
MIDI is the VLTI’s two telescope beam combiner thatoperates in the thermal IR (Leinert et al. 2003). W33Awas observed with MIDI on 16 September 2005, employ-ing the UT2-UT3 telescope configuration for a projectedbaseline of 45 . ◦ east of north;this is perpendicular to the star’s outflow direction. Aninterferometric calibrator star HD 169916 was observedduring the same night. Lunar occultation and indirectmethods reveal the calibrator diameter to be a few mas(Nather 1972; Pasinetti-Fracassini 2001), too small to beresolved by MIDI. Here we present visibilities between8–13 µ m that are spectrally dispersed at a resolution ofR = 30. The observations were executed in the High-Sens
MIDI mode (for details see Przygodda et al. 2003;Chesneau et al. 2005; Leinert et al. 2004). We summa-rize the main elements of the procedure that produces the
Fig. 1.—
Left:
MIDI (lower full line) and ISO flux calibrated spectra of W33A.
Right:
The correlated flux determined with EWS (dotswith errorbars), and MIA (full line). Note the ragged behaviour of MIA at low flux levels. EWS errorbars correspond to the variations invisibilities as function of reduction parameters. EWS and MIA are consistent except at correlated flux levels below 0.1Jy (dashed line),these values must be considered upper limits (see Jaffe et al. 2004; Matsuura et al. 2006). raw dataset. The two beams are interfered producingtwo complementary interferometric channels that havea phase difference of π radians. The dispersed fringesare modulated according to an introduced optical pathdelay (OPD). Interferograms are recorded for a rangein OPDs, a few millimetre around optical path lengthequalization, called a fringe scan. For both W33A andHD 169916 a total of 8000 interferograms, correspondingto 200 scans were recorded. Subtracting the two interfer-ometric beams eliminates the background and enhancesthe fringe signal. The coherent flux was extracted usingtwo different MIDI software reduction packages: EWS(Jaffe 2004 ) and MIA (K¨ohler et al. 2005). MIA esti-mates the amplitudes in the power spectra of the fringesignal Fourier transform (incoherent estimation). Theamplitudes are proportional to the correlated flux. EWSon the other hand aims at adding the fringes to maximizethe signal to noise (coherent estimation). This methodcorrects first the fringe spectra for their corrupted phasecaused by atmospheric and instrumental effects. EWSuses the fringe scan as phase reference to estimate thegroup delay due to the atmosphere. The observationsshow that the atmospheric group delay varied within arange of 50 µ m over 100s, indicating relatively good at-mospheric circumstances during the night. The Fouriertransform of the group delay function reveals the typicalsawtooth phase change due to the introduced OPD, in-dicating that fringes have been measured. Removing theatmospheric and the (known) instrumental group delaysconstitutes a linear correction to the dispersed fringe sig-nal, and straightens the dispersed fringe spectra, i.e. thephase is independent of wavelength. Next, the phase off-set due to varying water refraction index between thetime of recording of the fringe spectra has to be ac-counted for. In principle all spectra can then be added toa final fringe spectrum. Final visibilities are obtained bytaking two spectra of the source (one with each telescope)immediately after the interferometric measurement. Theaccuracy of the final visibility spectrum is limited by thesky brightness variation between the interferometric andflux measurement, and can amount to 10-15%.The flux spectra are corrected for sky contribution us-ing a median sky subtraction. Spectra are extracted bysumming the counts within 3 σ of the mean position, de-termined from a Gaussian fit to each column along the wavelength axis. This procedure results in a W33A spec-trum with an S / N ∼
100 and S / N ∼
300 for the stan-dard star. Four spectra are recorded, 2 for each telescopebeam and combined via a geometric mean. This quantityis also what is obtained for the correlated flux after beamcombination and thus ensures consistency when derivingvisibilities. The absolute flux calibration is done usingHD 169916, which has an average flux density of 30Jybetween 8 µ m and 13 µ m (Cohen et al. 1999). The dif-ference in airmass between the observation of the fluxcalibrator and W33A is 0.05 leading to a negligible cor-rection on the final fluxes. Errorbars on the spectrumin Fig. 1 indicate the systematic difference in flux levelsbetween the two telescopes beams. Absolute flux cali-bration is uncertain upto at least 35% due to differencein flux level of HD 169916 in each telescope beam. RESULTS
Fig. 1 presents the MIDI flux spectrum and correlatedflux spectrum from EWS and MIA. The flux spectrumis compared to the ISO-SWS spectrum, taken from theISO Data Archive maintained by ESA. Both spectra aredominated by a very strong silicate feature, that con-tains solid-state ammonia and methanol absorption fea-tures, at 9 and 9.7 µ m respectively (see Gibb et al. 2000).At the central wavelength of the feature no flux wasrecorded, and the actual depth is unknown. The MIDIand ISO spectra portray a similar overall shape but theMIDI spectrum has a flux level that is about a factor 2less. In addition to the uncertainties due to flux calibra-tion, we ascribe this difference to the much larger ISObeam (80 ′′ ) in comparison to MIDI (2 ′′ ).The corresponding visibility spectrum is obtained bydividing the correlated flux by the flux spectrum and ispresented in left panel of Fig. 2. Visibilities are not plot-ted when corresponding to correlated fluxes smaller than0.1Jy, a value below which the measurements become un-reliable (Jaffe et al. 2004, Matsuura et al. 2006). Thevisibility spectrum shows that the silicate wings are notstrongly affected by the decrease in flux but instead fol-low the declining trend of the continuum visibilities. Ifwe would represent the emission by a Gaussian emittingdistribution then the FWHM size increases from 30 masat 8 µ m to 60 mas at 13 µ m. RADIATIVE TRANSFER MODELLING
IDI observations of W33A 3
Fig. 2.—
DUSTY model fits to visibilities and SED. Model is described by a r − . density law, and A V = 100, star parameters are L = 4 10 L ⊙ , T eff = 10 000K, and a distance of 3800 pc. Left:
Calibrated MIDI visibilities of W33A (dots with errorbars). The level ofsystematic uncertainty is indicated by the bar in the right top corner. Black squares connected by the dashed line correspond to the best-fitmodel.
Right:
SED of W33A traced by ISO-SWS and LWS spectra and IRAS points (open circles). Full line is the model fit. IRAS/ISOlarge beam effect at the longer wavelength is illustrated by the improved photometry using MIPSGAL images (filled circle).
Previous model fits of MYSO SEDs have indicated thatsimple spherical radiative transfer models with roughlyconstant densities best fit the near-IR and far-IR data(Churchwell et al. 1990; G¨urtler et al. 1991). How-ever, evidence exists that the outer envelopes (10,000 AUscales) have steeper density profiles (Hoare et al. 1990;van der Tak et al. 2000). We explore here whetherthe visibilities produced by material on 100 AU scales (ascale previously unexplored), and the SED can be simul-taneously matched by simple 1D spherically symmetricdust models. Arguably, such models are inadequate forMYSOs which are likely to consist of circumstellar disks,bipolar cavities, etc. Despite this, we explore whetherthe basic levels and trends of the dispersed visibilitiescan be matched by a single unresolved star deeply em-bedded in a dusty envelope before introducing more freeparameters.For this purpose, we employ DUSTY, a code thatsolves the scaled 1D dust radiative transfer problem (seeIvezic & Elitzur 1997). We used a spherically symmet-ric dust distribution illuminated by a central, unresolvedstar. The only non-scaled parameter entering the codeis the dust sublimation temperature. Gas emission isnot taken into account by DUSTY. Solutions are inde-pendent of the central source luminosity, and the SEDand visibilities can be scaled accordingly. The luminos-ity is the prime stellar parameter that sets the inner dustsublimation radius and thus the size scale; an increasecauses the size of the emitting region to be larger, as r sub ∝ √ L . The canonical parameters for W33A are L = 1 10 L ⊙ , 3.8 kpc (Fa´undez et al. 2004) and we ini-tially chose T eff = 25000 K typical for 15 M ⊙ star. Thesevalues make the underlying star larger than a main se-quence B1 star, consistent with the idea that a star un-dergoing accretion at a high rate is swollen (e.g. Behrend& Maeder 2001). We tuned the model output to the ob-servations by varying DUSTY’s input parameters, viz. radial density distribution, dust properties, and source T eff . The challenge is to construct a single model thatfits the observed visibilities, silicate wings, and mid-IRSED up to 100 µ m. We restrict the SED fitting to thiswavelength interval as it is well-known that 1D modelshave particular difficulty in reproducing the short wave-length region because of the sensitivity to the viewing angle (cf. Yorke & Sonnhalter 2002).We adopt a dust sublimation temperature of 1500 Kand an MRN-DL dust mixture (Mathis et al. 1977;Draine & Lee 1984) with a typical interstellar graphite-silicate composition. The outer bound of the model isset at 1000 times the dust sublimation radius, where thetemperature corresponds to the presumed ambient tem-perature of the ISM, between 10 and 25 K. The precisevalue for the outer radius does not affect our conclusionshere.Previous studies have shown that W33A is best de-scribed by an A V between 100 and 200 (Capps et al.1978; G¨urtler et al. 1991). A first result is that the MIDIvisibilities cannot be reproduced for any density distri-butions with power exponents between − . . A V between 100 and 200 and L = 1 10 L ⊙ . The modelvisibilities at 10 µ m are too small by an order of magni-tude. Decreasing A V is not a solution for normal dust,because of W33A’s exceptionally strong silicate absorp-tion feature. Average model sizes at 10 µ m of ∼ µ m MIPSGAL data that we have obtained from theSpitzer archive. These data are taken at a vastly superiorresolution compared to both IRAS and ISO and revealthe presence of at least three point sources and strongdiffuse emission within the IRAS beam. For W33A, wemeasure a 70 µ m flux of 1 . Jy (20% uncertainty), asshown in Fig. 2. This indicates that the true luminosityis significantly less than deduced from IRAS data.Even with the luminosity reduced to 4 10 L ⊙ , ther-mally supported cores with r − dependency (Larson1969) are incompatible with the observed trend of de-creasing visibilities with wavelength and the SED. In-compatible models are also found if we adopt a constantinfall velocity type distribution with r − . (Shu 1977).Reasonable fits to the SED that produce sizes compa-rable to the MIDI visibilities are found for r − . loga-tropic distributions (e.g. McLaughlin & Pudritz 1996),but again they do not reproduce the observed decreas-ing visibility trend with wavelength. DUSTY producessmaller sizes, and thus larger visibilities, for longer wave-lengths. This affects especially the red wing of the sili-cate feature, due to the diminishing opacity further outin the silicate wing. This mismatch in visibilities is how-ever removed if we reduce the slope to a much shallowerdensity distribution with a r − . law or a constant den-sity law. These distributions fit the short wavelengthregion of the SED increasingly worse because of the lossof warm dust. These shallower density distributions alsorequire the luminosity of the central object to be fainterthan sB9observed. We thus find that for standard ISMdust, the SED and the sizes at scales of 100 AU limit thedust to follow distributions between r − . or r − . pow-erlaws. The required optical depths are relatively high,yet the best models do not fit the SED particularly well.Different dust compositions influence both the magni-tude and shape of the visibility spectrum. We now ex-plore dust made of cold and warm silicates by Ossenkopfet al. (1992) and amorphous carbon. Ossenkopf sili-cates have the strong advantage of a larger optical depthin the silicate feature than the DL grains. This prop-erty is advantageous in case of W33A, because it allowsmodels with relatively moderate A V to produce deep sili-cate absorption. Models with moderate extinction fit theshorter wavelength fluxes better. We find that the typicalISM ratio of 0.88 between graphite and silicates does notproduce enough silicate absorption to fit W33A’s strongabsorption feature. A better correspondence is reachedif this ratio is reduced to 0.50 (see also Churchwell etal. 1990). The warm Ossenkopf silicates bring the extraadvantage of somewhat reducing the sublimation radius.Although the reduction in A V , and a better correspon-dence to the silicate feature depth produces reasonablecorrespondence between the shallow density models withthe observation, models that fit the data better are foundby lowering the T eff of the star. The change of this pa-rameter produces a decrease of the size of the inner dustboundary, because of the lower dust opacity with tem-perature. We finally arrive at a simultaneous fit (Fig. 2)to the red silicate wing, SED peak and visibilities for acentral object with a relatively low T eff of 10 K (corre-sponding to a B9 supergiant, again consistent with thenotion that an accreting star may be swollen) and a r − . density law (Fig. 2). The dust sublimation radius for T d = 1500 K is found at 7 mas (26.5 AU).In summary, fitting 1D DUSTY models to the visi- bilities and SED of W33A has shown that for nominalstellar parameters, size scales of the emission region aretoo large and produce the wrong trend with wavelength.Shallow radial density distributions produce this trendof size with wavelength and in order to fit the depth ofthe silicate feature as well, dust models with warm Os-senkopf silicates with an increased silicate over graphiteratio reproduce the observation best, provided the lumi-nosity and T eff of the central object are reduced. CONCLUSIONS
We presented mid-IR high-resolution dispersed in-terferometric observations of the forming massive starW33A. The visibility spectrum indicates an equivalentGaussian FWHM of the emitting region of ∼ µ m, increasing to ∼ µ m. We interpretedthe interferometric data with simple spherically symmet-ric DUSTY models representative of a (unresolved) starembedded in a dusty envelope, aiming for a simultane-ous fit to the SED, silicate profile and visibility spectrum.For any radial dust distribution we found that the canon-ical value of W33A’s luminosity is not compatible withthe visibilities. This is supported by MIPSGAL data.We found that even for a reduced luminosity, the modelproduces too large emitting regions causing the visibil-ities and SED to be mutually incompatible. Changingthe dust composition improves the situation for an in-creased graphite to silicate ratio and using warm silicateoptical constants from Ossenkopf et al. (1992). We re-sorted to a substantial lowering of the T eff in order toobtain a satisfactorily match between observables andmodels. Further coverage of the uv-plane will be veryrewarding, constraining more appropriate, higher dimen-sionality models, which will eventually lead in a properdescription of the circumstellar environment of an accret-ing massive star. A full 2D axi-symmetric treatment ofthis and other MIDI data will be the subject of a futurepaper.It is a pleasure to thank Olivier Chesneau and RainerK¨ohler for discussion on the MIDI data reduction.in order toobtain a satisfactorily match between observables andmodels. Further coverage of the uv-plane will be veryrewarding, constraining more appropriate, higher dimen-sionality models, which will eventually lead in a properdescription of the circumstellar environment of an accret-ing massive star. A full 2D axi-symmetric treatment ofthis and other MIDI data will be the subject of a futurepaper.It is a pleasure to thank Olivier Chesneau and RainerK¨ohler for discussion on the MIDI data reduction.