AN ORDERED BIPOLAR OUTFLOW FROM A MASSIVE EARLY-STAGE CORE

Understanding how massive stars form is an important goal, since their radiative, mechanical, and chemical feedback play leading roles in regulating the interstellar medium, star formation activity, and the overall evolution of galaxies. Core accretion (e.g., McKee & Tan 2003, hereafter MT03) is one class of models for massive star formation, which involves the initial conditions of a high-mass, self-gravitating starless core, followed by relatively ordered collapse to a central disk and protostar (see Tan et al. 2014, p. 149 for a review). These models are scaled-up from those developed for low-mass star formation (e.g., Shu et al. 1987), but in the case of the MT03 Turbulent Core Model, they involve non-thermal forms of pressure support, i.e., turbulence and magnetic fields, for the initial core to be in approximate pressure and virial equilibrium.

Alternatively, competitive accretion models (e.g., Bonnell et al. 2001; Wang et al. 2010) involve a massive star gaining most of its mass by competitive, chaotic Bondi–Hoyle accretion in the center of a crowded protocluster of mostly low-mass stars. In these models, the initial conditions of massive star formation, i.e., the gas immediately surrounding the protostar that is destined to become a high-mass star, do not involve massive starless, self-gravitating cores, but rather low-mass cores, with most of the mass reservoir accumulating later from the protocluster clump.

To try and distinguish between these theories we have developed a method for searching for massive starless and early-stage core candidates. We first identify target regions in infrared dark clouds (IRDCs) using mid-infrared, i.e., $8\ \mu {\rm{m}}$ Spitzer-IRAC, extinction (MIREX) mapping (Butler & Tan 2009, 2012). We select regions that are peaks in the resulting mass surface density, Σ, map. We further check that these regions are dark in $24\ \mu {\rm{m}}$ (Spitzer-MIPS) and $70\ \mu {\rm{m}}$ (Herschel-PACS) images. We then search for ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) emission with ALMA, since the abundance of this species, i.e., the deuteration fraction ${D}_{{\rm{frac}}}^{{{\rm{N}}}_{2}{{\rm{H}}}^{+}}$ $\equiv \quad [{{\rm{N}}}_{2}{{\rm{D}}}^{+}]/[{{\rm{N}}}_{2}{{\rm{H}}}^{+}]$, is known to rise in cold ($T\lt 20$ K), dense (${n}_{{\rm{H}}}\gtrsim {10}^{5}\;{{\rm{cm}}}^{-3}$) conditions, especially when CO molecules are largely frozen-out onto dust grain ice mantles and the ortho-to-para ratio of ${{\rm{H}}}_{2}$ drops to low values (e.g., Kong et al. 2015). ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$ is known to be a good tracer of low-mass starless cores that are on the verge of collapse, i.e., pre-stellar cores (Caselli & Ceccarelli 2012), as well as early-stage low-mass Class 0 sources (Emprechtinger et al. 2009).

We carried out a pilot search of 4 IRDC regions (${\rm{C1,}}\;{\rm{F1,}}\;{\rm{F2,}}\;{\rm{G2}}$) with ALMA in Cycle 0 (compact configuration, 23 resolution), identifying $6\ {{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) cores (Tan et al. 2013, hereafter T13) by projection of their $3\sigma \ \ \ l-b-v$ space ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) contours. The two most massive cores were in  the IRDC G028.37+00.07 (hereafter Cloud C) C1 region: we refer to these as C1 north and south (C1-N, C1-S). We estimated masses in two ways: (1) from the MIREX map, finding that C1-N has $61\pm 30\;{M}_{\odot }$ and C1-S has $59\pm 30\;{M}_{\odot }$, with ∼50% systematic uncertainty due to distance ($5\pm 1$ kpc) and dust opacity (∼30%) uncertainties; (2) from millimeter dust continuum emission, finding that C1-N has ${16}_{7}^{33}\;{M}_{\odot }$ and C1-S has ${63}_{27}^{129}\;{M}_{\odot }$, with uncertainties mostly due to the adopted dust temperature of T = 10 ± 3 K, together with distance and dust emissivity uncertainties.

These ALMA observations resolve the cores with about three beam diameters. C1-S appears as quite round, centrally concentrated, and monolithic, while C1-N shows evidence of multiple fragments. Given their high levels of deuteration (Kong et al. 2016) and their dark appearances in Herschel-PACS images, even at wavelengths as long as $100\ \mu {\rm{m}}$, C1-S, and perhaps also C1-N, are among the best known candidates of massive starless or early-stage cores.

However, we note that Wang et al. (2006) reported a water maser detection in this area (just outside C1-S's lowest ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) contour), though at a different velocity (59.5 $\mathrm{km}\;{{\rm{s}}}^{-1}$) and in a single channel (0.66 $\mathrm{km}\;{{\rm{s}}}^{-1}$ wide). This water maser was not detected in the more sensitive observations of Chambers et al. (2009) and Wang et al. (2012). We also note that Pon et al. (2015) have detected CO(8-7) and (9-8) emission toward C1-N and S with Herschel-HIFI ($\sim 20^{\prime\prime}$ resolution) and argue that this emission results from turbulence dissipating in low-velocity shocks, which could be either driven by large-scale turbulent motions from the surrounding cloud or from protostellar outflow activity.

Here we search for potential protostars and outflow activity via ¹²CO(2-1) and other tracers using an ALMA Cycle 2 observation. Below we describe the observations (Section 2), present our results (Section 3), and discuss their implications (Section 4).

We use data from our ALMA Cycle 2 project (2013.1.00248.S, PI:Tan), which observed the C1 region in a compact configuration on 2015 April 05, yielding sensitivity to scales from $\sim 10^{\prime\prime}$ to $\sim 1^{\prime\prime}$. The position of the field center was R.A. = 18:42:46.5856, decl. = −04:04:12.361 (FK5 J2000 system) ($l=28.3230,b=+0.06750$). It was chosen to be between C1-N and C1-S, slightly closer to C1-S. Thus both cores are within the 27'' field of view.

The spectral set-up included a continuum band centered at 231 GHz with a width 1.875 GHz, i.e., from 230.0625 to 231.9375 GHz. The achieved sensitivity was 0.045 mJy per $1\_\_AMP\_\_farcs;51\times 0\_\_AMP\_\_farcs;84$ beam. In this continuum band, each channel has a width 1.129 MHz, i.e., a velocity resolution 1.465 $\mathrm{km}\;{{\rm{s}}}^{-1}$. The ¹²CO(2-1) line frequency is 230.538 GHz. C1-S's radial velocity from its ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) emission is $+79.40\pm 0.01$ $\mathrm{km}\;{{\rm{s}}}^{-1}$, with a one-dimensional (1D) dispersion of 0.365 $\mathrm{km}\;{{\rm{s}}}^{-1}$ (i.e., FWHM = 0.860 $\mathrm{km}\;{{\rm{s}}}^{-1}$), so the sky frequency of ¹²CO(2-1) from this source is 230.477 GHz. Thus we are sensitive to the presence of ¹²CO(2-1) emission with moderate velocity resolution. Note that ambient CO molecules in the core and even the wider scale IRDC are expected to be largely frozen-out onto dust grain ice mantles (Hernandez et al. 2011). However, ¹²CO(2-1) emission near ambient velocities is still likely to be very optically thick from this region. In our data the nine channels closest to the expected ambient velocity of C1-S have velocities centered at: 74.1, 75.4, 76.7, 78.0, 79.3, 80.6, 81.9, 83.2, 84.5 $\mathrm{km}\;{{\rm{s}}}^{-1}$.

The other basebands were tuned to observe ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2), ${{\rm{C}}}^{18}{\rm{O}}$(2-1), DCN(3-2), DCO⁺(3-2), SiO(v = 1)(5-4) and ${\mathrm{CH}}_{3}\mathrm{OH}\ v\ t=0\ 5(1,4)-4(2,2)$. These data will be presented and analyzed in full in a future paper, while in this Letter we focus mostly on the results of the broad continuum band and its ¹²CO(2-1) line, along with some results from ${{\rm{C}}}^{18}{\rm{O}}$(2-1) and DCN(3-2), that help probe denser gas of protostellar cores.

Figure 1(a) presents the MIREX Σ map (Butler et al. 2014) of the C1 region, together with contours of ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) integrated intensity (T13), which define C1-S and C1-N. Also shown are the locations of potential protostars defined by 1.3 mm dust continuum peaks detected in the ALMA Cycle 2 image (Figure 1(b)), which has been cleaned with natural weighting and had a primary beam correction applied. Figure 1(b) also shows integrated intensities of (continuum subtracted) ¹²CO(2-1) with red contours tracing ${v}_{{\rm{LSR}}}=85.8$ to 124.8 $\mathrm{km}\;{{\rm{s}}}^{-1}$, i.e., redshifted velocities up to 45.4 $\mathrm{km}\;{{\rm{s}}}^{-1}\;$with respect to C1-S's central velocity, and blue contours tracing from 33.8 to 72.8 $\mathrm{km}\;{{\rm{s}}}^{-1}$, i.e., blueshifted velocities up to 45.6 $\mathrm{km}\;{{\rm{s}}}^{-1}$ from C1-S. Figure 1(c) shows just high-velocity outflow gas that is $\gt 20$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ away from C1-S's ambient velocity. Figure 1(d) shows the integrated intensity of "ambient" ¹²CO(2-1), i.e., from 75.4 to 83.2 $\mathrm{km}\;{{\rm{s}}}^{-1}$. Figures 1(e) and (f) show integrated intensities of ${{\rm{C}}}^{18}{\rm{O}}$(2-1) and DCN(3-2), respectively, both potentially helpful for identifying dense gas associated with protostellar cores (e.g., Parise et al. 2009), and thus their radial velocities.

Zoom In Zoom Out Reset image size

From these images we clearly identify two protostellar sources that are spatially overlapped with the C1-S core (Table 1). We refer to the more central source as C1-Sa and define its spatial position as the location of the 1.3 mm continuum peak. This peak, with a flux density of 11 mJy/beam, is 131 from C1-S's center as defined by ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) (T13). Recall that for an estimated (kinematic) source distance of $5\pm 1$ kpc, 1'' corresponds to 5000 au, i.e., 0.024 pc, with ∼20% uncertainties. So the spatial location of C1-Sa is quite close to center of the C1-S core, which has a radius of 361 (i.e., 0.0875 pc; 18,000 au).

Table 1.  Cores and Protostars in the C1 Region

Name	l(°)	b(°)	${S}_{1.3\mathrm{mm}}\$(mJy)	${v}_{{\rm{LSR}}}\$ ($\mathrm{km}\;{{\rm{s}}}^{-1}$)	P.A.${}_{{\rm{outflow}}}\$(°)
C1-N	28.32503	0.06724	6.94 ± 0.72	81.18 ± 0.03a	...
C1-S	28.32190	0.06745	26.7 ± 0.77	79.40 ± 0.01a	...
C1-Sa	28.322093	0.067698664	19.5 ± 0.1	79.01 ± 0.12b	155
C1-Sb	28.321752	0.066847223	2.7 ± 0.1	81.36 ± 0.42c	113
C1a	28.324765	0.069543149	3.6 ± 0.1	80.2d	150
C1b	28.323272	0.069987301	3.5 ± 0.1	...	150

Notes.

^a ${v}_{{\rm{LSR}}}$ of C1-N and C1-S estimated from ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) (T13). ${v}_{{\rm{LSR}}}$ of protostars estimated from strongest ${{\rm{C}}}^{18}{\rm{O}}$(2-1) peak, but see the notes below for individual sources. ^bC1-Sa: ${v}_{{\rm{LSR}}}$ estimated from the strongest ${{\rm{C}}}^{18}{\rm{O}}$(2-1) peak. Secondary peaks at 75.08 ± 0.05 $\mathrm{km}\;{{\rm{s}}}^{-1}$, 81.08 ± 0.03 $\mathrm{km}\;{{\rm{s}}}^{-1}$, while DCN(3-2) has a single peak at 79.8 ± 0.2 $\mathrm{km}\;{{\rm{s}}}^{-1}$. ^cC1-Sb: ${v}_{{\rm{LSR}}}$ estimated from the strongest ${{\rm{C}}}^{18}{\rm{O}}$(2-1) peak. Secondary peak at 78.16 ± 0.05 $\mathrm{km}\;{{\rm{s}}}^{-1}$. DCN(3-2) is too weak to measure ${v}_{{\rm{LSR}}}$. ^dC1a: v_LSR tentatively estimated from weak N₂D⁺(3-2) emission (T13).

Download table as:  ASCIITypeset image

Figures 1 (g), (h), and (i) show the spectra of ${}^{12}$CO(2-1), ${{\rm{C}}}^{18}{\rm{O}}$(2-1) and DCN(3-2), also in comparison with the T13 observation of ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2), toward the protostars. We estimate the radial velocity of the C1-Sa protostar from the ${{\rm{C}}}^{18}{\rm{O}}$(2-1) and DCN(3-2) spectra toward the continuum peak. The ${{\rm{C}}}^{18}{\rm{O}}$(2-1) spectrum shows a main peak at +79.01 ± 0.12 $\mathrm{km}\;{{\rm{s}}}^{-1}$, while DCN(3-2) shows a single peak at +79.8 ± 0.2 $\mathrm{km}\;{{\rm{s}}}^{-1}$. Thus it seems very likely that C1-Sa is forming inside the C1-S ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) core, which has a mean velocity of +79.4 $\mathrm{km}\;{{\rm{s}}}^{-1}\;$and an FWHM of 0.86 $\mathrm{km}\;{{\rm{s}}}^{-1}$.

The Galactic coordinate frame position angle of C1-Sa's ¹²CO(2-1) outflow axis, which we define to the blueshifted axis, is ≃155°. The outflow is highly collimated and is seen to extend $\sim 12^{\prime\prime}$ (60,000 au, 0.3 pc), and is quite symmetric, i.e., the P.A. of the redshifted lobe is almost 180° greater than that of the blueshifted lobe. Also, the observed extent of the outflow is similar in each direction, although the highest velocity flow is more extended on the redshifted side. The outflow can be traced down to 3σ above the noise level, without bunching of the contours, so we expect the observed extent to simply be due to observational sensitivity and the actual extent could be much larger.

We re-checked our ALMA Cycle 0 data, which included a requested bandpass set to an intermediate frequency between DCN(3-2) and SiO($v\quad =$ 0)(5-4). However, this frequency was later mistakenly shifted to be closer to DCN(3-2), causing the SiO line center to be unobserved: only the potential blue wing up to ${v}_{{\rm{LSR}}}=+69$ $\mathrm{km}\;{{\rm{s}}}^{-1}\;$was observed. For this reason, T13 did not report the detection of any SiO emission toward C1-S. However, now we do see indications of the blue wing of SiO(v = 0)(5-4) overlapping with the central part of the blue lobe of the ¹²CO(2-1) outflow and extending to ${v}_{{\rm{LSR}}}\simeq +50$ $\mathrm{km}\;{{\rm{s}}}^{-1}$. We conclude that it is likely that the outflow from C1-Sa also emits strongly in SiO(v = 0)(5-4) across its full velocity range.

The second source, C1-Sb, has a much weaker 1.3 mm continuum flux of 1.3 mJy/beam and is located 20 from the center of the C1-S ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) core and 33 from C1-Sa. The ${{\rm{C}}}^{18}{\rm{O}}$(2-1) spectrum toward C1-Sb shows a main peak at $+81.36\pm 0.42$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ and a secondary peak (with about half the equivalent width) at +78.16 ± 0.45 $\mathrm{km}\;{{\rm{s}}}^{-1}$. The DCN(3-2) spectrum shows no particularly strong features, although a 3σ peak is seen in the integrated intensity map (Figure 1(f)). We tentatively assign C1-Sb's radial velocity to be that of the main ${{\rm{C}}}^{18}{\rm{O}}$(2-1) spectral feature. We discuss below that this assignment is potentially supported by examination of the channel maps of the ¹²CO(2-1) "ambient" gas. If this radial velocity is correct, then it would suggest that C1-Sb is not physically associated with the C1-S ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) core, and in fact may be part of a gas structure that is linked to the C1-N core. However, we cannot exclude the possibility that C1-Sb is also forming from the C1-S core.

The outflow from C1-Sb has a similar extent as that from C1-Sa and appears to have a wider opening angle. It has a P.A. = 113°, and on its blueshifted side the outflow spatially overlaps with that from C1-Sa. Analysis of lower intensity contours (down to 3σ) indicates a bunching close to the higher-intensity contours. This may indicate that, unlike for C1-Sa, we are seeing the full extent of C1-Sb's outflow.

In Figure 1(b) we identify two more candidate protostars that are away from the C1-S and C1-N cores. C1a is located in a region with faint ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) emission with a radial velocity of 80.6 $\mathrm{km}\;{{\rm{s}}}^{-1}$. There is relatively faint ¹²CO(2-1) emission, which may be driven from this source with P.A. ∼ 130°. There is no significant ${{\rm{C}}}^{18}{\rm{O}}$(2-1) emission associated with this source, and only a weak $3\sigma$ feature in the DCN(3-2) integrated intensity map. C1b is close to C1a and has a similar P.A. as its ¹²CO(2-1) emission of about 130°. There is no significant ${{\rm{C}}}^{18}{\rm{O}}$(2-1) or ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) emission at C1b, and only a very weak feature in DCN(3-2). Note that the protostellar nature of C1a and C1b is uncertain, especially since some of the observed ¹²CO(2-1) features may be affected by side-lobe contamination from C1-Sa and C1-Sb.

Figure 2 shows nine channel maps of "ambient" ¹²CO(2-1) emission from the C1 region, with high-velocity and low-velocity ends connecting with the "outflow" velocities plotted in Figure 1(b). The C1-Sa outflow lobes are visible in the low-velocity and high-velocity channels. Figure 2 also suggests that the C1-Sb outflow has a driving source with a ${v}_{{\rm{LSR}}}\simeq +81$ $\mathrm{km}\;{{\rm{s}}}^{-1}$, since the blueshifted lobe is already apparent in the 79.3 $\mathrm{km}\;{{\rm{s}}}^{-1}\;$channel, while the redshifted lobe appears to vanish by the 81.9 $\mathrm{km}\;{{\rm{s}}}^{-1}\;$channel.

Zoom In Zoom Out Reset image size

There are several other striking features seen in Figure 2. First, there is a very elongated "filament" that peaks in the 76.7 and 78.0 $\mathrm{km}\;{{\rm{s}}}^{-1}$ channels, but is visible from 74.1 to 79.3 $\mathrm{km}\;{{\rm{s}}}^{-1}$. The filament overlaps with the C1-Sa, C1-Sb, and C1b protostars and its orientation is almost perpendicular to their outflows. An interpretation of this filament as an ambient gas feature, rather than as a collimated bipolar outflow, is discussed below, considering its position–velocity diagram. Second, there is a relatively weak, but still highly significant, quasi-spherical "core" of gas seen in the 80.6 and 81.9 $\mathrm{km}\;{{\rm{s}}}^{-1}$ channels. Third, there is ¹²CO(2-1) emission in the vicinity of the C1-N ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) core. Fourth, there are additional emission features on the periphery of the image.

Figure 3 shows position–velocity diagrams of ¹²CO(2-1) emission along the outflow axes of C1-Sa and C1-Sb and along the axis of the "ambient filament." These are defined by rectangular regions of 3'' wide running length-wise along the P.A. of each outflow or the filament.

Zoom In Zoom Out Reset image size

The ambient nature of the filament is readily apparent in these figures. It has a very narrow velocity dispersion, which is comparable to other ambient gas tracers, like ${{\rm{C}}}^{18}{\rm{O}}$(2-1) and ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2), and does not show a significant gradient in radial velocity. The lack of such a gradient could be explained by an outflow that was precisely aligned in the plane of the sky, but this would still be expected to have a relatively broad velocity dispersion, which the filament feature does not exhibit.

Figure 4 shows ¹²CO(2-1) spectra and derived mass and momentum distributions of the C1-Sa and C1-Sb outflows, extracted from the same regions used for Figure 3. To derive the mass, one must assume an excitation temperature, with values of ${T}_{{\rm{ex}}}\simeq 10$−50 K typically being used. ${T}_{{\rm{ex}}}=17\;{\rm{K}}$ minimizes the mass estimate, while 50 K increases this by a factor of 1.5.

Zoom In Zoom Out Reset image size

For C1-Sa, integrating from $50\lt ({v}_{{\rm{LSR}}}/\mathrm{km}\ {{\rm{s}}}^{-1})\lt 140$ and assuming ${T}_{{\rm{ex}}}=17$ K, the blue/red lobes have masses ${m}_{w}=0.50/0.32\;{M}_{\odot }$ and momenta ${p}_{w}=3.6/6.6\;{M}_{\odot }$ $\mathrm{km}\;{{\rm{s}}}^{-1}$, respectively. Similarly, for C1-Sb, integrating from $45\lt ({v}_{{\rm{LSR}}}/\mathrm{km}\;{{\rm{s}}}^{-1})\lt 110$ the blue/red lobes have masses $0.73/0.38\;{M}_{\odot }$ and momenta $8.5/3.5\;{M}_{\odot }$ $\mathrm{km}\;{{\rm{s}}}^{-1}$, respectively. Note that the estimates for each blueshifted lobe are affected by the overlap of the sources, leading to a modest overestimation of their properties, especially relative to the redshifted lobes. However, in general the absolute values of the above estimates should be viewed as lower limits, not only because of the choice of ${T}_{{\rm{ex}}}=17$ K, but also because of inclination effects (which boost momentum estimates by $1/{\rm{cos}}(i)$, where i is the inclination of the outflow axis to the line of sight, with a random expectation value of i = 60°), and because of optical depth effects (both within the outflowing gas, which may boost momentum by factors of ∼6 (Zhang et al. 2016), and due to foreground absorption).

For C1-Sa, the mass-weighted mean velocities (${v}_{w}={p}_{w}/{m}_{w}$) of the blue/red lobes are ${v}_{w}=7.3/21\$ $\mathrm{km}\;{{\rm{s}}}^{-1}$, while for C1-Sb they are 12/9.2 $\mathrm{km}\;{{\rm{s}}}^{-1}$. Assuming that the lengths of all the outflows lobes are $\simeq 12^{\prime\prime}$ (60,000 au), the dynamical times for the blue/red lobes of C1-Sa are ${t}_{w}=(3.9/1.3)\times {10}^{4}$ year and for C1-Sb are $(2.4/3.1)\times {10}^{4}$ year. The correction factor for inclination is cos(i)/sin(i), i.e., 0.58 for i = 60°. The mass outflow rates for the C1-Sa blue/red lobes are $(1.3/2.5)\times {10}^{-5}\;{M}_{\odot }\;{{\rm{yr}}}^{-1}$, and for C1-Sb are $(3.0/1.2)\times {10}^{-5}\;{M}_{\odot }{{\rm{yr}}}^{-1}$. The correction factor for inclination is sin(i)/cos(i), i.e., 1.73 for i = 60°. Finally, the momentum injection rates for the C1-Sa blue/red lobes are ${\dot{p}}_{w}=(0.9/5.0)\times {10}^{-4}\;{M}_{\odot }$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ $\;{{\rm{yr}}}^{-1}$ and for the C1-Sb blue/red lobes are $(3.5/1.1)\times {10}^{-4}\;{M}_{\odot }$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ $\;{{\rm{yr}}}^{-1}$. The inclination correction factor is sin(i)/${{\rm{cos}}}^{2}i$, i.e., 3.46 for $i=60^\circ$. Allowing for a mass underestimation factor of 3 and assuming $i=60^\circ$, the overall estimates of the total momentum fluxes are boosted by a factor of 10, i.e., totals for C1-Sa of ${\dot{p}}_{w}\sim 5.9\times {10}^{-3}\;{M}_{\odot }$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ yr⁻¹ and for C1-Sb of $\sim 4.6\times {10}^{-3}\;{M}_{\odot }$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ yr⁻¹.

Outflow momentum flux is expected to be the most reliable direct tracer of protostellar properties, since it should be independent of the effects of ambient environment (unlike mass flux, which depends mostly on the mass that has been swept-up by the primary outflow). Models of massive protostar formation (Zhang et al. 2014) based on the Turbulent Core Model (MT03) for cores of $60\;{M}_{\odot }$ in a clump environment with ${{\rm{\Sigma }}}_{{\rm{cl}}}\simeq 0.4\;{\rm{g}}\;{\mathrm{cm}}^{-2}$ (relevant for C1-S) predict ${\dot{p}}_{w}\lesssim 6\times {10}^{-3}\;{M}_{\odot }\;\mathrm{km}\;{{\rm{s}}}^{-1}\;{\mathrm{yr}}^{-1}$, when the protostellar mass is ${m}_{*}\lesssim 3\;{M}_{\odot }$, rising to $\sim {10}^{-2}\;{M}_{\odot }\;\mathrm{km}\;{{\rm{s}}}^{-1}\;{\mathrm{yr}}^{-1}$ by the time ${m}_{*}\simeq 10\;{M}_{\odot }$. It is interesting that these estimates are comparable with the observed values of ${\dot{p}}_{w}$ from C1-Sa and C1-Sb. This suggests that if C1-Sa is a massive protostar in the process of formation, it is currently at a very early stage, i.e., it has yet to accrete most of its mass. This would be broadly consistent with the protostar having a relatively low luminosity such that it does not appear yet as an MIR source. High angular resolution, high-sensitivity MIR to FIR observations, e.g., with JWST, are needed to measure the SED of the protostar, which can then also constrain protostellar models. We conclude that we have detected protostars of relatively low current masses. C1-Sa appears to be embedded within the C1-S massive, cold core, as defined by ${{\rm{N}}}_{2}{{\rm{D}}}^{+}$(3-2) emission. It thus has a large mass reservoir from which to continue to grow: we speculate that it is destined to become a massive star.

As traced by ¹²CO(2-1), C1-Sa's bipolar outflow is highly collimated and has velocities extending to $\sim 50$ $\mathrm{km}\;{{\rm{s}}}^{-1}$. Similar (blueshifted) velocities are seen in SiO(v = 0)(5-4) emission. Using mass-weighted mean velocities, which are ∼10 $\mathrm{km}\;{{\rm{s}}}^{-1}$, the i = 60° inclination-corrected outflow timescale is $\sim 2\times {10}^{4}$ year. However, since C1-Sa's outflow is likely to extend to larger distances than what we observe, this is probably a lower limit on the duration of protostellar activity. Given the symmetric and linear morphology of the outflow lobes, it appears that C1-Sa (and C1-Sb) have not suffered significant dynamical disturbance from other nearby (proto)stars during the period they have been driving these outflows. This is consistent with the assumptions of core accretion models and is a constraint on competitive accretion models.

For constant instantaneous star formation efficiency from the core, ${\epsilon }_{c}$, the fiducial turbulent core model predicts ${m}_{*}={\epsilon }_{c}{M}_{c}{(t/{t}_{*f})}^{2}$, where the total star formation time is ${t}_{*f}\to 1.29\times {10}^{5}{({M}_{c}/60{M}_{\odot })}^{1/4}{({{\rm{\Sigma }}}_{{\rm{cl}}}/1{\rm{g}}{\mathrm{cm}}^{-2})}^{-3/4}$ year. So for the ${M}_{c}=60\;{M}_{\odot }$ and ${{\rm{\Sigma }}}_{{\rm{cl}}}=0.4\;{\rm{g}}\;{\mathrm{cm}}^{-2}$ case, then $t\gt 2\times {10}^{4}$ year implies ${m}_{*}\gt 0.36\;{M}_{\odot }$ (assuming ${\epsilon }_{c}\simeq 1$, expected during early stages when outflow cavity opening angles are small) and ${\dot{p}}_{w}\gtrsim 1\times {10}^{-4}\;{M}_{\odot }\;\mathrm{km}\ {{\rm{s}}}^{-1}\;{\mathrm{yr}}^{-1}$.

We conclude that C1-Sa is a good candidate for an early-stage massive protostar and as such it shows that these early phases of massive star formation can involve highly ordered outflow, and thus accretion, processes (see also Zhang et al. 2015). The massive C1-S core is potentially also forming a second protostar, C1-Sb: improved estimates of its radial velocity are needed to clarify its association with this core. The results presented complement work that has shown that collimated outflows can be launched from later-stage, more luminous massive protostars (e.g., Beuther et al. 2002). They also illustrate that there are similarities between high-mass protostars and their lower-mass cousins forming from lower-mass cores.

We thank Ke Wang and an anonymous referee for comments and discussions that improved the paper. J.C.T./S.K. acknowledge an NRAO/SOS grant and an NSF grant AST1411527. P.C. acknowledges the financial support of the European Research Council (ERC; project PALs 320620). This paper uses ALMA data: ADS/JAO.ALMA#2013.1.00248.S; ADS/JAO.ALMA#2011.0.00236.S. ALMA is a partnership of ESO, (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory (NRAO) is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.