ALMA Observations of the massive molecular outflow G331.512-0.103 II: physical properties, kinematics, and geometry modeling
Carlos Hervías-Caimapo, Manuel Merello, Leonardo Bronfman, Lars Åke-Nyman, Guido Garay, Nadia Lo, Neal J. Evans II, Cristian López-Calderón, Edgar Mendoza
DDraft version December 27, 2018
Preprint typeset using L A TEX style emulateapj v. 05/12/14
ALMA OBSERVATIONS OF THE MASSIVE MOLECULAR OUTFLOW G331.512-0.103 II: PHYSICALPROPERTIES, KINEMATICS, AND GEOMETRY MODELING
Carlos Herv´ıas-Caimapo , (cid:63) , Manuel Merello , Leonardo Bronfman , Lars ˚Ake-Nyman , Guido Garay , NadiaLo , Neal J. Evans II , , Cristian L´opez-Calder´on , and Edgar Mendoza Departamento de Astronom´ıa, Universidad de Chile, Casilla 36-D, Santiago, Chile Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M139PL, UK Universidade de S˜ao Paulo, IAG Rua do Mat˜ao, 1226, Cidade Universit´aria, 05508-090, S˜ao Paulo, Brazil Joint ALMA Observatory (JAO), Alonso de C´ordova 3107, Vitacura, Santiago, Chile Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712-1205, USA and Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon, 34055, Republic of Korea
Draft version December 27, 2018
ABSTRACTWe present observations and analysis of the massive molecular outflow G331.512-0.103, obtainedwith ALMA band 7, continuing the work from Merello et al. (2013a). Several lines were identifiedin the observed bandwidth, consisting of two groups: lines with narrow profiles, tracing the emissionfrom the core ambient medium; and lines with broad velocity wings, tracing the outflow and shockedgas emission. The physical and chemical conditions, such as density, temperature, and fractionalabundances are calculated. The ambient medium, or core, has a mean density of ∼ × cm − and a temperature of ∼
70 K. The SiO and SO emission trace the very dense and hot part of theshocked outflow, with values of n H ∼ cm − and T ∼ −
200 K. The interpretation of themolecular emission suggests an expanding cavity geometry powered by stellar winds from a new-bornUCHII region, alongside a massive and high-velocity molecular outflow. This scenario, along withthe estimated physical conditions, is modeled using the 3D geometry radiative transfer code
MOLLIE for the SiO(J= 8 −
7) molecular line. The main features of the outflow and the expanding shell arereproduced by the model.
Subject headings:
ISM: clouds — ISM: molecules — ISM: jets and outflows — stars: formation INTRODUCTIONThe formation of massive stars is an important topic instellar astrophysics that is open for debate. Despite thewealth of knowledge that has been obtained on youngmassive stellar sources, still there is no universal agree-ment on how massive stars are formed and evolve (Garay& Lizano 1999; Zinnecker & Yorke 2007; Kennicutt &Evans 2012).There is an accepted paradigm on how low-mass starsform and evolve (Shu et al. 1987). However, problemssuch as scarcity of young sources, large heliocentric dis-tances (several kpc in some cases) and the lack of ap-propriate spatial resolution (that has been addressedwith interferometers such as ALMA just in recent years)make the subject of formation of massive stars an openquestion. Two main theories compete on trying to ex-plain this mechanism: the monolithic gravitational col-lapse/turbulent core accretion model (McKee & Tan2003), which is basically the extension of the model offormation of low-mass stars to the massive ones, but withhigher accretion rates and energetics; and the competi-tive accretion model (Bonnell et al. 2001, 2004), whichstates that low-mass seeds that will eventually form mas-sive stars compete for the available gas in the gravita-tional potential of their native cluster.Since massive stars have very short lifetimes, sourcesthat are currently in the process of formation are veryvaluable. Molecular outflows are commonly detected to- (cid:63)
E-mail: [email protected] ward protostellar objects, and observed physical prop-erties such as outflow power, force and mass loss rate,seems to correlate over a large range of luminosities, sug-gesting that outflows found in massive star forming re-gions are a scaled-up version of those found in their low-mass star counterparts (Tan et al. 2014). The injectionof momentum and energy from outflows are importantin star forming regions at small and large spatial scales,although the role of feedback at different stages of proto-stellar evolution is still not clear (Frank et al. 2014; Bally2016). Thus, from an observational point of view, a prin-cipal objective is to identify and characterize examples ofyoung, massive and powerful outflow-associated sources,as is the case for the work here.The present work focuses on arcsecond-resolutionALMA observations of the massive molecular outflowG331.512-0.103 (Bronfman et al. 2008). This source islocated in the tangent of the Norma spiral arm, at aheliocentric distance of ∼ . mass of ∼ × M (cid:12) ), which shows evidenceof ongoing massive star formation (Merello et al. 2013b,which we will refer to as MM13b). There are several char-acteristics that make this object valuable and unique:the presence of very broad emission wings in CO, CSand SiO, which indicate the presence of a very powerfuloutflow; the compact emission (not resolved at 8 . (cid:48)(cid:48) a r X i v : . [ a s t r o - ph . GA ] D ec Herv´ıas-Caimapo et al.line of sight; an expanding bubble geometry, which indi-cates the presence of stellar winds possibly arising fromthe exciting star within the hyper compact HII region(Merello et al. 2013a, which we will refer to as MM13a);and high energetics, which makes this source one of themost powerful massive outflows discovered. For example,in the compilation of outflows in O-type young stellarobjects by L´opez-Sepulcre et al. (2009), only five objectshave higher bolometric luminosity.Part of the ALMA observations and the main re-sults are summarized in MM13a. The SiO(J= 8 − −
2) lines show broad velocity wings( ±
70 kms − ). The SiO line shows ring-like emission, sug-gesting an expanding motion, probably a cavity beingblown-out by the powerful stellar output from the centralsource. The H CO + (J= 4 −
3) line traces the systemicvelocity structure, as well as the surrounding core withits narrow emission.This work is the follow-up of MM13a. Here we investi-gate for the first time 18 newly analyzed lines, along withthe 4 already analyzed there, with the goal of character-izing the physical conditions, kinematics and morphol-ogy of the massive molecular outflow. In this new work,we analyze the SiO and SO molecules, widely recog-nized as outflow tracers (Schilke et al. 1997b); CH CCH,which thermalizes at densities of ∼ cm − and is avery good estimator of temperature (Fontani et al. 2002;Molinari et al. 2016); and SO isotopologues, which aretracing the hot-core chemistry (Charnley 1997; Wakelamet al. 2004). The paper is organized as follows: Section2 describes the ALMA band 7 observations, giving themain parameters of the observed spectral lines. Section3 presents the results of the observations, including in-tegrated emission maps and position-velocity plots. Sec-tion 4 presents the analysis of physical conditions, ge-ometry, and kinematics performed on the source. A 3Dradiative transfer model of the SiO(J= 8 −
7) line is alsoincluded. Our conclusions are summarized in Section 5. OBSERVATIONSThe observations were performed with the AtacamaLarge Millimeter/submillimeter Array (ALMA) duringCycle 0, as described in MM13a. The primary beam was17 . (cid:48)(cid:48) . (cid:48)(cid:48) × . (cid:48)(cid:48)
68, witha position angle of − ◦ .
6. The interferometric obser-vations miss the recovery of large spatial scales above ∼ − . (cid:48)(cid:48) (cid:48)(cid:48) × (cid:48)(cid:48) centered at α =16 h m . s , δ = − ◦ (cid:48) . (cid:48)(cid:48)
4. The identified linesare marked across both bands. The image shows severalblended lines that required a careful channel-by-channeldetermination of the emission mask during the clean-ing reduction process. The rich spectra exhibit sulphur-bearing, carbon chains and other complex molecules, characteristic of hot core line emission found towardother well-studied sources at similar frequency bands,such as Cepheus A East (e.g., Brogan et al. 2007, 2008),Orion KL (Schilke et al. 1997a), and G5.89-0.39 (e.g.,Hunter et al. 2008).For the present work, we focus on 18 (new) + 4(fromMM13a) lines, including, among others, the hot-corechemistry/shock tracer SO , a K-ladder of four CH CCHlines (commonly used as a source temperature test),the optically-thin line H CN probing densities up to10 cm − , and HC N, which is found toward warm anddense regions such as hot cores, where it is shieldedagainst destruction by photo-dissociation and C + ions(Rodriguez-Franco et al. 1998; Prasad & Huntress 1980).Figure 1 shows the integrated spectrum of each of theselines as a function of velocity (the systemic velocity ofthe source is −
89 kms − ). A couple of these lines ap-pear blended or contaminated with emission at similarfrequencies. We leave the analysis of the rest of the linesobserved in the spectra, mostly associated with complexorganic molecules such as dimethyl ether (CH OCH )and ethyl cyanide (C H CN), for upcoming studies.The H CO + (4 −
3) line shows a secondary peak ofemission at ∼ −
50 kms − . MM13a reports this as a“molecular bullet”, or gas expelled at very high velocityby the powerful energetics of the outflow. Their pres-ence has been reported previously in star formation re-gion outflows (Tafalla & Bachiller 2011). However, oncewe analyzed the full band and the other lines containedin it, this might be no longer the case. A line with anexcited vibrational level, HC N( J = 38 − v = 1)at 346.9491233 GHz, matches with this secondary peak.Other star forming regions, observed at similar frequen-cies, show the presence of this line (Brogan et al. 2008;Nagy et al. 2015). We also observed, in our band-7 data,the HC N( J = 38 −
37) line at 345.60901 GHz. Compar-ison between them in velocity and morphology might in-dicate that both belong to emission of the same molecule,so the presence of a molecular bullet can be discarded,instead appearing to be another line that was not recog-nized in the first analysis of the ALMA data. RESULTS3.1.
Moment 0 maps
Figure 2 shows imaging of the 0th moment mapsof emission, integrated in the velocity range − . − . − (systemic velocity range).In the integrated spectra, from Figure 1, we can iden-tify two groups of lines. The first group (labeled wide lines) are the ones with broad velocity wings. They arethe SiO, S O, HCO + , H CN, HC N, and most of theSO lines. At the systemic velocity, most of these lines,especially the SiO emission, trace an emission with aring-like shape (see Figure 2).The second group (labeled narrow lines) includes linesthat have narrow emission at systemic velocity. They arethe CH CCH, CH OH, and H CO + lines .Spatially, the emission from all lines comes roughlyfrom the same central region, shown in Figure 2, with adiameter of ∼ . (cid:48)(cid:48)
0. Most of the lines trace the ring-like H CO + also shows evidence of high-velocity wings, but witha signal-to-noise ratio (cid:46) LMA Observations of G331.512-0.103 II 3 OH(4 E ) 0246 SO (5 )051015 CH CCH(21 ) 0.02.55.07.5 SO (8 )04812 CH CCH(21 ) 0369 SO (15 )0369 CH CCH(21 ) 2.55.07.5 SO (15 )048 CH CCH(21 ) 051015 SO (20 )100 90 802345 SO (5 ) 100 90 800.01.53.04.5 SO (23 ) Narrow spectral lines
Velocity [kms ] F l u x [ J y ] CN(4 3) 0246 S O(9(9) 8(8))081624 HC N(38 37) 0306090 SO (19 )0153045 SO (7 ) 0255075 SO (20 )0153045 SO (8 ) 081624 SO (25 )0153045 SO (11 ) 0.01.53.0 SO (26 )100 500306090 SO (16 ) 100 500.00.61.21.8 SO ( v =1,19 ) Wide spectral lines
Velocity [kms ] F l u x [ J y ] Fig. 1.—
Integrated spectra of molecular lines discussed in this work and not presented in MM13a. The left group corresponds to lineswith narrow velocity emission, while the right group corresponds to lines with wide velocity emission. The dashed vertical lines shows thesystemic velocity of − . − . emission that is clearly seen at the systemic velocity.The 862 µ m continuum emission is shown in Figure 2as the inverted color map in the bottom right corner. Itis composed of a single peak that almost coincides withthe center of the SiO emission ring, and weaker emissionthat follows the edge of the ring feature.3.2. Position-velocity (PV) plots
To study the gas kinematics traced by the line emis-sion, we made PV plots of most of the observed linesalong two axes: parallel and perpendicular to the out-flow axis with a position angle of 102 ◦ .
5, as defined inMM13a. A “slit” of 11 pixels ( ∼ . (cid:48)(cid:48)
5) is used. Bothkinds of PV plots are shown in Figure 3. ANALYSISWe interpret the wide group of lines as tracers ofshocked high-velocity hot gas. The narrow group of linestrace the ambient core cold gas emission at the systemicvelocity. As reported in MM13a, we interpret the kine-matics seen in the PV plots of SiO as an expanding shell.We present our interpretation of the source, shown inFigure 4. The bipolar outflow is outlined by emissionfrom shock tracers, such as SiO and SO . We considerthe presence of an expanding shell, which we interpretas shocked gas material being blown out by the stellarwind from the ultra compact HII region at the center.The ring-like emission we see in the maps of the lineswould correspond to this expanding shell. The narrowgroup of lines traces emission at systemic velocities (e.g.CH CCH) and the ambient dense and warm gas core.4.1.
Physical conditions
Estimation of column density and temperature
SiO analysis -
Assuming a common T rot for the SiOtransitions and that they are optically thin, we calculatethe column densities of the outflow wings. See, for ex-ample, the analysis in the similar outflow source W51 North (Zapata et al. 2009). In this ALMA data set, weonly have observed one line of SiO at high resolution.Therefore, for the excitation temperature value, we con-sider the estimate from MM13b, using only two transi-tions of SiO, obtained with APEX, for each outflow wing:123 ±
15 K (blue wing) and 138 ±
29 K (red wing).The estimated averaged column densities are listed inTable 1. The significant values are the ones calculatedin the blue and red peak, 1 . × cm − and 1 . × cm − , respectively, which are taken as the columndensities of the outflow wings. Estimating the size ofthe emission from the 50% peak contour level in the 0thmoment maps in the blue and red peaks, and using anabundance of X SiO = 1 . × − (see Section 4.1.2), weestimated the masses to be ∼ M (cid:12) , for each outflowwing. This value is about half of the mass calculated forthe wings in Bronfman et al. (2008), with CO. SO analysis - One of the principal molecules observedin the ALMA data set corresponds to SO , since sev-eral transitions fall in the observed bandwidth. In thiscase, assuming optically thin emission, we use the rota-tional diagram technique (Linke et al. 1979; Goldsmith &Langer 1999) to obtain its column density and rotationaltemperature.The partition function for SO is (Claude et al. 2000) Q ( T rot ) = 5 . × (cid:115) ( T rot / K) ABC/
MHz , (1)where A =60778.5511 MHz, B =10317.96567 MHz and C =8799.80750 MHz.Since some SO lines are blended, a compromise mustbe reached between the range of the wings to be usedand the number of lines. In the ideal case, one wantsto include the complete width of the wings and use asmany lines as available. Using a narrow range of veloci-ties is not advisable since the column density is underes-timated. Two SO lines (J=5 , − , and J= 7 , − , )were discarded since they are blended. Another line Herv´ıas-Caimapo et al. SiO(8-7)
HCO + (4-3) H CN(4-3) S O H CO + (4-3) CH OH(4 E ) HC N(38-37) CH CCH(21 ) SO (26
9, 17
8, 20 ) SO (16
4, 12
3, 13 ) SO (8
4, 4
3, 5 ) (5 ) SO (23
4, 20
3, 21 ) SO (8
4, 4
3, 5 ) Fig. 2.—
Integrated emission maps at systemic velocity of some representative lines in this work. The peak integrated intensity is shownin each panel in units of Jybeam − . The bottom right panel shows the 862 µ m continuum emission map. The (0 . (cid:48)(cid:48) . (cid:48)(cid:48)
0) offset correspondsto 16 h m s . − ◦ (cid:48) . (cid:48)(cid:48) LMA Observations of G331.512-0.103 II 5 . S i O ( ) . H C O + ( ) . H C N ( ) . H C N ( ) . S O ( , ) . S O ( , , ) . S O ( ) . H C O + ( ) . C H O H ( E ) . C H CC H ( ) . S O ( , , ) . S O ( , , ) . S O ( , , ) . S O ( , , ) . S O ( , , ) . S O ( , , ) O ff s e t [ a r c s e c ] Velocity [kms ] . S i O ( ) . H C O + ( ) . H C N ( ) . H C N ( ) . S O ( , ) . S O ( , , ) . S O ( ) . H C O + ( ) . C H O H ( E ) . C H CC H ( ) . S O ( , , ) . S O ( , , ) . S O ( , , ) . S O ( , , ) . S O ( , , ) . S O ( , , ) O ff s e t [ a r c s e c ] Velocity [kms ] F i g . . — L e f t : P V p l o t s f o r s e v e r a l o b s e r v e d li n e s a l o n g t h e o u t fl o w a x i s . T h e n a m e o f t h e li n e a nd t h e p e a k e m i ss i o n i ss h o w n . R i g h t : P V p l o t s f o r s e v e r a l o b s e r v e d li n e s p e r p e nd i c u l a r t o t h e o u t fl o w a x i s . T h e n a m e o f t h e li n e a nd t h e p e a k i n t e n s i t y i ss h o w n ,i nun i t s o f J y b e a m − . Herv´ıas-Caimapo et al.
TABLE 1Physical conditions estimates on the averaged spectra.
Location Using SiO Using SO Column density Rot. temperature Column density Rot. temperatureblue/red wing blue/red wing blue/red wing blue/red wing[10 cm − ] [K] [10 cm − ] [K]50% emission peak 9 . ± . . ± . ± ± † ± . ± . ± ± . ± . . ± . ± ± ± ± . ± . . ± . ± ± ± ± . ± . . ± . ± ± ± ± Note . — † This values are estimated in Merello et al. (2013b) with APEX data, which covers thesource completely. Therefore, this is a rough estimate.
Ambient core OutflowOutflow Expanding shellCavity
Fig. 4.—
Working model of the G331 .512-0.103 massive coreand outflow. The main morphological features are shown. E u / k [K]252627282930 l n ( N u / g u ) Best linear fitLower/upper 80% confidence limitRed wingBlue wing
Fig. 5.—
Example rotational diagram of the SO lines, corre-sponding to the averaged spectra inside the 50% of the peak emis-sion contour. The blue and red points correspond to the spectrumintegrated in the blue and red wing. The solid straight line corre-sponds to the linear regression fit. The dashed lines are the 80%confidence intervals. (J= 16 , − , ) has a clear excess over the otherlines, and may be contaminated by an unknown line. Intotal, 7 lines were used in the blue and red wings.We consider the spatially averaged spectrum inside thecontour of the 50% of the peak emission, as well as inside Blue wing Blue wing Red wing Red wing log ( N SO /cm )
120 140 160 180 200 220 T rot /K log ( N SO /cm )
120 140 160 180 200 220 T rot /K Fig. 6.—
Column density (left column) and rotational tempera-ture (right column) maps made with SO rotational diagrams. Theblue and red wings velocity ranges are −
120 to − . − . −
67 km s − . a circle with 1 . (cid:48)(cid:48) σ errors are cal-culated from the covariance matrix from the linear fit.The blue and red peaks show the highest column den-sity, of 9 × and 8 × cm − , respectively. Thisis consistent with the same tendency observed in SiO.The average rotational temperature in the SO emittingregion is ∼
170 K. The red peak is hotter than the bluepeak, 202 and 146 K, respectively. As an example, in Fig-ure 5, we show the SO rotational diagram for the fluxaveraged inside the 50% peak emission contour. For thelinear fit, we consider a conservative error of 20% for theln( E u /g u ) values. We show the 80% confidence intervals.The estimated SO column density and rotational tem-perature maps for each outflow wing are shown in Fig-ure 6. The ring-like emission is traced in SO by the blueand red wing emission. The rotational temperature mapsshow an elongated peak at >
200 K on both wings. TheLMA Observations of G331.512-0.103 II 7bright peak of emission, seen in most of the SO lines,corresponds to the peak in temperature in the blue wing,and also seen in N SO in the red wing. This means thatSO is likely tracing the most dense, hot, and shockedpart of the proto-stellar core.4.1.2. Fractional abundances
One physical parameter that is hard to measure, yetneeded for modeling is the fractional abundance of eachmolecule. The fractional abundance X M , where M stands as a particular molecule, is measured as X M = N M N H where N M is the column density of the correspon-dent molecule.In order to calculate the fractional abundances, someevaluation of the mass is needed. The virial mass canbe used to estimate it (e.g. Shirley et al. 2003), or themass estimated from the thermal continuum emission,the dust mass , can also be used (e.g. Fontani et al.2002). The latter method is used in this work, and themass is calculated from the 862 µ m continuum emission.The equation to calculate the gas mass from dust ther-mal continuum, when the emission is optically thin, is M d = S ν D Rκ ν B ν ( T d ) , (2)where S ν is the observed flux, κ ν is the dust mass coef-ficient, D is the distance to the source, B ν is the Planckfunction, T d is the dust temperature and R is the gas todust ratio, which is uncertain, but usually set to 100in the literature. The following values will be used: κ ν = 1 .
89 cm g − (interpolated from model OH5 in Os-senkopf & Henning 1994) and D = 7 . M cd ) can be estimated from M cd = AN M µm H X M , (3)where µm H is the mass of the hydrogen atom times themean molecular weight (that accounts for a fraction ofhelium, equal to 2.29), A is the area of the emitting re-gion (which must be consistent with the area used tocalculate the dust flux S ν ) and N M corresponds to themeasured column density of the molecule M . Equating M cd to M d , and using the fact that Ω = A/D corre-sponds to the solid angle that the source subtends, wehave X M = Ω N M µm H κ ν B ν ( T d ) S ν R . (4)One quantity that is highly uncertain and sensitivefor the dust mass estimation is the dust temperature.MM13b estimates the dust temperature to be 35 K,based on the Spectral Energy Distribution (SED) of thesource. This value is for the entire clump though, so thedust temperature at the ALMA resolution is probablyhigher. A rotational diagram analysis of the CH CCHlines shows that the ambient core temperature is ∼
70 K.MM13a uses an equilibrium dust temperature of 400 Kwithin the central arcsec. We use these two limits to cal-culate a range of possible abundances, listed in Table 2.We compare our fractional abundances results withsimilar sources studied in the literature. As stated in MM13a, one of the most similar sources to G331.512-0.103 is G5.89-0.39, an UCHII region which exhibitsa powerful and compact bipolar molecular outflow.Klaassen et al. (2006) mapped this outflow with theJames Clerk Maxwell Telescope (JCMT) at ∼ . (cid:48)(cid:48)
0. Thefractional abundances for SiO and SO are listed in Ta-ble 2, column 4. The Orion-KL region (Genzel & Stutzki1989) has shown chemical richness. Tercero et al. (2010)and references therein, reported a spectral line surveyin the 90-300 GHz range. We list their fractional abun-dances, for comparison, in Table 2, column 5. Similarly,the UCHII region G34.26+0.15 (Garay et al. 1986; Wood& Churchwell 1989) is also a well studied massive starforming source. For comparison, the fractional abun-dances are listed in Table 2, column 6.4.1.3. Density and column density estimation with
RADEX
The radiative transfer code
RADEX (van der Taket al. 2007) is a non-LTE radiative code that uses ra-diative/collision rates information and basic geometriesto solve the statistical equilibrium and radiative transferequations, estimating the intensities of several transitionsof the most typical molecules observed in the sub-mmrange.Since the statistical equilibrium/population levelsand the radiative transfer/radiation field are coupled,
RADEX uses the Large Velocity Approximation (LVG)approximation by Sobolev (1960), where the mean inten-sity ¯ J is expressed as a function of the source function S and a photon escape probability β . All the knowledge ofthe geometry/optical depth goes into this probability β .An estimate is given by β ∼ (cid:104) exp( − τ ) (cid:105) = 1 τ (cid:90) τ exp( − τ (cid:48) ) dτ (cid:48) = 1 − exp( − τ ) τ ,(5)which coincides with the expression for a radially expand-ing sphere.The procedure consists of using the integrated inten-sity ratio of several SO lines to constrain the physicalparameters (e.g. Fu et al. 2012). Since RADEX doesnot include data on every SO transition, and also somelines of SO cannot be used (they are blended withneighboring lines in the blue and/or red wing), only4 lines were used for this analysis. The lines shouldbe as far away as possible in terms of energy, for thesame reason that in the rotational diagram case. Twoof them have low and two of them have high upper-level energies ( E u ), hence 4 line ratios are defined: (J=20 , − , )/(J= 8 , − , ), (J= 20 , − , )/(J=11 , − , ), (J= 25 , − , )/(J= 8 , − , ) and(J= 25 , − , )/(J= 11 , − , ).The code needs several input values. Density and col-umn density are left as free parameters. The kinetic tem-perature is fixed to 170 K, which is the value constrainedby SO . The line-width is set to 16 . − , which is theaverage of the FWHM of the 4 lines. The backgroundtemperature is set to 2 .
73 K. The main collision partnerthat
RADEX handles is H . The output of the code foreach pair of density and column density values are theline temperature T R , the excitation temperature T ex , theoptical depth τ of the line and the integrated intensityassuming a Gaussian shape. The latter is used to calcu- Herv´ıas-Caimapo et al. TABLE 2Range of abundances for the SiO and SO molecules. Fractionalabundances in other massive outflow sources. For reference, the spatialresolution of our ALMA data is ∼ kAU. G331, T d =70 K G331, T d =400 K G5.89-0.39 Orion-KL G34.26+0.15SiO 1 . −
9) 1 . −
8) 3( − a . − b ∼ − / d SO . −
7) 7 . −
7) 2( − a . − c . − . − e Note . — a Klaassen et al. (2006), JCMT, spatial resolution ∼
30 kAU b Tercero et al. (2011), IRAM 30m, s. r ∼ c Esplugues et al. (2013), IRAM 30m, s. r. ∼ d Hatchell et al. (2001), VLA, s. r. ∼
10 - 84 kAU e Mookerjea et al. (2007), BIMA, s. r. ∼ TABLE 3Points probed in SO lines to test the χ minimization withRADEX and the estimated properties. Location α δ log ( n ) log ( N (SO )) † J2000 J2000 [cm − ] [cm − ]SiO blue peak 16:12:09.91 -51:28:37.4 9 . . . . . . Blue Peak =0.02 =0.18 Red Peak =0.08 =3.46 Cavity Center =0.05 =1.90 log ( N ) [cm ] l o g ( n ) [ c m ] Fig. 7.— χ map in the density - column density parameter space. The dashed, solid, dotted and dashed-dotted lines represent the 0.5,1, 2 and 3 σ confidence intervals. The red and blue colors represent the fits for the red and blue wings emission, respectively. The labels foreach position are listed in Table 3. late the line intensity ratio.Our method is similar to the one used by Plume et al.(1997) and van der Tak et al. (2000). It consists of thefollowing: given N line ratios (in this case N = 4) andgiven the RADEX model that will predict line ratios asa function of density and column density, the quantitythat must be minimized is χ = (cid:88) i ( R i − R m,i ( n (H ) , N (SO ))) σ i , (6)where i stands for each of the defined 4 ratios, R i is theobserved line ratio, R m,i is the RADEX modeled lineratio and σ i is an estimation of the error in the line ratio.The wing emission in the SO lines must be usedto measure the line ratios. These ranges were set as − . − . − for the blue wing and − . −
50 km s − for the red wing. These limits are slightlydifferent than the ones defined in Section 4.1.1, since weare considering a subset of SO lines. We estimate the σ i errors in the following way: the error due to the RMSof the spectra is ∼ ∼ ∼ . (cid:48)(cid:48) χ =LMA Observations of G331.512-0.103 II 9 offset from mm peak [arcsec] F ( x ) [ J y b e a m ] n ( r ) r p F ( x ) x m Fig. 8.—
Diamonds: ALMA 862 µ m continuum radial profilemeasured with respect to the peak. The two vertical lines are thebeam size and the 5 σ contour radius from the map (the radius R ofthe source). Blue line: Fit of a power law to the flux radial profile.Only the points inside the two vertical lines were considered. Blackline: Flux radial profile calculated using an isothermal gas spherewith a density power law. See Section 4.1.4 for details. χ ( n H , N SO ) space, the confidence levels are estimatedas contours at chosen σ . The results for each of thelocations probed are shown in Figure 7. The contourscorresponding to 0.5, 1, 2, and 3 σ are shown as dashed,solid, dotted, and dashed-dotted lines, respectively. Theresults are tabulated in Table 3. The blue and red peakconditions are representative of the conditions in the out-flow. The measured density is quite high ( n (cid:38) cm − ).These conditions are representative of the dense andshocked gas in the outflow, as traced by SO .4.1.4. Density radial gradient, dust continuumobservations
One way to estimate the density distribution of theproto-stellar core is to use the thermal dust continuumobservations: in our case, a line-free band in our ALMAdata set. In the literature, this is usually performed as-suming that the density and temperature will follow apower law with the radius, i.e. n ∝ r − p and T ∝ r − q . Forexample, in Garay et al. (2007), 19 IRAS point sourceswere imaged with the SEST telescope in 1.2 mm. Theyfind that most of these sources are single peaked, with adensity power-law index in the range p ∼ . − .
2, whichmeans that massive dense cores can be described as cen-trally condensed. In Beuther et al. (2002), they charac-terize 69 massive star forming regions using the 1.2 mmcontinuum emission. They conclude that a unique power-law is not enough to characterize the radial distribution,but a flat center, followed by an inner and an outer ra-dial power-law is adequate. The mean value for the innerpower-law they found is p ∼ . ± .
5. However, all ofthese studies are using single-dish observations and areappropriate on clump scales of ∼ .
25 pc. The theorysays that assuming a density power-law n ∝ r − p , a tem-perature distribution T ∝ r − q and an observed flux dis-tribution F ∝ r − m , in the optically thin emission, thethree coefficients are related by m = p + Qq − Q is a frequency correction factor ( ∼ . µ m continuum ALMA observations, aradial profile is calculated using contours of equal flux.Given two contiguous levels, all the pixels in betweenthe two levels are averaged. This flux is associated to aradius ( b i + b i +1 ) / b i = (cid:112) A i /π , i.e., the averageof the two radii associated to the effective area betweencontour i and i + 1. As external radius of the clump, weuse 2 . (cid:48)(cid:48)
1, which is the radius of the 5 σ contour level inthe 862 µ m emission. The radial profile is presented inFigure 8. Assuming a density function n ( r ), the columndensity in a sphere can be calculated straightforwardlywith the following equation (Dapp & Basu 2009) N ( x ) = 2 (cid:90) Rx n ( r ) rdr √ r − x , (7)where x is the angle distance with respect to the peak(the impact parameter), r is the radial distance with re-spect to the center and R is the sphere radius. In the op-tically thin limit, the flux is given by (Kauffmann et al.2008) F ν ( x ) = B ν ( T d ( x )) κ ν µm H N ( x ), (8)where B ν is the Planck spectral law and T d ( x ) is the dusttemperature.A zeroth-order approximation is to consider the dusttemperature of the core to be constant. On the one hand,our estimation of temperature comes from the rotationaldiagram of CH CCH, which indicates that the core has atemperature around 70 K. On the other hand, a direct fitto the observed flux radial profile with a power-law givesan index of m = 1 . ± .
2. If we consider T d = 70 K inequation 8, we can fit a density power-law with an index p = 0 . ± . ∼
70 K and looking intothe literature for massive star-forming cores profiles. Weperform the following steps:1. We assume that the temperature profile follows apower-law with index q = 0 .
4, value used for high-mass star forming cores in the literature (e.g. vander Tak et al. 2000; Garay et al. 2010).2. We assume that the mean temperature is 70 K be-tween the radii of 0 . (cid:48)(cid:48) . (cid:48)(cid:48)
1. This is (cid:104) T d (cid:105) = 70K = 12 . (cid:48)(cid:48) − . (cid:48)(cid:48) (cid:90) r =2 . (cid:48)(cid:48) r =0 . (cid:48)(cid:48) Ar − . dr . (9)3. Under these conditions, the dust temperature pro-file is given by T ( x ) ∼ x/ . (cid:48)(cid:48) − . .In this more realistic approach, the fitted density modelis a power-law with index p = 0 . ± .
1. This measureddensity power law is less steep than values reported formassive dense clumps, which could be due to the out-flow moving gas outwards, therefore flattening the den-sity profile. We note the lack of large scale emission0 Herv´ıas-Caimapo et al.by the interferometric observations (missing scales above11 . − . (cid:48)(cid:48)
0) and the flattening at small scales towardsthe center.The fitted dust density power law is n ( r ) = 1 . × ( r/ . (cid:48)(cid:48) − . cm − , so assuming a gas-to-dust ratio of100, the mean gas density is ∼ . × cm − , whichagrees with the high-density environments reported formassive star formation sources in the literature ( > cm − in small scales (cid:46) .
05 pc, Evans 1999). This den-sity is ∼ µ m thermal dust continuum traces all the warmgas coupled with the dust. We would expect this to beless dense than the shocked and compressed outflow gas.4.1.5. SiO shocks: model and properties
Following the model of Gusdorf et al. (2008), we com-pare our observations of SiO(J= 8 −
7) with the rotationalspectra simulated in that work. They give a grid of mod-els that allows to constrain the properties of the shockand outflow.In this model, the release of SiO and related moleculesis through the erosion of charged dust grains and theirice mantles by collision with neutral particles driven by asteady-state C-type shock. They consider multiple vari-ables, such as the dynamics of the dust grains, the accu-rate description of the sputtering, the thermal balance,and a complex gas and solid phase chemical network. Fi-nally, the intensities of the emission of the SiO spectrumare calculated through a LVG code (Section 4.1.3).Leurini et al. (2013) applies this model to the MYSOIRAS 17233-3606, Klaassen et al. (2006) uses the relatedmodel of Schilke et al. (1997b) for the G5.89-0.39 out-flow. They conclude that the SiO rotational emissionspectra can be modeled in shocks located in high-massstar forming regions, using the tools developed for morequiescent low-mass star forming regions, but using highervalues for pre-shock densities. Gusdorf et al. (2016) usesthe shock model for Cepheus A, a massive star nursery,but for CO and OH line observations.The parameters of the model that can be constrainedand affect the model significantly are the pre-shock den-sity n H and the shock velocity v s . The transverse mag-netic field and the viewing angle are also considered inthe model, but found to have a minor influence.For this analysis, we use the single-beam ( ∼ . (cid:48)(cid:48)
0) ob-servations of SiO(J= 7 −
6) and (J= 8 −
7) from MM13bobtained with APEX. According to MM13a, the differ-ences between the single APEX spectrum and the inte-grated ALMA spectrum of SiO(J= 8 −
7) are less than10%. The (J= 7 − −
7) line ratio from theAPEX data is 0 . ± .
02 (red wing) and 0 . ± . J up , only two models can give an SiO(J= 8 −
7) inten-sity greater than the SiO(J= 7 −
6) one. Both have n H = 10 cm − and v s = 32 and 34 km s − . The mod-eled line ratios are ∼ .
93 and ∼ .
88, respectively. Weconclude that the shocks traced by the SiO observationshave a speed of v s ∼
34 km s − and pre-shock densities of n H ∼ cm − . Similar values are found by other worksin massive outflows (Leurini et al. 2013, 2014; Gusdorf ellipse azimuth [deg] v e l o c i t y [ k m s ] kms ] Fig. 9.—
Left: On grayscale, the azimuthal PV plot of theSiO(J= 8 −
7) line. The red crosses represent the velocities wherethe peak emission is for each ellipse azimuth. The blue horizon-tal line is the fitted V and the white line is the systemic velocity − . − . The blue line corresponds to the best fitted V Z ( θ ).Right: On grayscale, the SiO(J= 8 −
7) 0th moment map at am-bient velocities. The red crosses trace the observed ring. The blueline is the best fitted projected ellipse. et al. 2016). 4.2.
Geometry and kinematics
Radial and Azimuthal PV plots
PV plots along a straight line axis can give informationonly along that particular axis, so it is helpful to plotthe same information in a different fashion. We use thetwo natural directions of polar coordinates: radial andazimuthal. The radial PV plot is constructed by choosinga center point, and the pixels at the same radius withina ring are averaged, using the function kshell , part ofthe KARMA package (Gooch 1995). The azimuthal PVplot is constructed by fitting an ellipse to trace the ring-like emission in the integrated intensity map. At eachazimuthal angle around the perimeter of the ellipse, weaverage in a slit of 11 pixels, perpendicular to the tangentof the ellipse (L´opez-Calder´on et al. 2016). For a fulldescription, see Section 4.2.2.An expanding shell geometry should show an “invertedC” pattern in the radial PV plot (Purcell et al. 2009). Atthe center of the shell, there is little emission at the sys-temic velocity, and at ± the expansion velocity we seethe emission. The radial PV plots for several observedlines are shown in Figure 10. In this figure, the SiOpanel presents the described “inverted C” profile, whichhints an expanding shell. The shape is not sharply de-fined though, so estimating the expansion velocity onlyfrom this plot is difficult. The peak of emission is clearlydefined, peaking at a radius of ∼ . (cid:48)(cid:48) ∼ .
036 pc atthe source distance). The SO and HCO + lines have avery similar radial PV plot, however the latter is com-pletely self-absorbed, and we cannot draw further con-clusions. H CO + , CH OH and CH CCH show the sys-temic velocity emission, but the first also shows weakhigh-velocity wings. The H CN and H CO + lines showa small self absorption dip at the systemic velocity. TheCH OH and the CH CCH lines have a single peak ofemission, located at a radius of ∼ . (cid:48)(cid:48) ∼ .
044 pc).In H CO + , the emission peaks at ∼ . (cid:48)(cid:48) ∼ .
062 pc),which makes H CO + the line that traces the most ex-tended emission.LMA Observations of G331.512-0.103 II 11 (16
4, 12
3, 13 ) HCO + (4-3)0.0 0.5 1.0 1.5 2.0 2.5 3.010810410096928884807672 H CO + (4-3) 0.5 1.0 1.5 2.0 2.5 3.0H CN(4-3) 0.5 1.0 1.5 2.0 2.5 3.0CH OH(4 -3 E ) 0.5 1.0 1.5 2.0 2.5 3.0CH CCH(21 -20 ) offset [arcsec] v e l o c i t y [ k m s ] Fig. 10.—
Radial direction PV plots for several observed lines. The contours are at 10, 25, 40, 55, 70, 85 and 95% of the peak value.The center of the cavity (the 0 coordinate of the X-axis) is α J = 16 h m s . δ J = − ◦ (cid:48) . (cid:48)(cid:48) SiO(8-7) SO (16
4, 12
3, 13 ) H CO + (4-3) CH OH(4 -3 E ) CH CCH(21 -20 ) azimuth offset [arcsec] v e l o c i t y [ k m s ] Fig. 11.—
Azimuthal direction PV plots for several observed lines. The contours are at 50, 65, 75, 85 and 95% of the peak value.
An expanding perfect sphere appears isotropic in anazimuthal PV plot. However, a shell, which in principlecan have an ellipsoidal shape, should be sinusoidal in theazimuthal PV plot (e.g. L´opez-Calder´on et al. 2016).For the same lines shown in Figure 10, the azimuthal PVplots are shown in Figure 11. The SiO line can be morereadily interpreted. Around −
90 km s − , the sinusoidaloscillation can be seen in the peaks as a function of az-imuth angle. In Section 4.2.2, an analytical ring modelis fitted to this sinusoidal oscillation, shown in Figure 9. The SO azimuthal plot has a similar shape. However, inthe latter case, the oscillation is at an offset velocity withrespect to the SiO. A sinusoidal shape could be invoked,but the fact that it is centered at a slightly bluer velocitythan the systemic −
90 km s − makes this plot harder tointerpret. The H CO + , CH CCH and CH OH lines arevery similar, showing clumpiness at the systemic velocity.4.2.2.
Ring model for the SiO(J = 8 − ) azimuthal PV plot Figure 11 shows that the azimuthal PV emission of theSiO line seems to trace some sort of expanding motion,2 Herv´ıas-Caimapo et al.
TABLE 4Best fitted parameters for the ringmodel.
Parameter ValueCenter J2000 16:12:10.0 -51:28:37.4Line of sight angle α -8 ◦ Position angle β ◦ Radius R . (cid:48)(cid:48) V
21 km s − Systemic velocity V -90 km s − as mentioned above. Using the best-fit ellipse points,together with the azimuthal PV plot, a simple mathe-matical model of a ring is fitted. When referring to ring,we mean the ring-like emission seen in SiO. The methodfollows from L´opez-Calder´on et al. (2016).Consider a ring in the plane XY of the sky expand-ing in the outwards direction at a constant velocity. Theequations that describe each XY coordinate and the pro-jected velocity in the Z direction perpendicular to theplane of the sky are X ( θ ) = X + R (cos θ cos α cos β − sin θ sin β ) (10) Y ( θ ) = Y + R (cos θ cos α sin β + sin θ cos β ) (11) V z ( θ ) = V + V cos θ sin α , (12)where θ is the angle along the perimeter of the ellipse(the azimuth angle), R is the radius of the ring, V is thesystemic velocity, V is the constant expansion velocityof the ring, α is the inclination respect to the plane ofthe sky, and β is the position angle (rotation of the XYplane). Using the best-fit ellipse, the observed X and Y coordinates are chosen for each ellipse azimuth angle asan intensity-weighted mean of the 10 pixels, belonging tothe perpendicular “slit”, that are considered when aver-aging the spectra. In this way, the best-fit ellipse actsas a guide for the “real” XY points on the ring. In thesame way, the Z direction projected velocities are chosenas the velocity where the peak of the emission takes placefor each ellipse azimuth angle, in the PV plot shown inFigure 11.There is a list of observed coordinates X , Y and V Z for each ellipse azimuth angle. Further, there is a listof X , Y and V Z values, calculated with equations 10-12and θ between 0 and 2 π . A function calculating thesum of the distance in the 3D space ( X , Y , V Z ) betweenthe observed and modeled coordinates is defined. Then,this function is minimized for the best fit 7 parameters: X , Y , α, β, R, V and V . The results are shown in Fig-ure 9. On the left panel, the peaks at each correspondingellipse azimuth are shown as red crosses. These trace theproposed inclined ring with their sinusoidal shape. Themodel of projected V Z velocity is shown as the solid blueline. On the right panel, the modeled projected ellipse inthe plane of the sky is shown against the SiO 0th momentmap at systemic velocities.The best fit parameters for the model are listed in Ta-ble 4. The center of the ellipse naturally will coincidewith the measured cavity center in MM13a. The angle α is small, which supports the hypothesis that the cavity-outflow complex is closely aligned with the line of sight.The expansion velocity of ∼
20 km s − is very similar to
200 175 150 125 100 75 50 25 0
Velocity [kms ] F l u x [ J y b e a m ] Observations blue peakModel blue peakObservations red peakModel red peak
Fig. 12.—
Comparison of averaged SiO(J= 8 −
7) spectra, ob-servations and outflow model. Each spectrum is averaged insideone beam. The blue and red color indicate spectrum on blueand red peaks ( ±
40 km s − channels), respectively. The solidlines represent observations. Dashed lines are for the MOLLIEmodel. The outflow model accounts for the high velocity channels | v | >
40 km s − . the one determined by MM13a, derived by only analyz-ing the position-velocity plot along or perpendicular tothe outflow axis.To measure the significance of the model, we performeda χ -test. The estimation of the error is critical for thispurpose, since the value of χ can change greatly de-pending on the error. The error is estimated from theGaussian distribution of the measured velocities aroundthe velocity where the peak is located. The method to es-timate the error in the determination of the peak velocityis the following: given a level of noise in the measurementof flux (or temperature) of the spectrum, ∆ T RMS , we cal-culate the error as the width in velocity in the gaussian-fit maximum temperature to the gaussian-fit maximumtemperature minus ∆ T RMS . That is, we look for ( v − v )such that T (1 − exp( − ( v − v ) / /σ T )) = ∆ T RMS , where T and σ T are the maximum temperature and the stan-dard deviation of the Gaussian fit, respectively. We setas the error for the χ -test σ v = v − v for each spec-trum at each azimuth angle. Then, we calculate a χ statistic only based on the projected velocity fit, givenby χ = Σ i ( V Z ( θ ) − V iZ σ v ) , where i runs through each of theoffset angles in the perimeter of the ellipse and V Z ( θ ) isgiven by eq. 12. The resulting reduced- χ is χ ν = 0 . overfit the data and theconsideration of errors is conservative.4.3. Radiative transfer modeling of the SiO(J = 8 − )line data cube This section presents a radiative transfer model of theSiO(J= 8 −
7) line emission using all the knowledge gath-ered in the two previous subsections. The parametersthat we use as input in the model are listed in Table 5.This model is performed with the 3D radiative transfercode
MOLLIE . The full detailed description of how themodel is generated is in Appendix B. Here, we limit our-selves to present only the results.We simulated the emission in two steps: one for anexpanding shell, and one for a bi-conical outflow. Theresults of the outflow simulation are shown in Figure 12.LMA Observations of G331.512-0.103 II 13
200 175 150 125 100 75 50 25 0
Velocity [kms ] F l u x [ J y b e a m ] Observations cavityModel cavity
Fig. 13.—
Comparison of observed and modeled SiO(J= 8 − −
7) observed spectrum averaged over one beam insidethe shell. The green dashed line is the modeled spectrum averagedover one beam inside the shell. The vertical line corresponds tothe systemic velocity.
We show the spectrum averaged over 1 beam (1 . (cid:48)(cid:48) −
89 km s − , but rather ∼ − blue-shifted.The modeled spectrum shows this same behavior, eventhough the intensity might not reproduce exactly the ob-served levels. It should also be noted that the red peakof the outflow is closer to the cavity center than the bluepeak, therefore explaining the asymmetry between theblue and red wings in the spectrum in Figure 13.In order to compare the emission of the SiO line obser-vations with the modeled data cube, we add up boththe shell and outflow models, since they are at non-overlapping velocities. However, we are ignoring anypossible radiative interaction between the outflow andthe shell. We compare with the observations using anintegrated emission map at systemic velocities and theposition-velocity plot along the outflow axis. Figure 14shows both of these. In general terms, the model is ableto qualitatively reproduce the observations. Therefore,we can state that the following is consistent with theG331.512-0.103 massive outflow:1. A bipolar outflow with an opening angle of ∼ ◦ ,a maximum velocity around 85 km s − at the axisthat decreases with the cylindrical radius, as statedin Stahler (1994); and with an inclination of ∼ − ◦ with respect to the line of sight. The interiorof the outflow has low density ( ∼ cm − ) andhigh velocity. The outer layer of the outflow, thatis in contact with the ambient core gas, has high TABLE 5Summary of the properties derived for thesource.
Property ValueDistance 7.5 kpcMass of core ∼ M (cid:12) Mass outflow lobes ∼ M (cid:12) eachKinetic age ∼ ±
70 kms − Expansion velocity of cavity ∼
21 kms − Mean density 4 . × cm − Density power law index ∼ . ∼
70 KOutflow inclination ∼ ◦ density ( ∼ − cm − ) and low velocity.2. A cavity, with low density ∼ cm − andmedium-velocity ( ∼
20 km s − ), with a radiallyoutwards expansion. This structure has a radiusof ∼ . (cid:48)(cid:48) −
7) in the position ofthe ring-like emission at the systemic velocities isclearly stronger than the emission in the cavity cen-ter. This translates to a spectrum with a singlepeak, rather than a dip, along the ring-like emis-sion. The model is not able to reproduce this, sincethe dynamics of the expanding motion is present inthe entire central region of emission, hence a dou-ble peak spectrum with a dip at systemic veloc-ity is present in the entire emission region of themodel, both in the cavity and in the ring-like emis-sion perimeter.2. Peaks of emission and structures are present withinthe ring-like emission in the observations. For this,the observed integrated emission map at the sys-temic velocities shows a clear contrast between thering-like emission and the cavity center, i.e. theintegrated emission is about twice in the ring-likeemission compared with the cavity center. How-ever, the levels of emission inside the cavity centerare similar in both the model and the observations.To improve the model, possible inhomogeneitiespresent in the shell must be accounted for. The factthat the data show some level of clumpiness at the ring-like emission of SiO(J= 8 −
7) indicates that the shellis not homogeneous and there is some structure beyondour simple analytical model. SUMMARYThe G331.512-0.103 molecular core contains one ofthe most luminous and powerful massive outflows har-bored in a high-mass star forming region in our Galaxy.We analyzed several molecular spectral lines observed at ∼ . (cid:48)(cid:48) I n t e g r a t e d e m i ss i o n S i O ( - ) m o d e l [ J y b e a m k m s ] Integrated emission SiO(8-7) observations offset [arcsec] v e l o c i t y [ k m s ] PV plot along the outflowSiO(8-7) observations P V p l o t a l o n g t h e o u t f l o w S i O ( - ) m o d e l [ J y b e a m ] Fig. 14.—
Left: Integrated emission at ambient velocity range of the simulated data cube, on color scale. The dashed blue contourscorrespond to the 0th moment map at the same range in the SiO(J= 8 −
7) observations. The levels are 5, 10, 20, 30, 40, and 60 Jybeam − km s − . Right: Position-Velocity plot along the outflow axis with a “slit” of 1 . (cid:48)(cid:48)
54. On gray scale, the modeled SiO(J= 8 −
7) data cube.On dashed blue contours, the observations. The contour levels are 0.2, 0.86, 1.52, 2.18, 2.84, and 3.5 Jybeam − Based on the PV diagrams of the lines, there are twogroups. One that traces the broad high-velocity wings,and one that traces narrow systemic velocity emission.For the high-velocity outflow, we have determined itsproperties. The temperature of the most dense and hotgas in the outflow is ∼
150 K (blue lobe) and ∼
200 K(red lobe). The column densities of the SiO and SO areconstrained. These estimates are also supported by ananalysis with the code RADEX: a density of ∼ cm − is constrained for the outflow shocked gas. The fractionalabundance of SiO and SO is in agreement with valuesfound in other massive outflows in MYSOs.The properties of the systemic velocity emission, theambient core, are analyzed. The kinetic temperatureof the ambient core is T K ∼
70 K. The 862 µ m dustcontinuum emission can be well fitted with a densitypower-law with an index ∼ . n H = 4 . × cm − .The geometry and morphology of the ambient core ischaracterized by the peaks of the PV plots in the ra-dial direction: structures with radii of ∼ . (cid:48)(cid:48) . (cid:48)(cid:48) OH) and ∼ . (cid:48)(cid:48) CO + ). A ring model was fittedto the SiO(J= 8 −
7) azimuthal PV plot. The parame-ters of this model is an inclination angle of ∼ ◦ and anexpansion velocity of ∼
20 kms − .To model the source, composed of an outflow and anexpanding shell, we performed a radiative transfer modelof the SiO(J= 8 −
7) line using the code MOLLIE. Themodel is composed of two structures: a conical bipolaroutflow with a velocity field that scales with the cylin-drical radius, and an ellipsoidal cavity and expandingshell with a peak expansion velocity of ∼
20 kms − inthe spherically radial direction. The model was able to reproduce the main features at the wings (outflow) andambient velocity ranges, although the emission is notcompletely recovered. The strong peak that dominatesat systemic velocities, observed in a ring-like emissionsurrounding a cavity, is not totally accounted for in themodel.Our global scenario for the source is the following: Atthe center position, where a newborn massive star is lo-cated, there are structures proper to the ambient core,which is at systemic velocities. The proto-star shows amassive, bipolar, and high-velocity outflow with veloci-ties of ±
70 kms − , likely powered by a collimated jet andby the accretion of gas onto the source. The powerfulstellar winds and ionizing radiation from the proto-starpush against the ambient core gas, inflating a cavity andan expanding shell-like structure.C.H.C. acknowledges support by CONICYT Beca deMagister Nacional, folio 221220026, and partial sup-port by FONDECYT project 1120195. M.M. acknowl-edges support from the grant 2017/23708-0, S˜ao PauloResearch Foundation (FAPESP). L.B. and G.G. ac-knowledge support from CONICYT project Basal AFB-170002. We thank Al Wootten and the staff of NRAO fortheir help and assistance with the reduction of ALMAdata. This Paper makes use of the following ALMAdata: ADS/JAO.ALMA REFERENCESAdams, F. C. 1991, ApJ, 382, 544Bally, J. 2016, ARA&A, 54, 491 Beuther, H., Schilke, P., Menten, K. M., et al. 2002, ApJ, 566, 945
LMA Observations of G331.512-0.103 II 15
Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001,MNRAS, 323, 785Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 735Brogan, C. L., Chandler, C. J., Hunter, T. R., Shirley, Y. L., &Sarma, A. P. 2007, ApJ, 660, L133Brogan, C. L., Hunter, T. R., Indebetouw, R., et al. 2008,Ap&SS, 313, 53Bronfman, L., Garay, G., Merello, M., et al. 2008, ApJ, 672, 391Charnley, S. B. 1997, ApJ, 481, 396Claude, S. M. X., Avery, L. W., & Matthews, H. E. 2000, ApJ,545, 379Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092Draine, B. T., & McKee, C. F. 1993, ARA&A, 31, 373Esplugues, G. B., Tercero, B., Cernicharo, J., et al. 2013, A&A,556, A143Evans, II, N. J. 1999, ARA&A, 37, 311Fontani, F., Cesaroni, R., Caselli, P., & Olmi, L. 2002, A&A, 389,603Frank, A., Ray, T. P., Cabrit, S., et al. 2014, Protostars andPlanets VI, 451Fu, R. R., Moullet, A., Patel, N. A., et al. 2012, ApJ, 746, 42Garay, G., & Lizano, S. 1999, PASP, 111, 1049Garay, G., Mardones, D., Bronfman, L., et al. 2010, ApJ, 710, 567Garay, G., Mardones, D., Brooks, K. J., Videla, L., & Contreras,Y. 2007, ApJ, 666, 309Garay, G., Rodriguez, L. F., & van Gorkom, J. H. 1986, ApJ,309, 553Garc´ıa, P., Bronfman, L., Nyman, L.-˚A., Dame, T. M., & Luna,A. 2014, ApJS, 212, 2Genzel, R., & Stutzki, J. 1989, ARA&A, 27, 41Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209Gooch, R. 1995, in Astronomical Society of the Pacific ConferenceSeries, Vol. 77, Astronomical Data Analysis Software andSystems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 144Gusdorf, A., Cabrit, S., Flower, D. R., & Pineau Des Forˆets, G.2008, A&A, 482, 809Gusdorf, A., G¨usten, R., Menten, K. M., et al. 2016, A&A, 585,A45Hatchell, J., Fuller, G. A., & Millar, T. J. 2001, A&A, 372, 281Hunter, T. R., Brogan, C. L., Indebetouw, R., & Cyganowski,C. J. 2008, ApJ, 680, 1271Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, II, N. J., &Lee, C. W. 2008, A&A, 487, 993Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531Klaassen, P. D., Plume, R., Ouyed, R., von Benda-Beckmann,A. M., & Di Francesco, J. 2006, ApJ, 648, 1079Leurini, S., Codella, C., Gusdorf, A., et al. 2013, A&A, 554, A35Leurini, S., Codella, C., L´opez-Sepulcre, A., et al. 2014, A&A,570, A49Linke, R. A., Frerking, M. A., & Thaddeus, P. 1979, ApJ, 234,L139Looney, L. W., Mundy, L. G., & Welch, W. J. 2003, ApJ, 592, 255L´opez-Calder´on, C., Bronfman, L., Nyman, L.-˚A., et al. 2016,A&A, 595, A88 L´opez-Sepulcre, A., Codella, C., Cesaroni, R., Marcelino, N., &Walmsley, C. M. 2009, A&A, 499, 811McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap,K. 2007, in Astronomical Society of the Pacific ConferenceSeries, Vol. 376, Astronomical Data Analysis Software andSystems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127Merello, M., Bronfman, L., Garay, G., et al. 2013a, ApJ, 774, L7—. 2013b, ApJ, 774, 38Molinari, S., Merello, M., Elia, D., et al. 2016, ApJ, 826, L8Mookerjea, B., Casper, E., Mundy, L. G., & Looney, L. W. 2007,ApJ, 659, 447Nagy, Z., van der Tak, F. F. S., Fuller, G. A., & Plume, R. 2015,A&A, 577, A127Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943Plume, R., Jaffe, D. T., Evans, II, N. J., Mart´ın-Pintado, J., &G´omez-Gonz´alez, J. 1997, ApJ, 476, 730Prasad, S. S., & Huntress, Jr., W. T. 1980, ApJ, 239, 151Purcell, C. R., Minier, V., Longmore, S. N., et al. 2009, A&A,504, 139Rawlings, J. M. C., Redman, M. P., Keto, E., & Williams, D. A.2004, MNRAS, 351, 1054Rodriguez-Franco, A., Martin-Pintado, J., & Fuente, A. 1998,A&A, 329, 1097Schilke, P., Groesbeck, T. D., Blake, G. A., Phillips, & T. G.1997a, ApJS, 108, 301Schilke, P., Walmsley, C. M., Pineau des Forets, G., & Flower,D. R. 1997b, A&A, 321, 293Shirley, Y. L., Evans, II, N. J., Young, K. E., Knez, C., & Jaffe,D. T. 2003, ApJS, 149, 375Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23Sobolev, V. V. 1960, Moving envelopes of starsStahler, S. W. 1994, ApJ, 422, 616Tafalla, M., & Bachiller, R. 2011, in IAU Symposium, Vol. 280,The Molecular Universe, ed. J. Cernicharo & R. Bachiller,88–102Tan, J. C., Beltr´an, M. T., Caselli, P., et al. 2014, Protostars andPlanets VI, 149Tercero, B., Cernicharo, J., Pardo, J. R., & Goicoechea, J. R.2010, A&A, 517, A96Tercero, B., Vincent, L., Cernicharo, J., Viti, S., & Marcelino, N.2011, A&A, 528, A26van der Tak, F. F. S., Black, J. H., Sch¨oier, F. L., Jansen, D. J.,& van Dishoeck, E. F. 2007, A&A, 468, 627van der Tak, F. F. S., van Dishoeck, E. F., Evans, II, N. J., &Blake, G. A. 2000, ApJ, 537, 283Wakelam, V., Caselli, P., Ceccarelli, C., Herbst, E., & Castets, A.2004, A&A, 422, 159Wood, D. O. S., & Churchwell, E. 1989, ApJS, 69, 831Zapata, L. A., Ho, P. T. P., Schilke, P., et al. 2009, ApJ, 698, 1422Zhang, Y., Tan, J. C., & McKee, C. F. 2013, ApJ, 766, 86Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481
APPENDIX
A: INTEGRATED ALMA BAND 7 SPECTRA OF G331.512-0.103
Figure 15 shows the integrated spectra obtained with ALMA band 7 over the frequency range 345-348 GHz (SW3-SW2), and 356.5-359.5 GHz (SW0-SW1). Here we show all identified lines. The ones with blue labels correspond tothe lines presented in MM13a, the ones in red labels are the 18 lines newly analyzed in this study and the ones withblack labels are the rest of the lines that will be used in future studies.
B: DETAILED DESCRIPTION OF MOLLIE MODELING OF THE SOURCE
Estimation of physical conditions
Since physical conditions can vary across the source, and multiple combination of geometries and parameters canreproduce an observed spectrum, we need a rough estimate of these conditions in order to constrain how the modelis set up. In the following, we discuss broad estimates of the necessary physical conditions that must be specified ina MOLLIE grid model. For a 3D grid, we must specify for each voxel the density, the kinetic temperature, the locallinewidth, the 3 components of the velocity vector and the fractional abundance of the molecule that is being observed.The estimation of density is challenging. In Section 4.1, we estimated the density of the ambient gas with the 862 µ mdust continuum and the density of the region emitting in SO , tracing the outflow wings and therefore the high-velocityand shocked gas. In the following models, the “background” density, i.e the density of the ambient core as a functionof radius, is set to the power law estimated from the dust continuum flux in Section 4.1.4, n ( r ) = 2 . × [cm − ]( r/ . (cid:48)(cid:48) − . , (B1)6 Herv´ıas-Caimapo et al. . . . . . . . SO ( , − , ) SO ( , − , ) CH OCH ( , − , ) SO ( , − , ) SO ( , − , ) H CN( − ) SO ( , − , ) SO ( , − , ) SO ( , − , ) HC N( − ) CO( − ) CH OH( , − , ) SO ( , − , v ) SO ( , − , ) SO ( , − , ) H CO + ( − ) SiO( − ) . . . . . . . CH OCH ( , − , ) HCO + ( − ) SO ( , − , ) SO ( , − , ) SO ( , − , ) SO ( , − , ) CH OCH ( , − , ) SO ( , − , ) SO ( , − , ) Unidentified SO ( , − , ) SO ( , − , ) SO ( , − , ) SO ( , − , ) SO ( , − , ) SO ( , − , ) SO ( , − , ) Unidentified CH OH( − E ) S O( − ) Unidentified CH CCH( K =3 − K =3 ) CH CCH( K =2 − K =2 ) CH CCH( K =1 − K =1 ) CH CCH( K =0 − K =0 ) SO ( , − , ) SO ( , − , ) CH OCH ( , − , ) F r e q u e n c y [ G H z ] Flux [Jy] F i g . . — C o m p o s e dS W - S W (t o p ) a ndS W - S W ( b o tt o m ) b a nd s o f t h e i n t e g r a t e d s p ec t r ao b s e r v e d t o w a r d G . - . b s e r v e d w i t h A L M A . T h e m o s t p r o m i n e n t li n e s a r e m a r k e d a c r o ss t h e s p ec t r a . L i n e s u s e d i n t h i s w o r k a r ec o l o r e d . LMA Observations of G331.512-0.103 II 17where r is the radius expressed in arcsec. Note that this law corresponds to the “background density” in the model.It does not apply to the center of the model grid, where the shell or the outflow is to be modelled with a correspondingdifferent density law. For the shocked high density, i.e. the density of the most dense section of the outflow and thecavity/shell (where the stellar winds from the proto-star are impacting), we will use the order of magnitude valuesderived in Section 4.1.3 for the blue and red peaks as a reference. That is n = 10 cm − .We only have two estimates of temperature from 2 different molecules. One, the temperature estimated with theCH CCH molecule, tracing the core and systemic velocity environment. Since this molecule is a very good thermometer,we use it as the kinetic temperature of the ambient core. Then, we consider T K = 70 ± . Since this molecule is likely tracingthe most dense and shocked section of the gas, as stated before, we use it as the probe of outflow conditions. We usethe value T K = 150 K, as constrained by the rotational diagram.In the model, the line-width will increase from ∼ − at the edge of the structure emission at 4 . (cid:48)(cid:48) to ∼ − at the center of the cavity. This will be implemented with a linear velocity gradient.The fractional abundance of SiO is one of the most uncertain parameters. We will use the value estimated inSection 4.1.2, X SiO ∼ − . Since the SiO emission comes from the outflow and cavity regions, we cannot constrainthe fractional abundance of the cold ambient core gas with the same observations. A value of 10 − is used, consideringthat a value of ∼ − is cited as the abundance of SiO in cold dense/starless cores environments (Schilke et al. 1997b). Radiative transfer model
To simulate the response of the interferometer, the resulting simulated data cubes were processed through theCommon Astronomy Software Applications (CASA, see Section 2) tasks simobserve , to simulate a set of measuredvisibilities with the compact configuration of the ALMA array; and simanalyze , to produce synthetic deconvolvedimages from the visibilities.The high-velocity wings in the SiO spectra makes evident the presence of an outflow with ∼
70 kms − from thesystemic velocity of the cloud. The model that we will use for the outflow is a cone for each lobe. The estimation ofthe dimensions of this cone is performed in the following way: the spatial offset between the peaks of the ±
40 kms − spectral channels in the SiO emission is ∼ . (cid:48)(cid:48)
2. The model from Section 4.2.2 indicates that the cavity has an inclinationangle with respect to the line of sight of ∼ − ◦ . Assuming that the cavity and the outflow axis are aligned, thetotal extension of both outflow lobes is ∼ .
29 pc at the source distance along its axis. Therefore, the height of thecone is 0.145 pc. To estimate the opening angle (or equivalently, the base of the cone), we use the fact that the velocityat which the outflow starts in the red wing is − . − . At this velocity, the ring of emission is ∼ . (cid:48)(cid:48) . (cid:48)(cid:48)
2. This gives an opening angle of ∼ ◦ .On the inside of an outflow, the density is low and the velocity and temperature are high. On the contrary, thecore/envelope gas has high density and low temperature and velocity. Rawlings et al. (2004) modeled an outflowwith an inside-outflow density of 10 cm − . Zhang et al. (2013) simulated radiative transfer and SEDs of massivestar formation and considered the effects of outflows. Their densities inside the outflow are ∼ − cm − . Closeto the outflow axis, where the high-velocity collimated jet is located, the density can be somewhat higher. On theshocked region, close to the edge of the cone where the high-velocity flow is interacting with the quiescent core, weknow that the pre-shock density is close to 10 cm − for our outflow source, as calculated in Section 4.1.5. A shockhas an enhancement of 10-100 times the pre-shock density (Draine & McKee 1993). The size of a layer of very denseand shocked gas is taken as ∆ s = 5 × cm ∼ .
002 pc (Gusdorf et al. 2008), the typical length of a shock. Oneconstrain we have for the density is the total mass of the outflow. Each lobe has a total mass of ∼ M (cid:12) (Bronfmanet al. 2008), from lower resolution CO observations. The total mass will be estimated here by adding concentric disksapproximating the shape of the cone M T = µm H π (cid:88) z i ∆ z ( z i ) (cid:90) ψ cone n ( ψ ) ψdψ , (B2)where ψ is the cylindrical radius, ψ cone is the exterior cylindrical radius of the cone and ∆ z is the height of the diskat each particular z height.We will set a dense outflow region close to the edge of the cone. This will be limited by ψ = ψ cone − ∆ s for eachparticular z . The density is n ( ψ ) / cm − = (cid:40) . . ψ / ( ψ − ψ / if ψ < ψ . + 10 . enh − E ψ − ψ )) if ψ > = ψ , (B3) where enh is the enhancement factor of the density, from 10 to 100 in a shock length of 0.00162 pc. In this way,the density both close to the outflow axis and close to the cone limit is ∼ . cm − ; meanwhile, in the middle of theoutflow, it is nearly uniform and equal to 10 cm − . Finally, very close to the outflow limit there is a layer where thedensity is increased from 10 to 100 times the pre-shock density. Figure 16 (Left) shows the density profile. The totalmass of each lobe using this density is ∼ M (cid:12) .The kinematics of molecular outflows is discussed in Stahler (1994). In order to explain the PV plots of observedoutflows, the so called “outflow Hubble law”, Stahler proposes a series of equations characterizing the velocity distri-8 Herv´ıas-Caimapo et al. Z [ k A U ] l o g ( d e n s i t y ) [ c m ] | v e l o c i t y | [ k m s ] K i n e t i c t e m p e r a t u r e [ K ] Fig. 16.—
Left: Density distribution in cylindrical coordinates, used to simulate the outflow. The model is axysimmetric. Center: Thesame as the left, but showing the kinetic temperature profile. Right: The same as the left, but showing the module of the velocity field ateach location. The velocity vector has a ± Z direction. bution of an outflow. One conclusion is that the velocity field in an outflow depends on the cylindrical radius (if theZ-axis is the outflow symmetry axis). So the picture is that the highest velocities are found close to the outflow axisand then it decreases with the cylindrical radius. This picture is consistent with a high-velocity, highly collimated jet,close to the outflow axis. The velocity field will depend on the cylindrical radius (cid:107) (cid:126)v (cid:107) = 85 (cid:18) − ψ/ . (cid:48)(cid:48) . (cid:48)(cid:48) π/ (cid:19) [kms − ], (B4)where the cylindrical radius ψ is expressed in arcsec and 4 . (cid:48)(cid:48) π/
12) corresponds to the radius of the base of thecone. The velocity vector will have a direction in the ± Z-axis. The modulus of the velocity at each location is shownin Fig. 16 (Center).Using our estimates of temperature and the fact that outflows are somewhat hot compared with ambient core gas,we have a temperature of ∼ −
200 K for the shocked dense gas, and 70 K for the harboring core, estimatedfrom CH CCH. For comparison, the outflow in Rawlings et al. (2004) has a temperature of 50 K. The temperaturedistribution we use has ∼
70 K for most of the interior of the outflow and then increases to 150 K in the edge, wherethe gas is shocked, dense and hot. The temperature will be given by T ( θ ) = 70[K] exp(162( θ/rad ) ), where θ is inradians and corresponds to elevation angle (outflow opening angle). The distribution is shown in Fig. 16 (Right).A model for the shell emission is made independently. The expanding shell is modeled as two concentric ellipsoids.We know they are not spheres because the azimuthal PV plots are sinusoidal, and therefore there is an inclinationinduced asymmetry. The ellipsoids are oblate spheroids, with an aspect ratio of 1:1.5. The dimensions are 1 . (cid:48)(cid:48) . (cid:48)(cid:48) n ( r )cm − = (cid:40) . . r / ( r − r / if r < r . + 10 . enh − E r − r )) if r > = r , (B5) where r is the radius at which the most dense and shocked zone starts, located at the exterior limit of the innercavity. It is defined by r = r ellipsoid , inner − ∆ s . The density in the expanding shell, i.e. between the inner and outerellipsoids, is given by a power-law such that the border condition between the shell and the ambient core is fulfilled,that is, at r = r ellipsoid , outer , the density is given by n max ( r/r ellipsoid , inner ) − p shell = 2 . × [cm − ]( r/ . (cid:48)(cid:48) − . , where n max is the maximum density reached at the cavity r = r ellipsoid , inner . The power-law index of the density of the shellis therefore given byLMA Observations of G331.512-0.103 II 19 Z [ k A U ] l o g ( d e n s i t y ) [ c m ] | v e l o c i t y | [ k m s ] Fig. 17.—
Left: Density profile of the modeled cavity, shell and ambient core. The Z axis is the same axis as in Fig. 16, the symmetryaxis of the outflow. Right: The same as the left, but showing the magnitude of the velocity vector field. The direction of the velocity vectoris radially spherical. p shell = log(2 × ) − . r ellipsoid , outer / . (cid:48)(cid:48) − log( n max ) − log( r ellipsoid , outer /r ellipsoid , inner ) . (B6) The full density profile is shown in Fig. 17. The mass of the cavity plus the expanding shell using this densitydescription is ∼ M (cid:12) .Since the cavity is presumed to be blown-up by the stellar radiation output, the velocity will be set to the maximumvalue observed in the cavity, that is 20 kms − (see Section 4.2.2). The velocity is expected to diminish from themaximum value at r = r ellipsoid , inner to the systemic value, i.e. 0 kms − , within the expanding shell, that is at r = r ellipsoid , outer . This is because the ambient core is expected to be at the systemic velocity. To accomplish this, themodulus of the radial velocity is an exponential law, given by | v | = v max exp( − A ( r − r ellipsoid , inner )), (B7)where A is a constant that scales how fast the magnitude drops to zero. With this law, the velocity will be maximumat the edge of the cavity (or the inner ellipsoid of the expanding shell) and will approach zero (or the systemic velocity)at the outer edge of the expanding shell. The adopted value for A0