Table Of Content

DraftversionJanuary10,2011 PreprinttypesetusingLATEXstyleemulateapjv.08/22/09 AHIGHMASSDUSTYDISKCANDIDATE:THECASEOFIRAS18151-1208 C.Fallscheer1andH.Beuther MaxPlanckInstituteforAstronomy,Ko¨nigstuhl17,69117Heidelberg,Germany J.Sauter2andS.Wolf UniversityofKiel,InstituteofTheoreticalPhysicsandAstrophysics,Leibnizstrasse15,24098Kiel,Germany and 1 1 Q.Zhang 0 Harvard-SmithsonianCenterforAstrophysics,60GardenSt.,Cambridge,MA02138,USA 2 DraftversionJanuary10,2011 n ABSTRACT a Many questions remain regardingthe propertiesof disks around massive prototstars. Here we presentthe J observationsofahighmassprotostellarobjectincludinganelongateddustcontinuumstructureperpendicularto 6 theoutflow.SubmillimeterArray230GHzlineandcontinuumobservationsofthehighmassprotostellarobject IRAS18151-1208alongwithsingledishIRAM30mobservationsaffordushighspatialresolution(0.8′′)as ] A wellasrecoveryoftheextendedemissionthatgetsfilteredoutbytheinterferometer.Theobservationsof12CO confirmthe outflow directionto be in the southeast-northwestdirection, and the 1.3 mm continuumexhibits G anelongationinthedirectionperpendiculartotheoutflow.Wemodelthephysicalparametersoftheelongated . structurebysimultaneouslyfitting the observedspectralenergydistribution(SED)and thebrightnessprofile h p alongthemajoraxisusingthe3DRadiativeTransfercodeMC3D.Assumingadensityprofilesimilartothatof - alowmassdisk,wecanalsoreproducetheobservationsofthishighmassprotostellarobject. Thisisachieved o byusingthesamedensitydistributionandflaringparametersaswereusedinthelowmasscase,andscalingup r thesizeparametersthatsuccessfullymodeledthecircumstellardiskofseveralTTauristars. Wealsocalculate t s thataregionwithintheinner30AUofsuchahighmassdiskisstableundertheToomrecriterion. Whilewe a donotruleoutotherscenarios,weshowherethattheobservationsinthehighmassregimeareconsistentwith [ ascaledupversionofalowmassdisk. Implicationsonhighmassstarformationarediscussed. 1 Subjectheadings:stars: formation—stars: individual(IRAS18151-1208)—stars: earlytype v 6 1 1. INTRODUCTION suggested the presence and general morphology of outflow 3 intheregion. However,dataqualitywasnotgoodenoughto IncontrasttolowmassTTauristarsandevenintermediate 1 observespecificdetails. Narrow-bandH (λ=2.122µm)ob- mass Herbig stars of which numerousexamples exist where 2 . servations(Davisetal.2004)ofIRAS18151-1208showevi- 1 circumstellar disks have been detected (see review by Wat- denceoftwocollimatedmolecularoutflowsoriginatingfrom 0 son etal. 2007), the roleof disksin the case ofmassive star 1 formationisnotyetwelldefined. Inthe lowmasscase, itis two separate sources. In the case of low mass star forma- 1 clearthataccretiondisksandmolecularoutflowsareintegral tion, it is clear that accretion disks play a central role in the : partsoftheformationprocess. Whiletheoccurrenceofmas- existenceofoutflows. However,inthehighmassregime,ob- v servationsof scaled up versionsof these low mass accretion i sivemolecularoutflowsiswellestablishedincasesofmassive X star formation, evidence for the analogous disk component disksarerare(Cesaronietal.2007).Ourgoalistoshowthata scaled-upversionofadiskmodelingschemethatsuccessfully r has only tentatively been identified (see review by Cesaroni a etal.2007). reproducestheobservationsofdisksaroundTTauristarscan beappliedinthecaseofamassiveprotostellarobject,namely Inaconcertedefforttoshedsomelightonaccretiondisksin IRAS18151-1208. highmassstarformation,weobservedasampleofsevencan- TheT Tauristar IRAS04302+2247,nicknamedthe“But- didateswith the SubmillimeterArray(Zhangetal. in prep). terfly Star” because of its butterfly-shaped morphology, has ThehighmassprotostellarobjectIRAS18151-1208isoneof a well-defined circumstellar disk. Wolf et al. (2003, 2008) the sources included in this study. This source was also in- model the circumstellar environment of this source using a cludedinthesingledishsurveyofIRASsourcesconductedat MonteCarlo3D(MC3D)RadiativeTransfercode. Similarly, theIRAM30mtelescopebySridharanetal.(2002);Beuther thediskintheBokglobuleCB26(Sauteretal.2009)hasalso et al. (2002a), and has also been observed with the IRAM beenobservedandmodeled.CB26isalsoaTTauristardes- 30m and Mopra single dish telescopes by Marseille et al. tinedtobecomealowmassstarwithanedge-oncircumstellar (2008). IRAS18151-1208isatadistanceof3.0kpcandhas aluminosityofapproximately104L⊙(Sridharanetal.2002). disk. Additionally,ahandfulofotheredge-ondisksaroundT Tauristarshavebeenmodeledusingasimilartechnique[e.g. Previous single dish observations (Beuther et al. 2002b) HK Tau (Stapelfeldt et al. 1998), HV Tau (Stapelfeldt et al. Electronicaddress:[email protected] 2003), and IM Lupi (Pinte et al. 2008)]. On the theoretical 1Currentaffilliation:NationalResearchCouncil,HerzbergInstituteofAs- side, Whitney et al. (2003)use a similar Monte Carlo radia- trophysics,5071WestSaanichRd.,Victoria,BC,V9E2E7,Canada tive transfer approach to test various envelope/disk geome- 2Alsoat:MaxPlanckInstituteforAstronomy,Heidelberg,Germany triesforClass Isources. While theoreticaldescriptionssuch 2 Fallscheeretal. Forbothdatasets,thequasar3C273wasusedastheband- passcalibrator,andinbetweeneach10minutesourceobser- vation, 5 minute observations of the quasar 1743-038 were madeforphaseandamplitudecalibration. Thequasar1911- 201 was used as an additional phase and amplitude calibra- tor in the very extended data set. Fluxes for the calibra- tionsourceswereobtainedfromtheSubmillimeterCalibrator List4. Measured fluxes after the calibration were consistent withthe tabulatedvaluesto within 20%. A summaryof this informationcanbefoundinTable1. Forthemapping,a miriadrobustparameterof 0was used for the continuum while a more natural weighting scheme (largervalues of the robustparameter) was used for the line data. After combining the SMA data sets, the synthesized beam size of the continuum is 0.8′′×0.7′′. A summary of this and other observation parameters are given in Table 1. Spectralresolutionof0.55kms−1wasobtainedovera4GHz bandwidthdividedintoaloweranduppersidebandseparated by10GHz(seeFig.1)withcentralfrequenciesof219.6and 230.5GHzrespectively. ThelinesobservedbytheSMAare listedinTable2.ThesearediscussedinmoredetailinSection 3.2. 2.2. PicoVeleta30mTelescope The shortest baseline in our SMA datasets is 14 m which correspondstoaspatialscaleof19′′. TheSMAisnotsensi- tivetostructureslargerthanthis,sosupplementaryshortspac- inginformationwasobtainedintheon-the-flymodewiththe HERA receiver on the IRAM 30m telescope near Granada, Spain. Theseobservationsweremadeon2006November12 under very good observing conditions with τ(230) of 0.07. Thereceiverwastunedto230.5GHzcenteredonthe12CO(2– Fig.1.—ThespectrafromtheSMA.Upper:thelowersidebandspectrum 1)line. Aspectralresolutionof0.4kms−1wasobtainedwith atthecenterofthe1.3mmdustcontinuumpeakofthecombinedcompact andveryextendedconfigurationdata. Middle: thesameasaboveexceptfor a channel spacing of 0.3 MHz. The spectra were first pro- theuppersideband.Lower:theuppersidebandatapositionapproximately cessedwithCLASS90oftheGILDASpackage,thenfurther 7′′northeastofthecontinuumpeakonlyincludingthecompactconfiguration analyzedwithGREG.Whilethesingledishdatacaninprin- data. ciplebe combinedwiththe interferometerdata, theirquality is not good enough to do so in this case, so we analyze the as Whitney et al. (2003) and observational studies (Sicilia- datasetsindependently. Aguilar et al. 2006, for example) of low mass T Tauri disks arerelativelywelldeveloped,suchacompletepictureismiss- 3. OBSERVATIONALRESULTS inginthehighmassregime. 3.1. MillimeterContinuumEmission 2. OBSERVATIONS Combining the very extended and compact configuration 2.1. SubmillimeterArray interferometerdata,wemapthe1.3mmdustcontinuumata spatialresolutionof∼2500AU (Fig.2). Asidefromthepri- IRAS 18151-1208was observed at 1.3 mm with the Sub- marypeak,wedetectasecondarysource16′′tothesouthwest millimeter Array (SMA)3 in the compactand very extended of the primary. Zoominginto the primarypeak (Fig. 3) dis- configurationson2007July8and2007May17respectively. plays the elongationof the dust continuumin greater detail. Observations in the compact array were made under stable Thiselongationisperpendiculartothebipolarmolecularout- and excellent weather conditions, with a τ(230 GHz) of ap- flowshowninFig.4,andmaybeassociatedwithadiskcom- proximately 0.07. The very extended configuration obser- ponentinaroughlyedge-onorientation. Thepositionangles vations were made under observing conditionswith a τ(230 oftherelevantcomponentsarelistedinTable3. Tosimplify GHz) of 0.12. The combinationof these two configurations theanalysisofthe modeling,we haverotatedthezoomedin providedbaselinelengthsvaryingbetween14and520meters imageby61◦ suchthattheelongationliesalongthehorizon- correspondingto10kλand400kλat1.3mm. Thephaseref- talaxis. Thebrightnessprofilealongthisaxisisasymmetric, erence used was RA(J2000) = 18h17m57.1s and Dec(J2000) withthedifferenceinheightofthetwopeaksontheorderof =-12◦07′22′′withtheemissionpeakapproximately17′′east 1σ(seeFig.5). ofthesecoordinates. We adoptedarestvelocityv =32.8 LSR Wemeasureanintegratedfluxof0.4Jyintheregioncon- kms−1(Sridharanetal.2002). tainedwithinthe10,000AUx10,000AUboxshowninFig.3 andapeakfluxof0.042Jy. FollowingthemethodsofHilde- 3TheSubmillimeterArrayisajointprojectbetweentheSmithsonianAs- brand (1983) and Beuther et al. (2002b, 2005) we adopt a trophysicalObservatoryandtheAcademiaSinicaInstituteofAstronomyand Astrophysics,andisfundedbytheSmithsonianInstitutionandtheAcademia Sinica. 4http://sma1.sma.hawaii.edu/callist/callist.html IRAS18151-1208 3 TABLE1 SMAobservationparameters.Entriesincludethedateofobservation,thearrayconfiguration,thecalibrators,thesynthesizedbeamobtainedafter invertingandcleaningthedata,andthe1σcontinuumnoiselevel. Date Freq Configuration Calibrators beam rmscont [GHz] Bandpass Phase&Amplitude [′′] [mJy/bm] 2007May17 230 veryextended 3C273 1743-038&1911-201 2007Jul08 230 compact 3C273 1743-038 230 combined 0.8x0.7 2 Fig.2.—The1.3mmdustcontinuum. Thesizeofthebeamisindicatedinthelowerleftcorner. Contoursstartat3σincreasinginstepsof2σwhereσis2 mJy/beam. TABLE2 Observedlines Transition RestFrequency Eupper [GHz] [K] C18O2→1 219.560 16 SO65→54 219.949 35 13CO2→1 220.399 16 CH3OH8−1→70E 229.759 89 CO2→1 230.538 17 TABLE3 Positionangles(measuredEastofNorth)ofthecomponentsassociated withthesource. Component PA ◦ H2jet 128 SMA12COoutflow 130 30m12COoutflow 108 continuumelongation 29 dustopacityindexβ of 2.35(Weingartner& Draine2001) op which corresponds to the dust model used in our modeling describedinSection4. Usingareferencewavelengthof250 µm, this opacity index corresponds to an opacity κ of 0.19 cm2g−1 at 1.3 mm. Assuming a dust temperature of 30 K, Fig.3.—Azoomedinviewoftheprimarypeakofthe1.3mmdustcontin- and an optically thin system, our measured flux at 1.3 mm uum.Herewehaverotatedtheimageby61◦inordertoaligntheelongation withthehorizontalaxis. Thesizeofthebeamisindicatedinthelowerleft correspondsto a mass of 220 M⊙ and a beam averaged col- corner. Theintegratedfluxinthisareais0.4Jy.Onearcsecondcorresponds umndensityof1.1×1025cm−2. Themassandcolumndensity to3000AU. would of course be lower if we use different values for β , op but to maintain consistency in the comparison with T Tauri stars, we use β of 2.35 for the modeling. Using the Inter- op 4 Fallscheeretal. stellar Mediumvalueforβ of2.0(κ=0.35cm2g−1)(Hilde- surefluxes. Weusethesingledishfluxestocalculatetheout- op brand1983),themassandcolumndensitythenwouldbe120 flow mass, dynamical age and outflow energetics following M⊙ and6.4×1024 cm−2. Withavalueforβop of1.4(κ=0.93 the approachof Cabrit & Bertout (1990). Using an average cm2g−1)(Ossenkopf&Henning1994),themassandcolumn spatial extentof 30′′ for each 12CO outflow lobe, and maxi- density become 45 M⊙ and 2.4 ×1024 cm−2. These masses mumblueshiftedandredshiftedvelocitiesof10and17kms−1 are, of course, affected by missing flux in the interferometer respectively, we derive the properties listed in Table 4. Our data.Comparingourinterferometric1.3mmfluxtothesingle derivedvaluesmatchthoseofBeutheretal.(2002b)reason- dishpeakfluxatthiswavelengthwithinthe11′′ beamofthe ablywell,andarelikelybetterestimatesduetopoorweather MAMBO bolometer obtained by Beuther et al. (2002a), we conditionsatthetimeoftheearlierobservations. Theseval- estimate thatless than60%ofthe fluxwithin this11′′ beam uesconfirmtheoriginalclaimthattheregionwilleventually isfilteredout. formahigh-massstar. 3.2. LineData 3.3. SpectralEnergyDistribution Toward the continuum peak, there are only four spectral Forthepurposeofmodeling,wegathereddataforthespec- lines present in the 4 GHz bandwidth of the SMA data (see tral energy distribution (SED) of IRAS 18151-1208. The Fig. 1). Based solely on chemicalevolution, it is difficultto wavelengthsandassociatedfluxesthatweincludeinourSED determinewhetherthissourcehasnotyetreachedthehotcore (Fig. 6) are tabulated in Table 5. The 1.3 mm data pointfor phase, or whether it has already passed through this evolu- theSEDisdeterminedbymeasuringtheinterferometricflux tionaryphase. Threeofthefourlinespresenttowardthecon- ofIRAS18151-1208withinthecentral10,000AUby10,000 tinuum peak, 12CO(2–1), 13CO(2–1), and SO(6–5), are not AU region. We note that the flux measured with the inter- highdensitytracers,andC18O(2–1)appearstobeaffectedby ferometer is a factor of four smaller than the corresponding the outflow(see Sect. 5.3). Hence, with the data athandwe integratedsingledishfluxreportedbyMarseilleetal.(2008). are unable to infer the kinematics of the elongated structure Thediscrepancybetweenafactor4hereandafactorofabout detected in the continuum and discussed below. Methanol, 2 discussed at the end of Sect. 3.1 comes from the fact that a molecule linked to rotation in Fallscheer et al. (2009), is Marseilleetal.(2008)integrateoveramuchlargerareathan presentafew(∼7)arcsecondstothenortheastoftheprimary just the 11” beam of the MAMBO observations. Since we peakbutisextremelyweakdirectlyatthelocationofthepri- wouldalsoliketousethe450and800µmfluxesreportedin marysource. Marseilleetal.(2008),wealsodividethesefluxesbyafactor 4 since they were measured over the same spatial scales as 3.2.1. Outflow the1.2mmdata. Wejustifythisapproachbythefactthatthe SED approximately follows a power law in this wavelength PrevioussingledishobservationsbyBeutheretal.(2002b) regime. weretakenduringnon-idealweatherconditions,sothequality For the SED, the 24.8 µm point comes from Campbell ofourcurrentdatasetissignificantlyimproved.Nevertheless, the 12CO (2–1) emission is consistent with the outflow sug- etal. (2008)and the12.8µm and18.7µm data pointscome fromVISIRdata(ESOproposal075.C-0454,PIFuller). The gestedbythisolderdataset. TheSMAdatainFig.4demon- VISIRdataweretakenon5August2006withtheVLT.The stratethatourobservationsfollowthemorphologyoftheH 2 observationsweremadewith theQ2and[Neii]filterswhich jetsseenbyDavisetal.(2004). Inbothsetsofobservations, have central wavelengths of 18.72 µm and 12.81 µm and the two outflows observed are highly collimated. The two widthsof0.88µmand0.21µmrespectively.Integrationtime outflowswe detect are roughlyperpendicularto one another with the Q2 filter was 6 minutes and 5.5 minutes with the andappeartoemanatefromtwoseparatesourcesdetectedin [Neii] filter. The standard star HD169916was used for flux our 1.3 mm dust continuum. Outflows associated with low calibration(Cohenetal.1999). massstarformationdependonthepresenceofadisk. While Thenominalpeakpositionsvaryby∼1′′ betweenthemil- massive outflowsare likely also gas entrainedby jets driven limeterandmid-infrareddata.However,inthemid-infrared,a bymagnetohydrodynamicaccelerationfromadisk,adescrip- smallfieldofviewdoesnotallowhighpositionalaccuracyso tionofothermethodsisgiveninArceetal.(2007). weattributetheemissiontothesamesource.Alternativelythe Several instances of anomalous behavior are apparent in offsetmaybeduetoarealdifferenceinspatialpositionwhich Fig. 4. First, in the southeast-northwest outflow associated wouldleadtoaloweringoftheSEDatthesewavelengths.We with the primary continuum source, there is a blueshifted notethispossibility,butproceedassumingthatthedetections component on the northwest side which is otherwise asso- atdifferentwavelengthsarespatiallycoincident. ciated with redshifted emission (approximately located at ∆R.A.=4′′, ∆Dec=7′′). This might be an indication that the ThedataofCampbelletal.(2008)exhibitasilicateabsorp- tion dip at 10 µm. This dip is likely caused by absorption outflow hasa small inclinationangle relative to the plane of oftheenvelope. Sincethemillimeterinterferometerdatafil- thesky suchthatwe observebothrecedingandapproaching terouttheenvelope,we areprimarilyinterestedin thesmall sectionsoftheoutflow.Next,itisinterestingtonotethatboth scale disk-like structure, so we do not try to reproduce this sidesoftheoutflowassociatedwiththesecondarycontinuum feature. source are dominated by blueshifted emission. This region is dominated by redshifted emission in the single dish data (Fig.4)indicatingthatitislargescaleemissionandisfiltered 4. MODELING outintheinterferometricdata. Afterthe3Dradiativetransfermodelingofthecircumstel- AsdiscussedinSect.2.2,wedonotcombinethesingledish larenvironmentinthelowmassregimewasprovensuccess- datawiththeinterferometerdata.Instead,weusetheinterfer- ful (Wolf et al. 2003; Sauter et al. 2009), we wanted to test ometerdatatoaccuratelydeterminetheoutflowmorphology whetherasimilarmethodcouldworkinthehighmassregime. as described above, but rely on the single dish data to mea- We hypothesized that a scaled up version of a T Tauri disk IRAS18151-1208 5 Fig.4.—Theoutflowin12CO(2–1).IntegratedmapsofIRAS18151-1208overtheentirevelocityrangeoftheline.Blueshiftedemission(thicksolidcontours) integratedovervelocitiesof23-29kms−1andredshiftedemission(thickdashedcontours)integratedovervelocitiesof40-50kms−1.Left:SMA12CO(2–1)and continuum(grayscalewiththinsolidcontours)overlaidontheH2mapofDavisetal.(2004).ContoursforCOstartat3σandincreaseinstepsof3σwhereσis 0.15and0.16Jybm−1·kms−1fortheblue-andredshiftedemissionrespectively.Right:IRAM30m12CO(2–1).The1.3mmdustcontinuumpeakisindicatedby thetriangle,andtheoutflowdirectionandextentasdeterminedbytheSMAdataareindicatedbyarrows.TheregionplottedintheSMAfigureatleftisindicated bythedottedlines.Contoursstartat3σandincreaseinstepsof3σwhereσis11.4Jybm−1·kms−1fortheblueshiftedcomponentand10.4Jybm−1·kms−1for theredshiftedcomponent. TABLE4 OutflowpropertiesderivedfromCOobservations.TotaloutflowmassMt,momentump,energyE,size,outflowdynamicalaget,outflowrateM˙out, mechanicalforceFmandmechanicalluminosityLmaregiven.Inputsarediscussedinthetext. Mt[M⊙] p[M⊙km/s] E[erg] size[pc] t[yr] M˙out[M⊙/yr] Fm[M⊙km/s/yr] Lm[L⊙] 12 190 3.1×1046 0.42 32000 3.8×10−4 6.0×10−3 8.0 the11′′ beamof980mJy(Beutherpriv.comm.) issmoothly TABLE5 SourcesofSEDdatapoints. distributed, the single-dish flux per SMA beam is then 980 mJy/[11′′x11′′/(0.8′′x0.7′′)]∼3mJywhichisonly1.5times (µλm) fl(Juyx)a (+Jy/-) ape(′r′t)ure reference thermsof the SMA continuumimage (see Tab.1). Further- more, we can estimate the extinction at mid-infrared wave- 12.8 26.3 2 5.0 H.Linz,priv.comm. 18.7 41.0 2 5.0 H.Linz,priv.comm. lengthsduetotheenvelope. BasedonanextinctionAV of72 24.8 101 30 1.9 Campbelletal.(2008) mag(Campbelletal.2008),assumingR =5.0weestimate V 450 12.4 3.7 Marseilleetal.(2008) the extinction at 12 µm to be approximately 2 mag (Mathis 800 1.9 0.6 Marseilleetal.(2008) 1990). Therefore, although a componentis missing, model- 1300 0.4 0.1 3.3 thesedata aThe450and800µmfluxesareafactorof4smallerthanthesingledish ingonlythediskappearstobeareasonableapproximationfor fluxespublishedinMarseilleetal.(2008)inordertomakeamoredirectcom- comparisonpurposeswithlow-massdiskmodels. parisonwiththeinterferometerdata.FurtherdetailsaredescribedinSect.3.3. 4.1. TheMonteCarlo3DRadiativeTransferCode maybeabletoreproducetheobservationsofwhatmaypos- siblybeadiskaroundamassiveprotostellarobject. Usingtheself-consistent3DMonteCarloRadiativeTrans- While we are aware that the disk in IRAS 18151-1208is fer code, MC3D, developed by Wolf et al. (1999); Wolf still embeddedin a large-scaleenvelope,we preferto model (2003),wemodeltheelongationseeninthedustcontinuumof onlythediskstructure.Theinterferometricmapclearlyshows IRAS18151-1208.Inoursetup,MC3Dcalculatesthedensity anelongatedstructure,andwhileitmaystillhavesomecon- distributioninasphericalregionsurroundingthecentralstar tributionfromaflattenedenvelope,applyingadiskmodelto dividedinto91equallyspacedangulardivisionsand50log- this structure appears to be a reasonable approach. Further- arithmically spaced radial divisions such that the successive more,whilewehavemissingfluxonscalesoftheSMApri- grid cell is ∼1% longer than its inner neighbor. This distri- marybeam(55′′)ontheorderof75%,itisonlyontheorder bution is chosen to better resolve the density gradientin the of 60% when going to scales of the 11′′ beam of the IRAM moremassiveregionclosertothecentralstar. Using100log- 30mobservations(Sect.3.1and3.3). Goingtoscales ofthe arithmically distributed wavelengths between 0.05 and 2000 SMAsynthesizedbeam(0.8′′x0.7′′),themissingfluxiseven µm,MC3Dcomputesthethermalstructureofthedisk,derives less severe. Assuming that the single-dish peak flux within SEDs, and produces images at any of the 100 wavelengths. 6 Fallscheeretal. Forthepurposeofthiswork,weuseall100wavelengthsfor reasonably describes the observations of the high mass pro- the temperature structure and the SED, and we compute an tostellarobjectinquestion,wewouldexpectthevaluesforα imageat1.3mm. andβtoremainthesameasinthelowmasscase. Although thevalueforh betweenthetwocasesisnotdirectlycompa- 100 4.2. TheDisk rablebecause100AUrepresentsverydifferentregionsofthe respectivedisks,wecanscaleh toaproportionallysimilar We use a density profile of a parameterized rotationally 100 fractionofthedisk’ssizeandexpectthattobeproportionally symmetric disk to produce fits to the observed SED and larger than in the low mass case. In the ideal case, the disk brightness profile along the elongation. To maintain a more massshouldscalerelativetothevolumeofthesystemwhich even comparison, we assume the density profile used in the canbevariedviatheouterradius. lowmasscases(Wolfetal.2003;Sauteretal.2009)although Theparametersthatweincludeinthemodelingare: we cannotexcludeother density profileswith this approach. Thedensitydistributionweuseisdescribedbythefollowing • α, β, and h : the density distribution exponent, the 100 equation: disk flaring exponent, and the scale height of the disk at 100 AU respectively. These quantities describe the ρ (r )=ρ R∗ αexp −1 z 2 (1) dustdensitydistributionviaEquations1and2. Wetest disk cyl 0 r !  2 P !  13valuesofαbetween1and4and11valuesofβbe- wherezistheheightaboveocyrlbelowthemidhplane,R∗ isthe htwoeldenin0.t5heancdas3e.oTfhloewremlataisosndαisk=s3((Wβ-o21lf) eistsahl.ow20n0t3o; radius of the central star and rcyl is the radial distance from Sauteretal.2009). Weoriginallyvariedαandβinde- the centralrotation(z-)axis. The normalizationconstant, ρ0, pendently. However, after discovering that modifying isdeterminedbytheconditionsplacedonthetotalmassand αalonedidnotcausesignificantchangesinthemodel- volumeofthesystem,and ingresults,subsequentmodelsfollowedtheα=3(β-1) 2 rcyl β relation.Wetest8valuesofh100between1and30AU. P =h (2) h 100(cid:18)100AU(cid:19) • θ: the inclination of the disk with respectto the plane where h is the scale height of the disk 100 AU from the of the sky. The observations indicate the presence of 100 centralstar.Theshapeofthemodeleddiskisdescribedbythe aflattenedobject,andthiswouldonlybevisibleatin- parametersαandβwhicharetheradialdensityparameterand clination angles close to edge-on. However, we test 6 thediskflaringparameterrespectively. Ifαandβsatisfythe inclinationanglesrangingfrom45◦to90◦(edge-on). rtoelathtiaotnoαf=Sh3a(kβu-r21a),&theSnutnhyeaedven(s1i9ty73d)is.triTbhuetiodnepceonrrdeesnpcoendosf • Rin: the inner radius. We initially chose small values forR of10,20,and50AU,butalsotriedsignificantly thesurfacedensityonthediskshapeparameterscanthenbe in larger values such as 400, 1000, 2500, and 4000 AU obtained by integratingEquation 1 with respect to z: Σ(r ) cyl beforesettlingon20AUinourbestfitmodel. ∼rp wherep=β-α. cyl • R : theradiusatwhichthedensitygoestozero. The out 4.3. TheDustModel outerradiuscanroughlybeconstrainedbythe1.3mm continuum observations to a nominal value of 5,000 Under the assumption that the gas is optically thin in the AU,soweassumethisvaluethroughout. wavelength regime we are working in, it is the dust that is describedbythedensitydistributioninEq.1. Thedustgrains • M :thedustmass.UsingthedustmodelofWeingart- disk are assumed to be spherical, and we adopta power law size ner&Draine(2001)wedeterminearoughestimatefor distributionfollowingtheMRN(Mathisetal.1977)approach themasswithinthe5,000AUradiustobe220M⊙(See wherethenumberofdustgrainsofa givenradiusa isgiven detailsinSect.3.1).Fromthis,wederiveadustmassof by 2.2M⊙ basedonagas-to-dustratioof100. Takingthis n(a)∼a−3.5. (3) asourstartingestimate,wesubsequentlydeterminethe massbyensuringthattheSEDfitgoesthroughthe1.3 Dustgrainsizesavaryfrom5nmto250nmwhicharetypical mmdatapoint. oftheinterstellarmedium. We tested two different chemical compositions. The first • Central star: due to the luminosity constraint of followsthe graphiteand smoothedsilicate mixture of Wein- 20,000L⊙ describedinthenextsection(Sect.4.5),the gartner & Draine (2001) with relative abundances of 37.5% most massive star we consider is a B0.5[T=26300 K, (graphite)and62.5%(astronomicalsilicate). Theothercon- L=20000 L⊙, R=7.0 R⊙, M=15 M⊙]. We also test a tainsbymass24.2%ironandiron-poorsilicates,5.6%troilite, B1[T=25400 K, L=16000 L⊙, R=6.4 R⊙, M=13 M⊙] 25.6%refractoryorganics,4.4%volatileorganics,and40.3% and a B2 [T=22500K, L=6000L⊙, R=5.6 R⊙, M=10 water(Pollacketal.1994). M⊙]centralstar. 4.4. TheParameterSpace 4.5. ObservationalConstraints A motivating question behind this study is whether a low FortheSED,thefittingisconstrainedinastraightforward massdiskcanbescaleduptodescribeahighmassdisk.This mannerbytheobservedfluxes.WeattempttofitourSEDover wouldmanifestitselfintheformofhavingsimilarparameters the entire wavelength regime for which we have data points describingtheshapeofthediskandscaledupparametersfor available. thequantitiessuchasdiskmassandextent,accompaniedby WhileweusetheSEDtoguideoureffortssincetheeffects a more massive central star. In the scenario that scaling-up ofmodifyingmodelparametersaremorenoticeablethere(see IRAS18151-1208 7 Figs.7and8)thanintheimage,weaimtoreproducetheelon- gatedstructureseeninthe1.3mmdustcontinuumperpendic- ulartotheoutfloworientation.Amongotherinterpretationsof whatmaybetakingplaceintheregion,possibleexplanations includeadoublesourceandascaleduplowmassdisk. Here wefocusonwhetherascaleduplowmassdiskisconsistent withtheseobservations.Althoughwewereabletoreproduce asymmetricdoublepeakedstructurepriortosmoothingwith the PSF by making the inner radius larger, the convolution with the observed beam size smoothed out this effect in the imagebrightnessprofile. Wethenshiftedourfocustofitting theshapeoftheelongatedstructure’souteredges. Thisisnot asignificantoversightbecausethedifferenceintheheightof thetwopeaksisonlyontheorderof1σ. Itisworthmention- ingthatusinganextremelylargeinnerradiusflattensoutthe peakofthebrightnessprofilealongthemajoraxis(seeFig.8), butbecausewedonotknowtheoriginofthisslightasymme- try,wedonotincludeanythinginthemodelingtoaccountfor Fig.5.—Thedensityprofilealongtheelongationaxisofthebestfitmodel. thespecificasymmetrynorthedoublepeakandthereforedo Thesolidlineisthecutthroughthemidplaneoftheobservedelongatedstruc- notgiveitanyextraweightinourdeterminationofthebestfit ture. Thedottedlinesare+/-1σ,andthedashedlineisthemodel. Thedata pointsusedintheweightedleastsquaresfitareincludedintheplot. model. Forthemodelparameters,wewereabletoplaceroughcon- straintsonthecentralstar,theouterradius,andthemassofthe systembeforemodelingbegan. Sridharanetal.(2002)place theluminosityoftheIRAS18151-1208regionat20,000L⊙. We therefore use this as a limiting luminosity for the central star in our modeling. We use an outer radius of 5,000 AU which is the radius at which the midplane density goes to zero and correspondsapproximatelyto the 4σ contour of the elongation in the 1.3 mm dust continuum. Making the assumption that the region is optically thin, we can obtain a startingestimateforthemassofthesystemfollowingtheap- proach described in Section 3.1. Our measured flux at 1.3 mmcorrespondstoamassof220M⊙. Wetakethisasastart- ingestimate,butsincetherearemanysourcesofuncertainty in thiscalculationsuchasmissing fluxandunknownoptical depthtonameafew,wedidnotstrictlylimitthemassinour tests to thisvalue. Forthe test caseswith a gasmass of 220 M⊙, slightdifferencesontheorderofafewpercentbetween theobservedfluxandthecalculatedfluxat1.3mmsometimes Fig.6.—TheSEDofthebestfitmodel.Duetodifferencesintheobservation comeabout.Weattributethistothefactthatdespiteusingthe methods,the450and800microndatapointshavebeenadjustedforamore same dust model, MC3D takes optical depth into considera- directcomparisonwiththe1.3mmdatapoint(seeSect.3.3fordetails.) tionwhilethemassestimatedoesnot. bycalculatingthedifferencebetweenthetwocurvesatthesix 4.6. ComparingtheModelwithObservations pointsindicatedinFig.5. Weightingeachpointintheimage andSEDequally,wecombinethesedifferenceswiththedif- TheparametersofthebestfitmodelaresummarizedinTa- ferencebetweenthesixobservedandmodeledSEDpointsto ble 6 and for comparison,the best fit modelof the Butterfly arriveatthe finalvalueusedto determinethe bestfit model. Star (Wolf etal. 2003,2008)andCB 26(Sauteret al.2009) Thenumericalvalueoftheweightedleastsquarescalculation are also included in this table. The correspondingSED and cannotbe used on its own, but comparedto the valuesfrom brightness profile of the best fit model are shown in Figs. 5 othermodels,wehaveanumericalvalueonwhichtobaseour and 6. All data points in the SED as well as the six points decisionofthebestfit. Ourbestfithadachisquarevalueof showninFig.5wereusedtodeterminethebestfitmodel. 2.4. While manymodelshadchisquarevaluesless than10, Wedeterminethebestfitmodelbasedonconsideringboth somewereashighas30. theSEDandthebrightnessprofileofthe1.3mmimage. For each test we calculate a reduced chi-squared value as de- 5. DISCUSSION scribed in the last paragraph of the next section (Sect. 5.1). 5.1. TestingtheParameterSpace ThemodeledimageisconvolvedwiththeobservedPSFand thebrightnessprofilealongthemajoraxisisthenobtainedby As described in Sect. 4.4, a wide range of modelparame- summingthe flux fromthe pixelsin the columnsperpendic- tersweretestedbeforesettlingontheonesusedinthebestfit ular to the elongation. This is shown as the dashed line in model.Figures7and8containoverlaidSEDsandbrightness Figs.5and8.Thepixelsfromcolumnsalignedperpendicular profilesfromalongtheelongationthatdemonstratetheeffects totheelongationintheobservedimagearealsosummedand ofvaryingthemodelparameters.TheeffectsontheSEDsare thisisplottedasthesolidlineinFigs.5and8. Acomparison conspicuous,butthedifferencesinthebrightnessprofilesare betweentheobservedandmodeledbrightnesscurvesismade moresubtle. Wethereforequantifyeachfitbydetermininga 8 Fallscheeretal. Fig.7.—TheeffectontheSEDofvaryingindividualparameters. Althoughnumerouscombinationsofparametersweretested(seeSect.4.4foradescription oftheparameterspace),thesedemonstratehowindividualparametersaffecttheSED.Thesolidlineisalwaysthebestfitmodel(seeFig.6).SEDsvarying(from lefttoright,toptobottom)theflaringexponent(β),thescaleheightat100AU(h100),theinnerradius(Rin),thecentralstar’sspectraltype,theinclinationangle (θ),andthedustmass(Mdisk). IRAS18151-1208 9 TABLE6 Theparametersusedinthebestfitmodel. Object α β h0 Mdust θ Rin Rout R∗ T∗ L∗ (AU) (M⊙) (◦) (AU) (AU) (R⊙) (T⊙) (L⊙) IRAS18151a 2.4 1.3 200 2.2 60 20 5000 6.4 25400 16000 IRAS04302b 2.37 1.29 15 0.0007 90 0.07 300 2 4000 0.92 CB26c 2.2 1.4 6 0.003 85 45 200 2 4000 0.92 (Pθa=ra9m0◦etceorrsreinspclounddestthoeanraeddiagle-doennsdiitsyk)d,isthtreibinuntieornapnadroamuteetrerradαi,iflofartihnegdpeanrsaitmyedteisrtrβib,ustciaolneRhienigahntdoRfouthteanddistkheartad13iRuosu,ttehm0,pderuasttumrea,sasndMlduimski,nionsciltiynaotfiothneacnegnlteraθl star. aThiswork bWolfetal.(2003) cSauteretal.(2009) 10 Fallscheeretal. reducedchi-squaredvalue for all tests performed. Although havior. Smaller values of beta (β < 1.0) mean that wedonothaveenoughinformationabouttheenvelopecom- the disk is self-shadowedand thus less surface area is ponent to make quantitative statements about its effects, we availableforabsorptionandreemissionofradiation. A discussherehowitsomissionaffectsourmodelingresults. moreflareddisk(largervalueofβ)obscuresthecentral AsFigs.7and8demonstrate,someparameterscanbecon- starmorebecauseofitslargersurfaceareawhichinter- strainedbetterthanothers.Wefixtheouterradius,R based ceptsmoreofthestar’sradiation. Bothoftheseeffects out on the observationsand although we still vary the values of leadtoareductionofmid-infraredfluxintheSED.For the parameters for the disk’s dust mass m and the lumi- h , thicker disks (larger values) hide the central star dust 100 nosity, temperatureand radiusof the centralstar, we do this whilesmallervalueshavethesameamountofmassdis- mainlytoseehowvaryingtheseparametersaffectsthemodel. tributed over a smaller volume leading to higher den- Ultimately, the valuesfor these parameterswere determined sity whichincreasestheopticaldepth. Asaresult,the by the 1.3 mm observations rather than through the model- temperaturegradientbecomessteeperleavingasmaller ingprocess. By notincludinganenvelopecomponentin the areaforreradiationwhichinturnmeanslessfluxinthe model,allofthemillimeterdustemissionisattributedtothe micronregime. disk. This means that the disk dust mass is likely overesti- mated. The outer radius as defined by the millimeter image • The density distribution exponent α cannot be con- mayalsobeoverestimatediftheactualdiskcomponentisbe- strainedbythedatathemselves. Higherresolutiondata ing obscuredby an envelope. The inner radius, on the other would be necessary to providestronger constraints on hand,isnotwellconstrainedbyourmodels,andtheaddition α. After constrainingβ, we use the direct relation for ofanenvelopecomponentwouldnotaffectthesituation. Shakura&Sunyaev(1973)disksbetweenαandβ[α= Theflaringexponentβisconstrainedrelativelywellacross 3(β-1)]todetermineavalueof2.4forα.Thisapproach 2 the entire infrared and submillimeter wavelength regime as satisfactorilyreproducesourobservations. well as by the brightness profile, and we choose the corre- spondingvalueforthedensitydistributionexponentαbased • Inclination angle: the difference in the SEDs is small ontherelationknownforlowmassdisks(Shakura&Sunyaev forinclinationanglesbetween45◦ and60◦,butweuse 1973). Becausesomeofthediskmaybeobscuredbytheen- 60◦inourbestfitmodelbasedontheobservationaldata velope componentand especially the less dense outer flared thatindicateamoreedge-onorientation. TheSEDde- regions,ourvaluefortheamountofflaring,characterizedby viates further from the observations at inclination an- β,maybehigherthanifanenvelopehadbeenincludedinthe gles larger than 60◦. At very high inclination angles modeling.Becauseoftheβ-αrelation,αwouldbesimilarly (e.g.80◦),theupperlayersofthediskblockthecentral affected.However,theactualnatureoftheoverestimationde- star from view and the SED drops significantly in the pendsonthedensitydistributionintheenvelope,e.g.itsradial micronregime. power-law. The parameters for the disk scale height h , and the in- • Disk mass: inthemassregimewe tested, theonlypa- 100 clination angle θ are predominantly constrained only by the rameter that significantly affects the millimeter wave- SEDinthemid-infraredwavelengthregime.However,thein- lengthsideoftheSEDwithoutaffectingthemicronside clination angle is further limited by the millimeter image as isthedustmass. Figure7makesitapparentthatusing it is seen as a flattened elongated object and hence is closer thedustmodelofWeingartner&Draine(2001),adust to edge-onopposed to face-on. As the envelopecomponent massof2.2M⊙ isnecessarytofittheobserved1.3mm plays more of a role at these shorter wavelengths, the val- flux. Withagas-to-dustratioof uesderivedfortheseparametersshouldbetreatedwithsome M caution. Inthepresenceofathicksphericallysymmetricen- Gas ∼100, (4) velope, the elongation of the disk would still be observed, MDust providingsupportforthe limitationon the inclinationangle, but other envelope geometries may have different effects on thiscorrespondstoatotaldiskmassof220M⊙. Thisis muchlargerthanthemainsequencemassofaB1cen- theobservedcomponent. Usingsimilarreasoningtothedisk shape parameter β, the disk scale height would also be an tral star (13 M⊙) indicatingthat this system cannotbe Keplerianandshouldbequiteunstable.Toquantifythis overestimate because of the envelope component obscuring statement, we calculate the stability in the system un- theouterpartofthedisk. dertheassumptionthatitisadisk. FollowingToomre Varying the dust composition as described in Section 4.3 (1964),theconditionforstabilityis didnotsignificantlyaffectourresults,sowepresenthereonly theresultsoftheWeingartner&Draine(2001)compositionto c (r,z)κ (r) maintainamoredirectcomparisonwiththeButterflyStarand Q(r,z)= s k >1 (5) πGΣ(r) CB26. where c is the local sound speed, κ the angular fre- s k • The disk shape parameters β (flaring parameter), and quency of the disk derived from Kepler’s Law, G the h (scaleheightat100AU)stronglyaffecttheshorter gravitational constant and Σ the surface density. Fig- 100 wavelengths in the SED. As demonstrated in Fig. 7, ure9indicatesthatregionsonlywithinthecentral-most both lower and higher values of β and of h than 30AUofthesystemarestableundertheToomrecrite- 100 used in the best fit model significantly decrease the rion. These resultsare consistentwith whatKratter& flux in the shorter wavelength (micron) regime of the Matzner(2006);Vaidyaetal.(2009)findintheirmod- SED. Although the SED may be affected by a com- eling. Asexpectedforastar-disksystemthathasasig- plex interplay between several effects, it is likely that nificant fraction of the mass in the disk, the modeled the optical depth is primarily responsible for this be- diskinthissituationisunstableandislikelyfragment-