An Ordered Envelope-disk Transition in the Massive Protostellar Source G339.88-1.26

We report molecular line observations of the massive protostellar source G339.88-1.26 with the Atacama Large Millimeter/Submillimeter Array. The observations reveal a highly collimated SiO jet extending from the 1.3 mm continuum source, which connects to a slightly wider but still highly collimated CO outflow. Rotational features perpendicular to the outflow axis are detected in various molecular emissions, including SiO, SO2, H2S, CH3OH, and H2CO emissions. Based on their spatial distributions and kinematics, we find that they trace different parts of the envelope-disk system. The SiO emission traces the disk and inner envelope in addition to the jet. The CH3OH and H2CO emissions mostly trace the infalling-rotating envelope, and are enhanced around the transition region between envelope and disk, i.e., the centrifugal barrier. The SO2 and H2S emissions are enhanced around the centrifugal barrier, and also trace the outer part of the disk. Envelope kinematics are consistent with rotating-infalling motion, while those of the disk are consistent with Keplerian rotation. The radius and velocity of the centrifugal barrier are estimated to be about 530 au and 6 km/s, leading to a central mass of about 11 solar masses, consistent with estimates based on spectral energy distribution fitting. These results indicate that an ordered transition from an infalling-rotating envelope to a Keplerian disk through a centrifugal barrier, accompanied by changes of types of molecular line emissions, is a valid description of this massive protostellar source. This implies that at least some massive stars form in a similar way as low-mass stars via Core Accretion.


Introduction
Massive stars impact many areas of astrophysics, yet there is little consensus on how they form. One of the key questions is whether massive protostars accrete through rotationally supported, i.e., Keplerian or near-Keplerian, disks, as have been seen around low-mass protostars. High angular resolution observations, especially of gas kinematics with molecular lines, have provided a handful of candidates of Keplerian disks around massive protostars (Beltrán & de Wit 2016). Most of these disks are around B-type protostars (up to about 15 M e ; e.g., Sánchez-Monge et al. 2013;Beltrán et al. 2014;Ginsburg et al. 2018), with just a few examples reported from around O-type massive stars (e.g., Johnston et al. 2015;Zapata et al. 2015;Ilee et al. 2016;Cesaroni et al. 2017;Maud et al. 2018). More often, especially around O-type protostars, massive (∼100 M e ) rotating "toroids" of radii of 10 3 -10 4 au are found. It is unclear whether such toroids are feeding smaller Keplerian disks at their centers. Searching for Keplerian disks around massive protostars is challenging, due to the far distances and embedded, crowded environments. Often the kinematics of putative Keplerian disks and outer infalling-rotating envelopes are not easily distinguished.
In low-mass star formation, the transition from an infallingrotating envelope to a rotationally supported disk has been much better observed. We note that, in this paper, by the term "disk" we only refer to a rotationally supported disk, i.e., where rotation is the dominant form of support against infall, and we also do not distinguish an infalling-rotating envelope from an infallingrotating "pseudo-disk" (see also discussions in Section 4.3). As the material in the core collapses, the infalling motion is gradually converted to rotation. The innermost radius that such material can reach without losing angular momentum is the centrifugal barrier (Sakai et al. 2014b; see also Stahler et al. 1994), inside of which a rotationally supported disk is expected. The centrifugal barrier is accompanied by not only a change of kinematics from infalling rotation to Keplerian rotation but also a change of chemical compositions due to the accretion shock and different temperature and density conditions in the disk and envelope (e.g., Sakai et al. 2014a;Oya et al. 2015Oya et al. , 2016Oya et al. , 2017. For example, in the low-mass sources studied by these authors, the infalling-rotating envelopes outside of the centrifugal barriers are often traced by molecules such as CCH, c-C 3 H 2 , CS, and OCS. The accretion shocks associated with the centrifugal barriers are often highlighted by volatile species such as SO and saturated organic species. The Keplerian disks can be traced by molecules such as C 18 O and H 2 CO. The observations of different molecules help to disentangle the disk and envelope and highlight the transition region, which provide a powerful diagnostic tool to understand the whole picture of disk formation. There are two main theories for massive star formation: core accretion, which is a scaled-up version of low-mass star formation (e.g., McKee & Tan 2003), and competitive accretion (e.g., Bonnell et al. 2001;Wang et al. 2010), in which stars chaotically gain their mass via the global collapse of the clump without passing through the massive core phase (see Tan et al. 2014 for a review). In the core accretion scenario, a transition from relatively ordered rotating infall of a massive core at scales of ∼10 3 au to a rotationally supported disk on scales of several× 10 2 au or even smaller is expected, and these different components may be highlighted by emissions of different molecules, similar to the case of lowmass star formation, though the particular molecules tracing the various components may be different from the case of lowmass star formation, due to different temperature, density, and shock conditions. In the competitive accretion model, disks are also expected, but they are likely to be much smaller, due to close protostellar interactions. In addition, since the collapse is more disordered, a simple envelope-disk transition, accompanied by clear change of kinematic and chemical patterns, is not expected. Therefore, searching for Keplerian disks and understanding how and where the envelope-disk transition happens by using multiple molecular lines is important to test different theories of massive star formation.
Recent observations have shown indications of a centrifugal barrier between the envelope and disk around the massive protostar G328.2551-0.5321 from the enhancement of highexcitation CH 3 OH emissions (Csengeri et al. 2018). However, so far, there are no simultaneous confirmations of centrifugal barriers around massive protostars using both kinematic and chemical features. In this paper, we report observations by the Atacama Large Millimeter/Submillimeter Array (ALMA) of the massive protostellar source G339.88-1.26, revealing the transition from an infalling-rotating envelope to a Keplerian disk, identified by both kinematic and chemical patterns.

The Target and Observations
2.1. The Target G339. 88-1.26 Our target G339.88-1.26 (aka IRAS 16484-4603; hereafter G339) is a massive protostellar source at a distance of d= 2.1 kpc (Krishnan et al. 2015). Interferometric radio continuum observations revealed an elongated structure of ∼10″ with a position angle of 45°east of north (Ellingsen et al. 1996;Purser et al. 2016), which is believed to be tracing an ionized jet. Also, 6.7 GHz CH 3 OH maser spots are found to be linearly distributed in both space (within a scale of ∼1″) and velocity, with the spatial distribution approximately perpendicular to the radio jet (Ellingsen et al. 1996), which led to an explanation that the masers are tracing a rotating disk. Mid-Infrared (MIR) observations at 10 and 18 μm revealed emissions on a scale of 4″, elongated in about the east-west direction, which is resolved into three peaks (1A, 1B, and 1C;De Buizer et al. 2002). De Buizer et al. (2002) further argued that there are two stellar or protostellar sources present. One source, which corresponds to the MIR emission peak 1B, is an embedded high-mass source driving the radio jet. Another source (about 1″ to the west of 1B, and not corresponding to any MIR emission peak) is a massive star slightly in the foreground and less obscured, which is creating an extended H II region (however, see Section 6.2). In such a scenario, the masers are believed not to trace the accretion disk but to be generated by shocks associated with the embedded source. Such an origin of the maser emissions is also supported by the polarization observations of the CH 3 OH masers (Dodson 2008). Based on spectral energy distribution (SED) fitting using models of single massive protostars (Zhang & Tan 2018), the total luminosity of the source is estimated to be about (4-6)×10 4 L e , and the protostellar mass is estimated to be about 12-16 M e (Liu et al. 2019).

Observations
The observations were carried out with ALMA in Band 6, covering a frequency range from 216 to 232 GHz, on 2016 April 4 with the C36-3 configuration and on 2016 September 15 with the C36-6 configuration. The total integration time is 3.5 and 6.6 minutes in the two configurations. A total of 36 antennas were used in both configurations. The baselines ranged from 15 to 462m in the C36-3 configuration and from 27m to 3.1km in the C36-6 configuration. J1427-4206 was used for bandpass calibration, J1617-5848 and Titan were used for flux calibration, and J1636-4102 and J1706-4600 were used as phase calibrators. The source was observed with single pointings, and the primary beam size (half-power beamwidth) was 22 9.
The data were calibrated and imaged in CASA (McMullin et al. 2007). After pipeline calibration, we performed selfcalibration using the continuum data obtained from a 2 GHz wide spectral window with line-free bandwidth of about 1.6 GHz. We first performed two phase-only self-calibration iterations with solution intervals of 30 and 6 s, and then one iteration of amplitude self-calibration with the solution interval equal to the scan interval. We applied the self-calibration phase and amplitude solutions to the other spectral windows. The peak signal-to-noise ratio of the continuum image is increased by a factor of 2 by the self-calibration. The data of two configurations were then combined after the calibration and self-calibration of the individual data sets. The resultant largest recoverable scale is about 11″. To image the data, the CASA task clean was used, using robust weighting with the robust parameter of 0.5. For the continuum imaging, in addition to the 2 GHz wide spectral window (line-free bandwidth of 1.6 GHz), we also combined the line-free channels of other spectral windows, making a total continuum bandwidth of about 2.3 GHz. The synthesized beam of the 1.3 mm continuum data is 0 29×0 21 with P.A.=28°. 2. The continuum peak position is derived to be a d  Figure 1 shows the 1.3mm continuum emission of G339. This reveals one compact main source at the center with extended emission, as well as two separate sources: one in the northwest (about 11″, i.e., 2.3×10 4 au, from the main source), and one in the north (about 8″, i.e., 1.7×10 4 au, from the main source). The extended emission associated with the main source appears to have substructures. It is elongated mostly in the north-south direction. Slightly to the south and connected to the central source, there is an elongated structure in the northwest-southeast direction, which may be affected by the outflow (see Section 3.2). The source in the northwest is also compact. We identify it as a protostar with outflow activities (see Section 3.2). The source in the north, however, appears to be not so compact and may be part of the extended emission associated with the main source, especially as the even more extended emission is resolved out by the interferometric observation.

1.3 mm Continuum Emission
The total flux of the continuum emission above 3σ within 6″ from the central source is 0.71Jy. To estimate the gas mass from the continuum flux, we adopt the dust opacity of Ossenkopf & Henning (1994;k = - ) and a gas-to-dust mass ratio of 141 (Draine 2011). While a temperature of 30 K is typically used for molecular cores, radiative transfer simulations show that the average temperature of the envelope within ∼10,000 au around a 12-16 M e protostar (see Section 2.1) is ∼70 K . The resultant gas mass is 58 M e assuming a temperature of 30K, or 23 M e assuming 70K. This is most likely to be a lower limit for the gas mass, due to the resolving out of more extended emission, which suggests that there is enough material for future growth of the massive protostar. The total flux of the continuum emission within 0 5 (1000 au), i.e., the immediate compact structure around the central protostar, is 0.16Jy, which corresponds to a gas mass of 13 M e assuming a dust temperature of 30 K, 5 M e with 70 K, or 1.6 M e with 200 K, considering even higher temperatures close to the protostar . The continuum fluxes of the sources in the north and northwest are 0.032 and 0.027Jy, respectively, which correspond to gas masses of 2.2 and 2.6 M e assuming a dust temperature of 30 K, or 0.9 and 1 M e assuming a temperature of 70 K. ). In the integrated emission map (Figure 2(a)), the eastern and western outflows appear to have different opening angles and directions. However, as the 12 CO channel maps ( Figure 13 in Appendix A) show, in low-velocity channels (e.g., = -V 43 lsr and −39 km s 1 ) the eastern and western outflow lobes appear to have similar opening angles and are quite well aligned along a common axis. Therefore, we estimate that the large-scale CO outflow has a half-opening angle of about 20°with a position angle of about 120°east to the north on both sides (indicated by the red dashed lines in Figure 13 in Appendix A). Furthermore, although the eastern outflow is dominated by redshifted emission and the western outflow is dominated by blueshifted emission, there are both blue-and redshifted emissions on both sides at low velocities, which are best seen in the channels of = -V 43 lsr and −23 km s 1 (i.e., » -| | V 10 km s out 1 ). Detections of both red-and blueshifted emissions on both sides suggest a near-edge-on view of this outflow, i.e., the inclination angle between the outflow axis and the plane of sky is likely to be small (similar to or smaller than the half-opening angle, i.e., i20°). An outflow from this source in roughly the east-west direction was not detected before. The direction of this CO outflow is actually similar to that of the elongated MIR emission (De Buizer et al. 2002; see Section 6.2). The direction of this outflow is also similar to that of the CH 3 OH maser distribution, which suggests that these masers may trace the shocks produced in outflow activities (see discussions in Section 6.1). There is a tentative indication of a second molecular outflow from the main source with position angles of 18°(blueshifted) and −135°(redshifted) (labeled with dashed arrows in Figure 2(a)). In the channel maps at low velocities (e.g., = -V 43 lsr and −39 km s 1 ), there are also small-scale emissions close to the main source with a position angle of about 45°. These emissions may belong to the molecular outflow associated with the radio jet previously observed (Ellingsen et al. 1996;Purser et al. 2016), which has a position angle of 46°(indicated by yellow arrows in Figure 2(a)). If this is the case, it indicates that there is an embedded proto-binary system with individual outflows almost perpendicular to each other. However, at current resolution (∼0 3), this binary is still unresolved. For the continuum source located in the northwest of the main source, a small collimated bipolar outflow is detected. There is no clear outflow emission associated with the continuum source north of the main source. Figure 2(b) shows a zoom-in view of the CO outflow emission. The SiO (5−4) emission shows a highly collimated jet extending from the continuum peak to about 5″ away. The jet structure in the SiO only shows emission in the redshifted velocities ranging from = -V 33 lsr to −7 km s 1 (see Figure 14 in Appendix A). The position angle of the SiO jet is about 110°, slightly different from the large-scale 12 CO outflow with a position angle of 120°. However, the redshifted 12 CO emission on the same scale of the SiO jet coincides very well with the SiO emission. This may be caused by the jet precession. The width of the SiO jet is about 0 6, and there is no apparent change in the jet width with distance to the source. 1,1 ) lines. The SiO emission not only traces the jet but also is strongly peaked at the continuum source. The SiO emission peak coincides very well with the continuum peak, with an elongation in the direction of the jet. The SiO emission at the continuum peak is detected within a wide velocity range from = -V 60 lsr to 6 km s 1 (i.e., up to about 30-40 km s 1 relative to the systemic velocity). The spatial distributions of the other molecular line emissions can be categorized into three types. First, the SO 2 and H 2 S emissions are strongly peaked at the positions of the continuum source and the SiO emission peak. The SO 2 emission is only seen within ∼1″ (2100 au) from the central source, while the H 2 S emission also traces the extended structure up to ∼4″ (8400 au) from the central source. The emissions of these two molecules at the continuum peak are detected in a velocity range smaller than the SiO emission (up to about 20 . Their emission peaks are offset from the continuum peak and the SiO emission peak. On the larger scale, they generally follow the morphology of the extended continuum emission. In addition to that, they also show strong emissions associated with the outflow  cavity, which is most clearly seen in the CH 3 OH emission, showing a cone-like structure surrounding the SiO jet. Third, the C 18 O emission is widely spread. Its emission peak is ∼1″ to the north of the continuum peak, corresponding to a separate peak in the extended continuum emission. The large-scale C 18 O emission follows the morphology of the continuum emission, especially to the north. However, it is much more widely distributed than the continuum emission. The C 18 O emission is only detected within a narrow velocity range up to about 5 km s 1 relative to the systemic velocity (see Figure 15 in Appendix A). Figure 4 shows the intensity profiles of the integrated emissions of these lines, as well as the continuum emission, along a cut passing through the continuum peak and perpendicular to the SiO jet (i.e., =  P.A. 20 ). This further confirms the features discussed above. We see that SiO, SO 2 and H 2 S emissions are highly peaked at positions close to the central source, while the peaks of the CH 3 OH and H 2 CO emissions are more offset from the central source (∼0 3). On the other hand, the C 18 O emission is quite widespread across the region. Figure 5 shows the moment 1 maps of SiO (5 − 4), C 18 O (2−1), CH 3 OH (4 2,2 −3 1,2 ; E), H 2 CO (3 2,1 −2 2,0 ), SO 2 (22 2,20 −22 1,21 ), and H 2 S (2 2,0 −2 1,1 ) emissions in a region within 1″ (2100 au) from the central source. Velocity gradients are seen in all of these molecular line emissions across the central source approximately in the north-south direction. The magnitudes of the velocity gradients are different in these molecular emissions. The SiO and SO 2 emissions appear to have the strongest velocity gradients confined in a region close to the central source (<0 5, 1000 au). The H 2 S emission has a strong velocity gradient on similar scales, but also a smaller gradient on larger scales. On the small scale, the velocity gradients in CH 3 OH and H 2 CO are smaller than those in SiO or SO 2 , but velocity gradients can also be seen on larger scales. The velocity gradient in the C 18 O emission is the weakest among these lines. The moment 1 maps of these molecular lines for the whole region are shown in Figure 20 in Appendix B. Besides the varying velocity gradient levels, the detailed directions of the velocity gradients are also different in these molecular lines. The velocity gradient in the SiO emission, as well as those in the SO 2 and H 2 S emissions, is mostly perpendicular to the SiO jet axis ( i.e., P.A.=20°; black arrow in Figure 5), which is consistent with rotation. The velocity gradients in the CH 3 OH and H 2 CO emissions, however, also have components along the direction of the outflow axis, which may indicate infalling motion in the envelope and/or outflow motion, in addition to rotation.

Velocity Structure
To better show the velocity structures, in Figure 6 we show the position-velocity (PV) diagrams of these molecular line emissions, along a cut passing through the continuum peak and perpendicular to the SiO jet (i.e., P.A.=20°). Signatures consistent with rotation across the jet axis are seen in all these molecular lines, with the southern side dominated by the blueshifted emissions and the northern side dominated by the redshifted emissions. Note that, from the PV diagram of the SiO emission (panel (a)), the offsets are symmetric to a position slightly south of the continuum peak position by 0 05 (marked by the horizontal lines in Figure 6), which is much smaller than the resolution beam size (0 3). However, while the SiO emission probes material with high rotation velocities up to ∼30 km s 1 close to the position of the central source, the highest velocities detected in the CH 3 OH and H 2 CO emissions are <10 km s 1 and offset from the central source in space. The kinematics of the SO 2 and H 2 S emissions shown in the PV diagrams appear to be in between those of SiO and those of CH 3 OH and H 2 CO. In velocity space, they show velocities higher than those of CH 3 OH and H 2 CO, but not as high as those of the SiO emission. Spatially, the main parts of the SO 2 and H 2 S emissions are more confined than the CH 3 OH and H 2 CO emissions. Note that the H 2 S emission also contains an extended component with very low velocities. The C 18 O emission only shows a very weak velocity gradient. The positional offsets of the most blue-and redshifted C 18 O emissions are also consistent with those seen in the CH 3 OH and H 2 CO PV diagrams. Figure 7 shows the spectra of these lines within a radius of 1″ from the continuum peak position. The SiO line clearly shows strong high-velocity wings, which are not seen in other molecular lines. The SO 2 and H 2 S lines are narrower than that of SiO, but wider than the other lines. Note that the emission at --V 16 km s lsr 1 in the SO 2 spectrum is due to the CH 3 CHO (11 1,10 −10 1,9 ; E) emission. The CH 3 OH and H 2 CO lines, despite being very bright, are even narrower. The C 18 O line is the narrowest among these lines. The different velocity ranges of these molecular lines are from genuinely different velocity components, as we discussed above, rather than the result of different line intensities.
The different behaviors of these molecular lines in their velocity and spatial distributions suggest that they trace different structures around the protostar. In the CH 3 OH and H 2 CO emissions, the highest velocities are detected offset from the position of the central source. In addition, on the northern side, while most of the emission is redshifted, there is significant blueshifted emission as well. The opposite is true for the southern side. Meanwhile, velocity gradients are also seen in the PV diagram perpendicular to the disk direction (see Figure 8). Such behaviors are consistent with infalling-rotating motion in the envelope (e.g., Sakai et al. 2014a;Oya et al. 2015Oya et al. , 2016Oya et al. , 2017, rather than pure rotation in the disk (see also Appendix C). The locations of the most blue-and redshifted emissions then correspond to the innermost radius of such an envelope (indicated by the blue arrows in Figure 6), where all the kinetic energy is converted to rotation (i.e., the centrifugal barrier; Sakai et al. 2014b). The rotation reaches its maximum velocity at this point without losing angular momentum. Higher rotational velocity is not seen in these lines since they do not trace the disk inside of the envelope. The SiO, SO 2 , and H 2 S emissions, on the other hand, show higher velocities close to the central source, indicating that they are tracing the disk. However, the SO 2 and H 2 S emissions have their highest velocities lower than those of the SiO emission, suggesting that they may only trace the outer part of the disk. Note that the low-velocity components of the SiO, SO 2 , and H 2 S emissions in the PV diagrams are consistent with the CH 3 OH and H 2 CO emissions. This suggests that they trace not only the disk but also some parts of the envelope. The C 18 O line has much lower critical density than the other lines (Table 1), so its emission is dominated by the outer low-density material and does not show high velocities.
To summarize, the spatial distributions and kinematics of these molecular emissions are consistent with a scenario in which there is a transition from an infalling-rotating envelope to a disk, and different molecular emissions trace different components (see also Appendix C). The SiO emission traces the disk and the inner envelope, the CH 3 OH and H 2 CO emissions trace the envelope, and the SO 2 and H 2 S emissions trace the outer part of the disk, as well as the inner envelope. Based on such a scenario, we can approximately estimate the velocity and radius of the centrifugal barrier (marked by the red dashed lines and blue arrows in Figure 6 , where i is the inclination angle defined with i=0°for an edgeon disk and i=90°for a face-on disk. Note that the estimated dynamical mass * m d includes both the protostellar mass and disk mass. The ratio between the disk mass and the protostellar mass f d is uncertain. Theoretical modeling suggests that the disk mass can be a significant fraction of the protostellar mass (f 1 3; d e.g., Kratter et al. 2008), while some observations provide estimates of f d =1/30−1/15 (e.g., Ilee et al. 2016). The estimated dynamical mass is roughly consistent with the mass estimate of 12-16 M e based on SED fitting (Liu et al. 2019). With such an estimated central mass of 11 M e , the red curves in Figure 6(a) show the rotation curve of a Keplerian disk inside of the centrifugal barrier, assuming a negligible disk mass. Note that the rotation velocity of the Keplerian disk at the centrifugal barrier is a factor of 2 lower than the rotation velocity of the infalling-rotating envelope at the centrifugal barrier (see Equation (3)). The high-velocity components of SiO emission are consistent with such a rotation curve, suggesting a rotationally supported disk inside the centrifugal barrier. However, it is difficult to obtain the detailed rotation profile in the disk with the current angular resolution of the data.

Kinematic Model of the Rotating Envelope-Disk
In order to better illustrate the changes of kinematics and types of molecular line emissions in the transition from envelope to disk, we construct simple models to compare with the observed PV diagrams. In order to separate the kinematic feature of an infalling-rotating envelope and a disk, we construct models for these two components separately. In the model, the envelope starts from the centrifugal barrier at a radius of r CB as its inner boundary and extends to an outer boundary r out , with the rotation velocity v j and infall velocity v r described as Such motion conserves both angular momentum and mechanical energy. r CB is the innermost radius that such infalling gas can reach with the angular momentum conserved. At = r r CB , v r =0 and = j v v CB . Since we only focus on the kinematic features in this paper, we adopt simple geometry and density structures in the model. We assume that the envelope has a height of h(r)=0.2r on each side of the midplane, and the density distribution follows r µ -( ) r r 1.5 . For simplicity, we assume that the emissions are optically thin and the excitation conditions are universal across the region. We also fix =  r d 0. 8 out ( =ŕ 1.7 10 au out 3 ) based on the observed PV diagrams. Therefore, we have in total three free parameters for the model of infalling-rotating envelope: the radius of the centrifugal barrier, r ; CB the rotation velocity at the centrifugal barrier, v ; CB and the inclination angle, i. In such a model, the central mass is . For the disk inside of the envelope, we assume that it is Keplerian and it has the centrifugal barrier as its outer boundary. The rotation velocity v j and infall velocity v r in the disk are Here we assumed * * = m m d . We assume that the disk also has a height of h(r)=0.2r on each side of the midplane and a density profile of r µ -( ) r r 2.5 . Similarly, we also assume that the emissions are optically thin and the excitation conditions are universal across the region.
To obtain the best-fit model, we compare the model PV diagrams of the infalling-rotating envelope with the PV diagrams of the CH 3 OH, H 2 CO, and the outer part of the  SiO emissions, and we compare the model PV diagram of the Keplerian disk with the PV diagram of the inner part of the SiO emission. We explore the inclination angle i with values ranging from 0°to 40°with an interval of 10°, the angular radius of centrifugal barrier r CB /d in a range of 0 1-0 4 with an interval of 0 05, and the projected centrifugal barrier velocity v i cos CB in a range of 3-8 km s 1 with an interval of 1 km s 1 and approximately determine the best-fit model by eye (given the approximate nature of the modeling).
The best-fit model has an inclination angle of i=10°b etween the line of sight and the disk midplane, a radius of the centrifugal barrier of r CB /d=0 25, i.e., r CB =530 au, and a velocity at the centrifugal barrier of = v i cos 6 km s CB 1 . The central mass is therefore derived to be about 11 M e . The best model is shown in Figure 8. It shows that the CH 3 OH and H 2 CO emissions indeed can be explained as an infallingrotating envelope, without any high-velocity component associated with the disk rotation. There are extended emissions in the direction of the outflow axis (panels (d) and (f) in Figure 8), which may be affected by the outflow cavity structure, as mentioned in Section 3.3. The high-velocity SiO emission is consistent with a Keplerian disk, while the lowvelocity SiO emission in the outer region is consistent with the infalling-rotating envelope. Especially, the velocity gradient seen along the outflow axis direction across the central source (panel (b)) cannot be explained by pure rotation in the disk.
The SO 2 and H 2 S emissions show some high-velocity components associated with the disk, but the highest velocity is not as high as the SiO emission, suggesting that they only trace the outer part of the disk. Therefore, we construct another disk model with its inner part truncated at a certain radius r in for modeling the SO 2 and H 2 S emissions. We explore the inner radius of the disk r in /d in the range of 0 01-0 15 with an interval of 0 01. The other parameters of this model are the same as the best model obtained above. By comparing the PV diagrams of the model and observation, we estimate r in /d=0 02, i.e., r in =40 au. Meanwhile, there are some extended low-velocity SO 2 and H 2 S emissions that may be associated with the infalling-rotating envelope.
To summarize, the model supports our hypothesis that the different molecular emissions trace different parts of the transition from an infalling-rotating envelope to a Keplerian disk. The CH 3 OH and H 2 CO emissions trace the infallingrotating envelope outside of the centrifugal barrier at a radius of 530 au. The SiO emission also traces the Keplerian disk inside of the centrifugal barrier, in addition to the envelope. The SO 2 and H 2 S emissions trace the centrifugal barrier and the disk inward but outside of a radius of about 40 au. In Appendix C, we discuss the possibility that all these emissions trace different parts of a single rotationally supported disk without infall motion.

The Model Caveats
We note that the model presented above is only a simple example that is designed to be illustrative. It is focused on explaining the kinematic features seen in different molecular lines. In reality, the geometry of the structures, including the density and temperature distributions, is likely to be more complicated than our model assumptions. The observed change of the types of molecular line emissions in the transition from the envelope to the disk is not only a result of the change of chemical composition but also affected by the change of excitation conditions (see Section 6.1). Also, unlike the model in which the molecules are uniformly distributed in the two components (envelope and disk), these molecules are likely to have more complicated abundance distributions, including vertical differentiations.
In the model, we also ignore the motions of the infalling material in the z-direction (i.e., perpendicular to the envelope/ disk midplane). In ballistic models of infall (e.g., Ulrich 1976; Cassen & Moosman 1981;Stahler et al. 1994), material streams land on the midplane at their centrifugal radii. These streams will collide and form a disk structure, with its outer boundary set by the centrifugal radius of the material infalling along the midplane. Material in this disk structure will continue to spiral inward with both rotation and infall until reaching the centrifugal barrier, inside of which a rotationally supported disk forms (Stahler et al. 1994). Therefore, a pseudo-disk dominated by infalling-rotating motion can form between the envelope and the rotationally supported disk (e.g., Lee et al. 2014). In this paper, however, by only considering the motion along the envelope/disk midplane, we do not distinguish between an infalling-rotating pseudo-disk and an infallingrotating envelope. Furthermore, the disk formation process can be strongly affected by the magnetic field, by removing angular momentum via the magnetic braking effect (e.g., Li et al. 2014;Zhao et al. 2016). Such effects are also beyond the scope of our simple model.
In the model, the velocity gradient seen in the PV diagram along the outflow axis (i.e., perpendicular to the disk direction) is caused by the infalling motion of the envelope. It is slightly blueshifted on the eastern side and redshifted on the western side. In such a case, the back side of the envelope is on the east side and the front side of the envelope is on the west side. That is to say, if an outflow is perpendicular to this envelope, it will be blueshifted on the east and redshifted on the west, which is not consistent with the observed jet. However, this misalignment is actually small considering the near edge-on view of both the envelope and the outflow (see Section 3.2). It is possible that there is a small misalignment between the envelope and the disk/jet. Changes of angular momentum direction during the accretion process are possible considering that the core is highly substructured and expected in the collapse of a turbulent core. It is also possible that the direction of the jet may have changed modestly over time, as discussed in Section 3.2. In fact, there are some distinct blueshifted SiO outflow emissions to the east of the central source. which can be seen in Figure 8(b). If that is the most recently launched jet, it is consistent with the inclination of the envelope. . Based on the spatial and velocity distributions of the SiO emission, the SiO jet can be divided into two parts. The first part is the extended structure starting from about 1″ from the central source up to the end of the jet. It is composed of several knots, with the brightest one located at about 1 8 from the central source. There are widespread velocities associated with these knots, which is typical for the SiO emission from shocks along the jet pathway.

The SiO Jet
The second component is within about 1″ from the central source. As discussed in Section 4.1, the SiO emission at the continuum peak position is dominated by the disk and inner envelope, but the SiO emission peak is elongated up to 0 5 from the central source on the east side. As Figure 9 , despite that the rest of the SiO jet is redshifted. As discussed in Section 3.2, the outflow is likely to be close to an edge-on view, i.e., lies in the plane of the sky. In such a case, small changes in jet direction may cause a change from redshifted to blueshifted, i.e., from being inclined away from the observer to being inclined toward the observer. This blueshifted emission may represent the inclination of the most recent jet activity. As discussed in Section 4.2, this inclination is consistent with how the envelope is inclined with respect to the plane of the sky.
Another possibility is that these emissions at 0 3 from the central source trace the material directly launched from the disk, and the wide velocity range is due to the rotation of this disk wind. The SiO molecules in the jet structure can form via two different mechanisms. The first is that they are released to gas by dust sputtering or grain-grain collisions in shocks during the interaction between the jet and the surrounding material, or the interaction between different components in the jet. In such a case, the SiO emission marks the location of strong shocks along the jet pathway. The second mechanism is that the SiO molecules are released to the gas phase in the disk by shocks or strong radiation and then launched to the outflow. In such a case, the SiO emission may trace the motion of the gas that is directly launched from the disk, i.e., a disk wind. The SiO emission tracing the directly launched material from the disk is expected only in regions very close to the central source (Hirota et al. 2017;Lee et al. 2017). In this source, since there is already strong SiO emission associated with the disk, it is natural that some SiO emission also traces the wind launched from the disk. In order to better explore such a possibility, we compare the PV diagrams of the SiO emission along two cuts perpendicular to the jet direction in Figure 10. The distances of these two cuts to the central source are 0″ (along the disk midplane, i.e., same as Figure 6 , which is similar to that of the other parts of the jet. This may suggest that at 0 3 above the disk plane, the material has been accelerated to outflowing velocities. At this height, the PV diagram (Figure 10(b)) also shows a small level of velocity gradient across the jet, with the most blueshifted emission slightly to the south and the most redshifted emission slightly to the north, which is consistent with the rotation direction in the disk (panel (a)). Therefore, it is possible that the widespread velocities at 0 3 above the disk midplane are due to the rotation in the launched material. The maximum rotation velocity detected at 0 3 above the disk midplane is about = -V 33 km s rot 1 , which is higher than the maximum rotation velocity detected along the disk midplane, which is about = -V 27 km s rot 1 . This is possible considering that the specific angular momentum in the launched material is a factor of several times that of material at the launching point on the disk, according to magnetocentrifugal outflow models (e.g., Ferreira et al. 2006) and observations (e.g., ) in lowmass star formation. Transition from disk rotation to outflow rotation in massive protostellar sources has also been inferred in Orion Source I (Hirota et al. 2017). Higher-resolution observations are needed to further confirm such a picture in G339. Figure 11 shows the mass and momentum distributions of the outflow measured from the 12 CO (2−1) emission. To obtain the gas mass, we assume optically thin emission and adopt an abundance of 12 CO of 10 −4 relative to H 2 and a gas mass of 2.34×10 −24 g per H 2 molecule. An excitation temperature of Figure 10. PV diagrams of the SiO (5−4) emission along two cuts perpendicular to the SiO jet direction (P.A.=20°). The distances of these two cuts to the central source are 0″ (i.e., along the disk midplane; panel (a)) and 0 3 (630 au above the disk midplane; panel (b)). The red bar in the lower right corner of each panel indicates the resolution beam size.

Properties of the Large-scale CO Outflow
» -T 10 50 K ex is typically used for deriving the mass from low-J CO transitions (Dunham et al. 2014). Within this range, an excitation temperature of = T 17.5 K ex minimizes the mass estimate from the CO (2−1) line. An excitation temperature of 50K would increase the mass estimate (and therefore momentum estimate) by a factor of 1.5. In each velocity channel, we only include the primary-beam-corrected emissions above 3σ within the region with the primary beam response >0.2. We exclude the emissions that are not related to the main outflow on the east-west direction by applying masks that vary with velocity channels. We also exclude any emissions at velocities -< -| | V V 5 km s lsr sys 1 , due to the confusion of the low-velocity outflow emission with the core material emission. Since the blue-and redshifted emissions appear in both the eastern and western outflow lobes, we also divide the two lobes by a line passing through the central source with a position angle of =  P.A. 30 (i.e., perpendicular to the direction of the large-scale outflow; see Section 3.2) and show the blue-and redshifted masses of the individual lobes. 1 . There are several factors affecting these estimates. First, a different excitation temperature will increase these estimates, e.g., = T 50 K ex would increase the mass estimate (and therefore the other estimates) by a factor of 1.5. Second, the correction factors for inclination on the momentum, dynamical timescale, mass outflow rate, and momentum injection rate are i 1 sin , i i sin cos , i i cos sin , and i i cos sin 2 , respectively, where i is the inclination angle between the outflow axis and plane of sky. Assuming an inclination of i=20°, which is likely to be an upper limit (see Section 3.2), the correction factors for inclination on the momentum, dynamical timescale, mass outflow rate, and momentum injection rate are 3.0, 0.36, 2.7, and 8.0, respectively. Third, the optical depth effect and the missing low-velocity outflow emissions will further increase the mass and momentum estimates. In one example of a low-mass protostellar outflow, this factor is about 10 for mass estimation and about 8 for momentum estimation . Adopting these correction factors for the mass and momentum estimation, and combining the correction factors for the inclination, the total correction factors for the mass, momentum, dynamical timescale, mass outflow rate, and momentum injection rate are 10, 24, 0.29, 34, and 82, respectively, which gives We note that the estimated outflow timescale is likely to be at least an order of magnitude smaller than the formation time of the protostar in the turbulent core model of McKee & Tan (2003), so that only a part of the outflow history is being traced by these observations. The derived mass and momentum rates are consistent with the outflows of other protostellar sources with total luminosities of several´ L 10 4 (e.g., Beuther et al. 2002;Zhang et al. 2005;Maud et al. 2016;Yang et al. 2018). Figure 11(c) shows the mass spectra of the outflow. The mass spectra can be characterized by two power laws, with sudden changes of the spectrum slopes at V out =34 km s 1 in the blueshifted outflow and at = -V 24 km s out 1 in the redshifted outflow. The low-velocity mass spectra have power-law indices of −2.3 and −2.7 for the blue-and redshifted outflows, respectively, and the high-velocity mass spectra have power-law indices of −13.4 and −15.0 for the blue-and redshifted outflows, respectively. The break in the mass spectrum slope around 20-30 km s 1 may be related to molecular dissociation caused by jet shocks (Downes & Cabrit 2007). Such a cutoff at high velocities is also seen in some other massive protostellar outflows, such as the G028.37 +00.07 C1-Sa outflow (Tan et al. 2016), which also has the cutoff at » --V 20 30 km s out 1 . The low-velocity mass spectrum slope index around −2.5 is consistent with other massive outflows (e.g., Maud et al. 2016) and also low-mass outflows (e.g., Richer et al. 2000). The low-velocity mass spectra may be steeper if optical depth effects are considered.

Origin of the Chemical Change across the Centrifugal Barrier
The change of the type of molecular line emissions across the centrifugal barrier is possible from a chemical point of view. The CH 3 OH and H 2 CO molecules are released from icy grain mantles to gas phase in the warm envelope and strongly enhanced around the centrifugal barrier (e.g., Aota et al. 2015). The enhancement around the centrifugal barrier may be due to the accretion shock, but more importantly, it may be due to the broadened inner edge of the infalling-rotating envelope, which can be directly irradiated by the central source (e.g., Sakai et al. 2017). They may be destroyed during the accretion shock at the same time, or reduced as they slowly move inward in the disk. There is also a lack of continuous supply of CH 3 OH and H 2 CO molecules to the gas phase in the disk, since they are mostly formed on the ice mantle of dust grains and have already been released into the gas phase around the centrifugal barrier. Therefore, the CH 3 OH and H 2 CO emissions trace the envelope, have their emission peaks around the centrifugal barrier, and do not show high rotation velocities associated with the disk. In a similar manner, the SO 2 and H 2 S molecules are also enhanced around the centrifugal barrier and gradually reduced as material moves toward the inner part of the disk. However, the observed transitions of these molecules have higher upper energy levels than the CH 3 OH and H 2 CO transitions (see Table 1), so that they tend to trace the region further inward of the CH 3 OH and H 2 CO emissions, i.e., the outer part of the disk inside of the centrifugal barrier. On the other hand, the accretion shock, internal shocks, or a strong radiation field may have destroyed some fraction of the dust grains to liberate SiO, which can remain in the gas phase in the disk (or disk surface). Meanwhile, destruction of dust grains in the disk provides a continuous supply of SiO molecules. Therefore, the SiO emission can reach higher rotational velocities in the inner disk.
As discussed in Section 3.3, the CH 3 OH and H 2 CO emission is also enhanced along the outflow cavity. In fact, Figure 12(a) shows that the locations of the methanol masers in this source coincide very well with the positions of strong thermal methanol emissions, which are to the south of the central source and elongated along the outflow direction. The methanol masers have velocities around = -V 39 lsr to −34 km s 1 , which are also consistent with the elongated structure in the outflow direction shown in the thermal methanol emissions (see Figure 5(c) and Figure 16 in Appendix A). It is possible that the shocks due to the outflow have enhanced the thermal emission and also excite the maser emission of methanol along the outflow cavity wall. Note that the transition region between the envelope and disk is also expected to be the base of outflow cavity. Therefore, in this source, the thermal methanol emissions appear to trace both the envelope (along the midplane) and the outflow cavity walls. Higher angular resolution observations are needed to clearly separate these components.
There is also an asymmetry in the distributions of the CH 3 OH and H 2 CO emissions, with the southern blueshifted emissions much brighter than the northern redshifted emissions. One reason is that the northern redshifted emission is affected by the self-absorption of the foreground infalling material in the envelope, which has a slightly redshifted velocity. This effect is especially strong in the CH 3 OH and H 2 CO emissions, as they only trace the low-velocity components. Similar behaviors have also been seen in lowmass sources for molecular lines tracing the infalling-rotating envelope (e.g., Sakai et al. 2014a;Oya et al. 2016). Another reason is that the shocks may be stronger on the southern side, which enhances the CH 3 OH and H 2 CO emissions. The CH 3 OH maser emissions are also seen on the southern side (see Figure 12(a)), supporting this scenario. The asymmetric shock conditions may be a result of a clumpy distribution of material in the infalling envelope.  Figure 12 compares the previously observed 10 μm, 18 μm, and 9 GHz continuum emissions (De Buizer et al. 2002;Purser et al. 2016) with the ALMA 1.3 mm continuum and molecular line emissions. The 1.3 mm continuum peak does not coincide with any of the three peaks identified in the 10 μm emission (panel (a)), which challenges the claim by De Buizer et al. (2002) that the MIR peak 1B is an embedded protostar driving the radio jet. The 9 GHz radio continuum emissions reveal three knots, with the central knot lying very close to the 1.3 mm continuum source (panel (b)). The ALMA observation shows 12 CO emissions at the positions of the northern and southern radio continuum lobes at velocities close to the systemic velocity (panel (b); see also Figure 13). Other molecular species observed by ALMA, which trace higher densities than 12 CO, do not show emissions concentrated toward these northern and southern lobes, unlike protostellar HC/UC H II regions that are associated with dense molecular structures. These features support a scenario in which the northern and southern radio continuum lobes are actually tracing an outflow from the main source, rather than being from separate massive protostars.
The extension of the MIR emission is in the same direction as the CO outflow (panel (d)), suggesting that the MIR emissions are coming from the outflow and/or outflow cavity. This is consistent with previous observations of other massive protostellar sources and theoretical models, which have shown that the 10-20 μm continuum emissions are strongly affected by the outflow cavity structures in massive young stellar objects (e.g., De Buizer 2006; Zhang et al. 2013Zhang et al. , 2014De Buizer et al. 2017). The driving source of this outflow, however, is not located at the MIR peak 1B as speculated previously. In fact, the 1.3 mm continuum source is located close to the gap seen in the 18 μm emission (panel (c)), which can be naturally explained by the high extinction of the dust concentrated around the central source. The 1.3 mm continuum peak also coincides with the peak of the dust color temperature distribution derived from the 10 and 18 μm emissions (De Buizer et al. 2002), which was believed to indicate a massive star slightly in the foreground and less obscured (see Section 2.1). The ALMA observation suggests that it is the heating from the embedded driving source that is responsible for this temperature peak.
The prediction of MIR emission from outflow cavities is that the blueshifted (near-facing) side should be brighter in the MIR, as the redshifted (far-facing) side is more obscured by the envelope ( we see the opposite in the CO and MIR emissions for G339. The eastern outflow, which is mostly redshifted, has brighter MIR emission than the western outflow, which is mostly blueshifted. One possible reason for this is the almost edge-on view of this outflow. In such a case, the brightness distribution of the MIR emission is not dominated by the overall inclination of the outflow cavity, but strongly affected by the detailed distribution of the extended cold dust in the region. Figure 12(c) shows that cold dust traced by millimeter continuum emission is distributed north, south, and west of the central source. If this dust is in the foreground of the MIR emission, extinction would explain why the MIR peak 1A (where there is no 1.3 mm emission) is so much brighter than source 1C (which is behind the extended millimeter emission). In fact, the MIR emission from 1A appears to be constricted by the extended millimeter emission from the south and north.
To summarize, the ALMA observations of 1.3 mm continuum and molecular line emissions, combined with the MIR and radio continuum observations, support the following scenario for G339. The central source is likely to contain an unresolved massive protostellar system, which may be a binary, but is being fed in a relatively ordered way from the infall envelope. The main source drives an outflow in the east-west direction, seen in 12 CO and MIR continuum. A second outflow may be present in the northeast-southwest direction, seen in radio continuum and also 12 CO emission. The MIR emission is dominated by the main outflow cavities, with its brightness distribution affected by the distribution of the extended cold dust.

Implications for Massive Star Formation
The estimated radius of the centrifugal barrier in G339 (r CB =530±100 au) is similar to the radius of the centrifugal barrier recently identified in the massive protostellar source G328.2551-0.5321, which is 300-800 au (Csengeri et al. 2018). In massive source G17.64+0.16, indication of a change of kinematics from rotation with radial motion to pure rotation is found at a radius of ∼200 au, which may also suggest a centrifugal barrier at that radius (Maud et al. 2018). If the centrifugal barrier is the outer boundary of the disk, the measured radius of the centrifugal barrier in G339 also indicates a smaller disk than most disks or circumstellar structures identified around massive protostars, which typically have sizes of several× 10 3 au (Beltrán & de Wit 2016), although most of these structures are likely to be pseudo-disks, which are expected to extend farther out of the centrifugal barrier, rather than Keplerian disks. Still, the G339 disk is smaller than some recently identified Keplerian disks around massive protostars, e.g., G35. On the other hand, the estimated radius of the centrifugal barrier in G339 is larger than those found around low-mass protostars, which are typically smaller than 200 au (e.g., Sakai et al. 2014b;Oya et al. 2016Oya et al. , 2017Alves et al. 2017). The specific angular momentum measured at the centrifugal barrier in G339 is » r v 3000 au km s CB CB 1 , which is about an order of magnitude higher than those found in low-mass systems (e.g., Sakai et al. 2014b;Oya et al. 2016). The total mass, momentum, and momentum injection rate in the G339 outflow are also significantly higher than those of typical low-mass protostellar outflows (e.g., Dunham et al. 2014;Zhang et al. 2016). However, in spite of these quantitative differences, G339 is found to be qualitatively very similar to lower-mass protostars, including possessing a highly collimated outflow, and ordered transition from an infalling-rotating envelope to a Keplerian disk, which is also accompanied by change of types of molecular line emissions.
As mentioned above, our analysis is based on a simplified model in which the effects of the magnetic field are not considered. It is likely that the magnetic field is playing at least two essential roles in this source. The first is magnetic braking, not only affecting disk formation but also suppressing the fragmentation of the core. Despite the massive envelope containing hundreds of thermal Jeans masses, the lack of fragmentation into small compact sources implies efficient support against collapse and/or transportation of angular momentum by magnetic braking (e.g., Commerçon et al. 2011;Seifried et al. 2011). The second role of the magnetic field is driving the outflow. Recently, both observations (e.g., Hirota et al. 2017) and magnetohydrodynamical (MHD) simulations (e.g., Matsushita et al. 2017;Staff et al. 2018) showed that magnetocentrifugal disk winds can arise in massive star formation in a similar way to that in low-mass star formation. The observed outflow motion in G339 is also consistent with the magnetocentrifugal outflow (see Section 5.1). Higher-resolution observations will provide more detailed information about the processes of magnetic braking and magnetocentrifugal wind launching, as both are most efficient inside the centrifugal barrier, where the magnetic field becomes twisted and tightly wrapped.
Recently, Liu et al. (2019) compared MIR observations of this source with the SED model grid of massive star formation (Zhang & Tan 2018). The model grid is based on the turbulent core model of massive star formation (McKee & Tan 2003) and includes evolutions under various initial and environmental conditions. The best-fit models also show narrow outflow cavities with half-opening angles of 10°-20°, which is consistent with our observations. Among the five best models, there is one model that has a close edge-on view with inclination of i≈20°b etween the outflow axis and the plane of sky. This suggests that the ALMA observations of the outflow can help to break the degeneracies in the SED fitting of infrared fluxes. That model has a protostellar mass of 12 M e , consistent with the dynamical mass estimated from the gas kinematics. This particular model gives a total envelope mass of ∼300 M e within a radius of ∼30″ and an envelope mass of about 17 M e within 10,000 au. The latter is close to the total mass measured from the continuum emission associated with the central source, which is 23 M e assuming a dust temperature of 70K (see Section 3.1). The outflow mass estimated from 12 CO emissions is about 5M e (within ∼15″), including a factor of ∼10 for corrections of optical depth and missing low-velocity emissions (see Section 5.2). If the majority of the 12 CO outflow is entrained material, and if the outflow entrainment is the reason for outflow cavity widening, the mass in the envelope should be about q q -» ( ) cos 1 cos 15 w w times the outflow mass, i.e., ∼75 M e , where we assume a half-opening angle of the outflow of θ w =20°and isotropic distribution of material in the core. This mass is somewhat larger, but still on similar levels as, the measured envelope masses from the 1.3 mm continuum (23 M e ) and SED modeling (17 M e ), considering that the outflow mass is measured on a larger scale (15″) than the 1.3 mm continuum (10,000 au). This result supports a scenario in which the outflow cavity opens up as material is gradually entrained into the outflow (i.e., Arce & Sargent 2006), and it is also expected in the outflow feedback model for massive star formation based on the turbulent core accretion model Tanaka et al. 2017;Staff et al. 2018;Zhang & Tan 2018).
Overall, our results strongly indicate that, at least in this particular example, massive star formation can be considered to be a scaled-up version of low-mass star formation, i.e., forming via relatively ordered core accretion. Furthermore, there is quantitative consistency in many parameters of the system as derived by comparing to semianalytic models based on the turbulent core accretion theory of McKee & Tan (2003). The larger centrifugal barrier radius and specific angular momentum compared to lower-mass sources can be explained by collapse of a rotating massive core on larger scales. On the other hand, in competitive accretion, the fragmentation of the gas clump into many small low-mass interacting cores and protostars, including outflows, is likely to prevent the formation of an ordered disk and rotating infall envelope on these ∼500 au scales. The disk sizes in the competitive accretion model are expected to be much more compact than in core accretion from a massive core (e.g., Kratter & Matzner 2006). The main outflow seen in the 12 CO emission is highly collimated and extends to >15″ (0.15 pc). This suggests that the outflow has not experienced significant disturbance from interactions with other protostars in the region during the past ∼10 4 yr, which is not expected in competitive accretion models. The lack of multiple compact continuum sources in the region, i.e., an apparent relatively low multiplicity, also supports the core accretion scenario rather than the competitive accretion scenario.

Summary
We have presented ALMA observations of the envelope, disk, and outflow system in the massive protostellar source G339.88-1.26 (G339), including 12 CO, SiO, C 18 O, CH 3 OH, H 2 CO, SO 2 , and H 2 S emissions. Our main conclusions are as follows: (1) The 12 CO (2−1) emission reveals a collimated bipolar outflow driven by the G339 protostar. The SiO (5−4) emission also reveals the redshifted jet at the base of the large-scale 12 CO outflow.
(2) The envelope/disk system is traced by SiO (5−4), SO 2 (  -22  22 2,20 1,21 ), H 2 S ( -2 2 2,0 1,1 ), CH 3 OH (  -4  3 ; 2,2 1,2 E), and H 2 CO ( -3 2 2,1 2,0 ) emissions. Based on their spatial distributions and kinematics, we found that these molecular lines trace different parts of the envelope-disk system. The SiO emission traces the disk and inner envelope. The CH 3 OH and H 2 CO emissions trace the infalling-rotating envelope outside of the disk. The SO 2 and H 2 S emissions appear to be enhanced around the transition region between the envelope and disk, i.e., the centrifugal barrier, and trace the outer part of the disk. Therefore, the transition from an envelope to a disk is seen not only in the change of kinematics but also in the change of types of molecular line emissions.
(3) The kinematics of the envelope can be well fit by a model of infalling-rotating motion. In such a model, the envelope collapses with the angular momentum conserved, and the kinetic energy is completely converted to rotation at its inner boundary, i.e., the centrifugal barrier. Based on our model fitting, we estimate the radius of the centrifugal barrier to be about 530±100 au and the rotation velocity at the centrifugal barrier to be about 6±1 km s 1 , leading to a central mass of about -+  M 11 5 6 . Inside of the centrifugal barrier, the rotation appears to be consistent with Keplerian rotation, but higher-resolution observations are needed to confirm it.
(4) We found that the SiO emission slightly above the disk plane (∼0 3) may trace the outflowing material launched from the disk. This emission shows a signature of rotation, which smoothly connects to emission in the disk midplane, indicating angular momentum transfer from the disk to the outflow. (5) We estimate a total mass of 0.46 M e and a total momentum of 5.9 M e km s 1 in the outflow based on the 12 CO emission, without any corrections for inclination or optical depth. After correcting for these effects, we estimate, very roughly, that the total mass, momentum, mass outflow rate, and momentum injection rate are 4.6 M e , 1.4×10 2 M e km s 1 ,´- 1 , respectively. (6) Combined with previous MIR and radio continuum observations, the ALMA observations suggest that the central source drives an outflow in the east-west direction, seen in 12 CO and MIR continuum. A second outflow may be present in the northeast-southwest direction, seen in radio continuum and also 12 CO emission. The MIR emission is dominated by the main outflow cavities, with its brightness distribution affected by the distribution of the extended cold dust. (7) The envelope-disk-outflow system detected in the massive protostellar source G339 appears to be highly ordered and qualitatively very similar to those observed around low-mass protostars, even though the disk size, angular momentum, and outflow strength are much larger than in the lower-mass cases. Our results imply that at least some massive stars form in a way that is similar to those of low-mass stars, i.e., via core accretion. Channel Maps of the Molecular Lines     (  -4  3 ; 2,2 1,2 E), H 2 CO (  -3  2 2,1 2,0 ), SO 2 (  -22  22 2,20 1,21 ), and H 2 S ( -2 2 2,0 1,1 ) emissions on large scales.

Appendix C Kinematic Model Fitting of Rotating Disk without Radial Motion
In Section 4.2, we compared the observed PV diagrams of SiO, CH 3 OH, H 2 CO, SO 2 , and H 2 S emissions with a model composed of an infalling-rotating envelope on the outside and a   Figure 8, but showing models of only the Keplerian disk with pure rotation but no radial motions. All the models (shown in contours) have the same outer boundary but different inner boundaries (see the text for details). The solid contours show the best-fit disk model. The dashed contours show the disk model in the best-fit disk+envelope model presented in Section 4.2 (Figure 8) but extended to outer regions. The left column shows the PV diagrams along a cut perpendicular to the outflow axis (along the disk direction). The positive offsets are to the north of the source. The right column shows the PV diagrams along a cut perpendicular to the disk direction. The positive offsets are on the east side of the source. The model contours are at levels of 0.1, 0.3, 0.5, 0.7, and 0.9 of the peak intensities.