Multidirectional Mass Accretion and Collimated Outflows on Scales of 100–2000 au in Early Stages of High-mass Protostars

As part of the ALMA cycle 3 program 2015.1.01596.S, we observed two fields centered on W51north (α(J2000) = 19^h23^m4005, $\delta ({\rm{J}}2000)=+14^\circ 31^{\prime} 05$5) and W51e2/e8 (α(J2000) = 19^h23^m4391, $\delta ({\rm{J}}2000)=+14^\circ 30^{\prime} 34$6) with long baselines in Band 6 (216–237 GHz). The project was carried out in three executions on 2015 October 27 and 30. Between 37 and 42 antennas in the 12 m array were employed and provided baselines ranging from 85 to 16,196 m. The precipitable water vapor was between 0.85 and 1.98 mm for the observations, and there was reasonable phase stability (see below). We employed the Band 6 sideband-separating receivers in dual polarization mode to record nine spectral windows (SPWs). Seven SPWs had an ∼234 MHz bandwidth, were recorded with 960 channels, and were Hanning smoothed by a factor of 2, achieving a frequency resolution of 564 kHz (corresponding to ∼0.75 km s⁻¹ at Band 6). This setup enabled us to cover a large number of spectral lines from different molecular species. The remaining two SPWs had a broader bandwidth of 1.8 GHz and a frequency resolution of 1.13 MHz (corresponding to ∼1.5 km s⁻¹) to obtain sensitive continuum measurements at 217 and 235 GHz, as well as to cover additional spectral lines. We spent 2.2 hr on-source during a total observing time of 5.3 hr.

The data calibration and imaging was carried out in the Common Astronomy Software Applications (casa) package (version 4.5.1) and followed standard procedures. Fast referencing was used for these long-baseline observations; the phase calibrator was observed every ∼70 s. The phase calibrator–to–target field separation angle was small at <12. The combination of fast switching times and a close phase calibrator is imperative for accurate phase calibration interpolation for these long baselines. Additionally, the water vapor radiometer (WVR) scaling algorithm (Maud et al. 2017) was implemented to improve the short-term (3–6 s) phase stability. This is a stand-alone modular package that runs inside the casa environment and provides a user with the optimal value to scale the water vapor corrections in the WVRGCAL task within casa. The phase calibration of the data was then further improved by phase self-calibrating on the continuum emission down to an integration timescale of 50 s (using the casa tasks tclean and gaincal). The combined effect of WVR scaling and phase self-calibration improved the dynamic range of the final continuum maps by about 35%. The self-calibration solutions from the continuum were also applied to the full spectral data set (using the casa task applycal) in order to improve the signal-to-noise ratio of the line cubes.

The calibrated visibilities were transformed from the Fourier plane (uv domain) into the image domain using the tclean algorithm. Three sets of images were produced using uniform, Briggs (robust parameter R = 0.5), and natural weighting schemes that achieved angular resolutions of ∼19, 28, and 35 mas, respectively. Since interferometric observations lead to spatial filtering of large-scale structures, the maximum recoverable scale in our images is 04 (about 2000 au at a distance of 5.4 kpc). This implies that we do not detect any emission that is extended over scales larger than 2000 au. In our images, we only include data from baselines longer than 300 m to mitigate striping in the images from large-scale emission detected on the shortest baseline that could not be properly imaged (baselines with a length below this value being predominantly in one direction).

2.2.1. Continuum Emission

2.2.2. Line Emission

Figure 2 shows the 1.3 mm continuum emission, in units of brightness temperature (T_B ), tracing warm dust in the W51e2e, W51e8, and W51north hot cores. The dust emission is resolved into complex structures and displays multiple components on scales of 100–3000 au. In the following, we estimate sizes (Section 3.1.1) and masses (Section 3.1.2) of different components, as well as luminosities (Section 3.1.3) of the three dusty sources.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.1. Sizes

3.1.2. Masses

For the mass estimates, we make a distinction between optically thick and optically thin emission.

We first assume that in the highest column density regions (i.e., toward the source centers), the sources are just barely optically thick (τ = 1). Therefore, N(H₂) = τ/κ_ν  = 1/κ_ν , where κ_ν is the dust opacity coefficient. We use the same coefficient, κ_{227 GHz} = 0.0083 cm² g⁻¹, as in Ginsburg et al. (2017), resulting in N(H₂) = 2.6 × 10²⁵ cm⁻². The mass in a beam is then M_beam = N(H₂) × A = A/κ, where A is the beam area (0033 × 0024) in cm², which yields dust masses of each source M_beam = 0.36 M_⊙.

Although the highest column density peak positions have increased optical depth, the emission surrounding the peak position is likely to be optically thin. We can then use the flux density measurements to compute gas masses using the following equation:

$$\begin{eqnarray*}&&M=\displaystyle \frac{{S}_{\nu }{d}^{2}}{{\kappa }_{\nu }{B}_{\nu }({T}_{d})},\end{eqnarray*}$$

where S_ν is the source flux density, d is the distance to the source, κ_ν is the dust opacity coefficient (which includes the gas-to-dust ratio of 100), and B_ν (T_d ) is the Planck function for a blackbody at dust temperature T_d . The latter is not known, but we can estimate lower limits to T_d in two different ways. The first estimate is provided by the peak brightness temperature of the millimeter continuum emission at the source center (where we assume that the emission is optically thick). The peak intensities of the three sources in the robust 0.5 maps are S_{227 GHz}(e2e) = 19.4, S_{227 GHz}(e8) = 14.7, and S_{227 GHz}(north) = 13.2 mJy beam⁻¹, which correspond to brightness temperatures of ${T}_{B,\max }\,=$ 575, 435, and 390 K, respectively. These provide the first lower limits on T_d . A second estimate comes from LTE modeling of CH₃OH lines imaged with the ALMA cycle 2 program (Ginsburg et al. 2017), which provided temperatures in the range 200–600 K inside 5000 au. Therefore, we can assume $T={T}_{{\rm{gas}},\min }=200$ K as a second lower limit on T_d .

In the following, we use these two approaches (assuming thin and thick emission) and independent T_d lower limits to obtain mass estimates of the spatially resolved cores.

W51e2e.—The central core object spans 27 beam areas (radius r ≈ 500 au), implying a mass ${M}_{{\rm{e}}{\rm{2}}{\rm{e}},\mathrm{core}}=0.36\times 27=9.7$ M_⊙ if it is entirely optically thick. If it is optically thin, on average, using the expression above for ${T}_{d}={T}_{B,\max }=575$ K, the total is 4 M_⊙, while for ${T}_{d}={T}_{\mathrm{gas},\min }=200$ K, the total is 12 M_⊙. In the surrounding filamentary structures, we estimate the mass assuming T_d  = 200 K and obtain M_{e2e,filaments} = 25 M_⊙ (excluding the inner core).

W51e8.—For the filament extending to 1200 au from the central peak along the NW, we measure a mass of 4 M_⊙ for ${T}_{d}={T}_{B,\max }=435$ K and 9 M_⊙ for ${T}_{d}={T}_{\mathrm{gas},\min }=200$ K. For extended structure toward the SW, we obtain M(200 K) = 16 M_⊙.

W51north.—The central peanut-shaped structure occupies 14 beams, implying a lower-limit mass ${M}_{\mathrm{north},\mathrm{core}}=0.36\times 14=5$ M_⊙ if it is optically thick or M = 4.2 M_⊙ if it is optically thin and isothermal at ${T}_{d}={T}_{B,\max }=390$ K. For the larger surrounding structure, we obtain a mass range M(390 K) = 14 M_⊙ and M(200 K) = 28 M_⊙.

Table 1 summarizes the sizes and masses of the different components: centrally peaked core, surrounding compact core, and filaments.

Table 1. Properties of the 1.3 mm Dusty Sources Associated with HMYSOs W51e2e, W51e8, and W51north

Source		Central Beam		Compact Core		Extended Structure		Dust		Protostellar
		Radius	Mass	Radius	Mass	Radius	Mass	TB Peak		Luminosity a
		(au)	(M⊙)	(au)	(M⊙)	(au)	(M⊙)	(K)		(L⊙)
W51e2e b	⋯	70	0.36	500	4–12	2700	25	575		>8.2 × 103
W51e8 c	⋯	75	0.36	⋯	⋯	2500	16	435		>7.5 × 103
W51north d	⋯	⋯	⋯	350	4–5	1200	14–28	390		>2.8 × 103

Notes. The sizes and masses are estimated for different components of the dust continuum emission.

^a These protostellar luminosities are lower limits, as inferred in Section 3.1.3. ^b In W51e2e, the mass estimate in the extended filamentary structure does not include the inner core. ^c In W51e8, there is no core-like symmetric structure surrounding the bright central source. ^d W51north does not contain a central bright point source.

Download table as:  ASCIITypeset image

We compare the recovered flux in the long-baseline maps to that in the ALMA cycle 2 maps at the same frequency (Ginsburg et al. 2017). In the robust 0.5 maps, in W51e2e, 13% of the flux is recovered; in W51e8, 12%; and in W51north, 23%. These low recovery fractions indicate that most of the dust emission is smooth on  ≳2000 au scales. It is therefore likely that the cores reside in a larger envelope of dense material containing up to 5–10 times as much mass within r ≲ 2000 au as we have estimated above (Ginsburg et al. 2017).

3.1.3. Protostellar Luminosities

Each of the three target sources is a typical hot molecular core and extremely rich in spectral lines (see Ginsburg et al. 2017). Within the plethora of identified spectral lines, we identified a subset of lines that are unblended and strong enough to be used for kinematic analysis. The molecular lines observed directly toward the dusty structures are only seen in absorption against the continuum, with the exception of SiO, which is seen in the outflow (see Section 3.3).

In order to study gas kinematics, we used intensity-weighted velocity maps and velocity dispersion maps of different dense gas tracers, which achieve a typical spatial resolution of about 150–200 au. In Figure 4, we display velocity field maps for selected lines of CH₃CN ⁸ for W51e2e, W51e8, and W51north. A remarkable feature seen in these velocity maps is the nearly perfectly flat velocity field in all three sources near their systemic velocity (shown in green). This property is best seen toward W51north, where different CH₃CN lines from the K-ladder (with lower energy levels between 122 and 515 K) show exactly the same remarkably flat velocity profile regardless of their excitation (Figure 4, bottom row). This result does not depend on the molecular tracers or probed scales and thus mirrors real physical properties of the targeted HMYSOs. The only exception is W51e2e (Figure 4, top left panel), which shows redshifted velocities along the NW dust lane (see Section 4.1 for an interpretation).

Zoom In Zoom Out Reset image size

In Appendix A, we review possible explanations for the lack of velocity structure toward these cores, and we conclude that it is an observational effect caused by the high optical depth in the dust continuum, as predicted in Krumholz et al. (2007).

We detected bright SiO (J = 5–4 v = 0 line) emission that traces compact bipolar outflows stemming from the three HMYSOs. These small-scale SiO outflows are barely detected beyond r ∼ 2000 au, but they connect to larger-scale outflows (2000–30,000 au) detected in ¹²CO (J = 2–1) ⁹ by Ginsburg et al. (2017). The outflows from W51e2e and W51north are shown in Figures 5 and 6 (see Appendix B for a detailed description of their structures).

Zoom In Zoom Out Reset image size

Both CO and SiO outflows display high speeds (with peak velocities ±100 km s⁻¹). Such high speeds coupled with the compact sizes imply that these outflows are dynamically very young. In particular, the CO outflows have a dynamical age of a few thousand years, whereas the SiO outflows have a dynamical age of roughly 100 yr (see Appendix C.1 for dynamical age estimates).

While the W51e2e outflow has a simple bipolar morphology, the W51north outflow shows a sharp change from small to large scales. The inner SiO outflow is simple and bipolar at a position angle (P.A.) of −70° out to r < 1300 au, at which point the redshifted flow is sharply truncated (the blueshifted lobe extends out to r ≲ 1600 au). Starting at about 800 au from the center, the redshifted SiO emission suddenly turns by 80° (at a P.A. of ∼10°) and continues to the north out to ∼5000 au, where it meets the large-scale ¹²CO outflow, which extends out to ∼20,000 au at a P.A. of −13°. The outflow is also detected in H₂O maser emission, which consistently traces the SiO morphology and kinematics (Imai et al. 2002; Figure 6, lower panel).

Zoom In Zoom Out Reset image size

We discuss possible scenarios to explain this peculiar outflow structure in Section 4.4.

We can use the SiO emission structure to estimate the momentum and ejection rates of the outflows. One straightforward method would be to determine the mass of the outflowing gas from SiO emission and divide it by the outflow dynamical age. While this may provide the correct order-of-magnitude answer in some cases, molecular outflows, such as those seen in CO and SiO emission, are driven by an underlying primary wind/jet and probe mostly entrained ambient gas moving at lower velocity (this applies especially to CO).

Since in protostellar jet/outflow systems, momentum is conserved (we assume momentum-driven outflows; e.g., Masson & Chernin 1993), one can set a tight constraint on the mass-loss rate using $\dot{{m}_{j}}{v}_{j}=\dot{{m}_{o}}{v}_{o}={\dot{P}}_{o}$, where $\dot{m}$ is the mass-loss rate, v is the speed, and $\dot{P}$ is the momentum rate (j and o stand for the primary jet and molecular outflow, respectively).

The outflow momentum rate can be obtained by dividing the outflow momentum by its age. In Appendix C.1, we estimate the dynamical age of the outflows from the ratio of their projected length to their maximum speed by assuming an inclination i = 45°. In Appendix C.2, we estimate the momentum of the outflowing gas by computing its mass from the SiO emission integrated in the full velocity range. In order to convert from the SiO column to the H₂ column, we need to fix the abundance of SiO, X_SiO = N(SiO)/N(H₂). Toward a sample of 14 high-mass star-forming regions, Leurini et al. (2014) used the APEX single-dish telescope to target the SiO (J = 3–2 and J = 8–7) lines, obtaining X_SiO = (0.7–4.8) × 10⁻⁸, while Sánchez-Monge et al. (2013b) used the IRAM 30 m telescope to target the same sample in the SiO (J = 2–1 and J = 5–4) lines, obtaining X_SiO = (0.1–3) × 10⁻⁸. We assume a conservatively high abundance of SiO X_SiO = 10⁻⁷ (a lower abundance of SiO would imply higher mass and momentum in the outflow). The ages, momenta, and resulting momentum rates derived with this analysis are reported in columns 5, 7, and 8 of Table 2, respectively. The momentum analysis was conducted on the blueshifted lobes (which display a simpler structure; see Appendices B and C.2); therefore, the total outflow mass and momentum rate are found by multiplying those numbers by a factor of 2 (to account for both outflow lobes), yielding M_o  = 0.72, 0.36, and 0.22 M_⊙ and ${\dot{P}}_{o}=0.42,0.24\,\mathrm{and}\,0.06$ M_⊙ yr⁻¹ km s⁻¹ for the outflows driven by W51e2e, W51north, and W51e8, respectively. Using the ¹²CO J = 3–2 line, Shi et al. (2010b) found a larger outflow mass of 1.3 M_⊙ and a lower momentum rate of 0.04 M_⊙ yr⁻¹ km s⁻¹ for W51e2 (both estimates use only the blueshifted lobe). Their larger mass is not surprising, since ¹²CO probes swept-up ambient gas on much larger scales (≳2'' or  ≳10,000 au versus ≲2600 au), while a lower momentum rate is expected for a much older fossil outflow (with an age of about 1600 yr versus 115 yr for the SiO outflow).

Table 2. Outflow Parameters Derived from the SiO (J = 5–4) Line Emission in W51north, W51e2e, and W51e8

Source	ΔV	Vmax	Rmax	Tdyn	Iint	M	P	$\dot{P}$	${\dot{M}}_{\mathrm{ejection}}$
	(km s−1)	(km s−1)	(au)	(yr)	(K km s−1)	(M⊙)	(M⊙ km s−1)	(M⊙ yr−1 km s−1)	(M⊙ yr−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
W51north	[−36, 56]	95	1540	77	1.39 × 104	0.18	9.2	0.12	2.4  × 10−4
W51e2e	[−56, 53]	105	2560	116	1.34 × 104	0.36	24.8	0.21	4.3  × 10−4
W51e8	[0, 60]	42	390	44	0.94 × 104	0.11	1.5	0.03	0.7  × 10−4

Note. The parameters are estimated from the blueshifted lobe. The total outflow mass and rate are twice the quoted values. To convert from the SiO column to the H₂ column, we assume a conservatively high abundance of SiO, X_SiO = N(SiO)/N(H₂) = 10⁻⁷. A lower X_SiO would imply higher mass and momentum in the outflow. The age estimates assume an inclination angle of 45°. For an inclination i, the values need to be multiplied by a factor $\cos (45^\circ )/\cos (i)$.

Download table as:  ASCIITypeset image

At this point, we can convert the computed momentum rates of the molecular outflows into the ejection rates of the primary winds/jets. This requires knowledge of their speed, which, however, we do not measure. Jet speeds measured in low-mass outflows via spectroscopy of atomic lines in the optical and near-IR wavelengths are typically a few hundred kilometers per second (e.g., Frank et al. 2014). For massive outflows, there are fewer measurements of jet velocities, but a handful of proper-motion studies of the radio continuum jets provide speeds of about 500 km s⁻¹ (e.g., HH 80/81; Marti et al. 1995). Using the total momentum rates computed above, we derive the following mass ejection rates for the primary jet: $\dot{{m}_{j}}=\left[\displaystyle \frac{{v}_{j}}{500\,(\mathrm{km}\,{{\rm{s}}}^{-1})}\right]\times \left[\tfrac{\cos 45^\circ }{\cos (i)}\right]\times 8.6\times {10}^{-4}$ M_⊙ yr⁻¹ in W51e2e, 4.8 × 10⁻⁴ M_⊙ yr⁻¹ in W51north, and 1.4 × 10⁻⁴ M_⊙ yr⁻¹ in W51e8, where i is expressed in units of degrees and v_j is the jet speed in kilometers per second. We emphasize that the estimate for W51e8 is less reliable owing to its high apparent inclination (we detect red- and blueshifted emission on both sides of the central protostar, indicating i < 45°).

The dust emission does not show a simple flattened structure at the center of the outflows, as one would expect if there is a large accretion disk (see Ginsburg et al. 2018); instead, the emission is resolved into multiple "streamers" or "dust lanes" that converge onto the central cores. These dust lanes appear elongated in several directions, both perpendicular and parallel, relative to the SiO outflows in all three sources. They typically extend across 01–04 (∼500–2000 au) from the compact sources; have relatively high surface brightnesses, indicating high intrinsic temperatures (≳50–150 K; see Figures 2 and 3); and carry a significant amount of mass (several solar masses; see Section 3.1.2).

Lack of kinematic signatures of accretion.—Having identified these dusty streamers, we attempted to determine their kinematics to assess their nature, in particular to distinguish between material that is plunging into the cores (infall) or being flown out (outflow). Unfortunately, we do not identify any clear kinematic signatures toward them (Section 3.2 and Figure 4). The only exception is W51e2e (Figure 4, left panel), which shows redshifted velocities along the NW dust lane. Since these lines are seen in absorption against the dust continuum, one would be prompted to conclude that they probe dense gas infalling toward the compact core. However, the redshifted velocities are observed away from the dust continuum peak and only along the redshifted lobe of the SiO outflow, which is pointing away from us, suggesting that this material may actually be entrained by the outflow. In this scenario, the lower opacity expected along the outflow cavity may explain why dense gas tracers like CH₃CN display kinematics only along that specific line of sight.

Are the dust lanes outflow rather than accretion flows?—The presence of CH₃CN along the outflow raises the question of whether the dust lanes are actually probing warm dust heated by the outflows (e.g., by shocks), rather than accretion filaments. There are, however, at least three counterarguments to this scenario. First, in W51e2e, there are at least four dust lanes with similar sizes and brightnesses, which appear elongated in several random directions relative to the SiO outflows (especially when considering we see a 2D projection from a 3D structure); therefore, it is unlikely that they all identify dust in the outflow cavities. Similar arguments apply to W51north, where there are three different lanes with no common point of symmetry. Second, if the outflow were responsible for the observed dust emission, the lack of bright emission toward the blueshifted lobe in the SE from W51e2e would imply a much lower density of material in that direction, yet the blue lobe truncates sharply at a distance of about 2000 au at a point where no hot dust emission is observed. The hot dust emission is therefore almost certainly driven by (1) the presence of (much) more gas+dust toward the NW than the SE and (2) the proximity of that gas+dust to the central star. Even if the outflow is interacting with this material, it is clear that there is a major asymmetry in the total mass. Third, we suspect that heated dust in an outflow cavity would be too weak to be seen at the distance of W51, while the measured high surface brightness along the filaments indicates a high intrinsic temperature.

Since these features are not part of an outflow, we suggest that they are accretion filaments feeding those central sources. Indirect evidence for infall comes from the fact that these dust lanes have a symmetric structure around protostellar sources with very high accretion rates (see Section 4.3). More indirect evidence comes from a previously published magnetic field study with ALMA, as we detail in the next paragraph.

Magnetic fields and filamentary accretion.—The magnetic field and gravity are the two main forces that influence the dynamics of dense cores where star formation takes place. We compare our observations with the polarization maps of Koch et al. (2018), who conducted a multiscale analysis of magnetic fields in W51. The plane-of-the-sky component of the magnetic fields inferred from the polarized dust emission exhibits a different morphology at angular resolutions of 2'', 07, and 026 for W51e2e, W51e8, and W51north, respectively (see Figure 2 in Koch et al. 2018). At a resolution of 2'', the magnetic field in W51e2e displays a nearly uniform distribution in the east–west direction. The higher-resolution images reveal a radial distribution of the field centered at W51e2e. The magnetic field orientations align well with the major axes of the streamers presented in Figure 1. The alignment of the dust filaments with the radially distributed magnetic fields within the W51e2e core supports a scenario where the core collapse leads to accretion flows along the filaments that drag the magnetic fields, rending a radial distribution of the field lines. The dust emission toward the W51e8 and W51north cores is less filamentary compared to the W51e2e core. Nevertheless, the magnetic field maps in these two cores also exhibit radial distributions that are suggestive of the field being dragged by the mass inflow. This scenario has been proposed to explain magnetic field properties in dust polarization studies of high- and low-mass star formation (e.g., Zhang et al. 2014; Hull & Zhang 2019; see also Section 5). Future observations of kinematics in the W51 cores are needed to confirm the accretion flows in these regions.

The compact cores at the centers of the dusty streamers (i.e., toward the dust continuum peaks) are where one would expect to see kinematic signatures of rotation (e.g., velocity gradients) if large accretion disks were present. Although we cannot directly estimate their actual size, we can measure upper limits on the potentially disk-containing regions based on the observed extent of the optically thick dust emission (Section 3.1.1). This yields diameters of 150 au for W51e8, 700 au for W51north, and 1000 au for W51e2e (indicated with black arrows in Figure 1).

The putative compact disks are surrounded by a morphologically flattened envelope of size of order of 1000–2000 au, where the filaments appear to converge. Zapata et al. (2009, 2010) used the SMA to reveal a large rotating structure on scales  ≳10,000 au. Although the maximum recovery scale in our study is only 04, which prevents us from recovering the large scales probed with the SMA, no such large-scale rotation was identified in the cycle 2 study with a maximum recovery scale of 9'' (Ginsburg et al. 2017, Figures 31 and 32).

Despite the lack of (large) disks, the presence of fast (∼±100 km s⁻¹), collimated bipolar outflows in SiO emission provides a clear indication of ongoing accretion onto the three HMYSOs. Lacking a clear signature of accretion, we cannot directly estimate the mass accretion rate. In jet–protostar systems, however, ejection rates are expected to correlate with accretion/infall rates. Therefore, we can use the ejection rates calculated in Section 3.4 as a proxy for the infall/accretion rates. Here we use the term "infall rate" to describe gas in the core infalling onto the disk/protostar system and the term "accretion rate" to describe gas accreting onto the protostar itself. Therefore, we can define ${\dot{M}}_{\mathrm{infall}}={\dot{M}}_{\mathrm{ejection}}+{\dot{M}}_{\mathrm{accretion}}$ (assuming that the fraction of the infalling gas that is not accreted is actually reejected via protostellar winds/jets; e.g., Frank et al. 2014). We define ${f}_{j}={\dot{M}}_{\mathrm{ejection}}/{\dot{M}}_{\mathrm{accretion}}$ as the fraction of accreting gas that is launched in the jet. Therefore, ${\dot{M}}_{\mathrm{infall}}=(1+{f}_{j})/{f}_{j}\times {\dot{M}}_{\mathrm{ejection}}$. This fraction is observationally uncertain, but models predict f_j  = 0.2–0.5 (e.g., Offner & Arce 2014; Kuiper et al. 2016). Adopting these values implies ${\dot{M}}_{\mathrm{infall}}\sim 3\mbox{--}6\ {\dot{M}}_{\mathrm{ejection}}$ and ${\dot{M}}_{\mathrm{accretion}}\sim 2\mbox{--}5\ {\dot{M}}_{\mathrm{ejection}}$. Using these assumed fractions, we infer infall rates of (2.6–5.2) × 10⁻³ and (1.4–2.9) × 10⁻³ M_⊙ yr⁻¹ and accretion rates of (1.7–4.3) × 10⁻³ and (1.0–2.4) × 10⁻³ M_⊙ yr⁻¹ onto the W51e2e and W51north protostars, respectively. ¹⁰

In Section 3.3, we show that the W51north outflow appears to have changed direction substantially, as measured with a change of P.A. (projected onto the plane of the sky) from ∼−10° to  ∼−70° from scales ≳2000–20,000 au down to the smallest scales, ≲200–1000 au. Here we argue that this "multicomponent" morphological structure is best explained as a single protostellar outflow that has changed direction over a short period of time (but see possible alternative scenarios in Section 4.4.1). In this scenario, the prominent and compact SE–NW (P.A. ∼ −70°) faster component represents the recently launched outflow, while the faint and diffuse north–south (P.A. ∼ −10°) slower component represents a fossil outflow. In particular, assuming the outflow was driven at a constant velocity v_outflow = 30 km s⁻¹ (the average velocity of the bulk of the SiO emission), we infer the following history of the outflow.

• 1.  
  "Old" flow. From about 22,000 to 1900 au, as traced by CO (and from about 4300, as traced by SiO), the outflow was consistently pointed at a P.A.  of −13° (3500–300 yr ago).
• 2.  
  "Transition" flow. From 1900 to 800 au, the outflow changes directions from a P.A. of −13° to  −70° (300–125 yr ago, duration 175 yr).
• 3.  
  "Young" flow. From 800 au to the resolution limit, the outflow has been consistently pointed along  −70° (∼125 yr ago to now). This outflow history is sketched in Figure 7 for the SiO portion of the outflow. (See also Appendix C.1 for details on the dynamical age estimates.)

The outflow morphology fits this scenario; there is emission detected in both the red and blue lobes at intermediate P.A.s (−10° > P.A. > −70°). This emission is not as cleanly symmetric as the more linear "old" and "young" outflows, which may be caused either by asymmetries in the surrounding medium or by changing the speed of the outflow launch. A point in favor of this scenario is the lack of detected outflow emission in either SiO or CO (Ginsburg et al. 2017) beyond 1300 au along the current SE–NW outflow axis. If the W51north outflow has not changed direction in the last ∼100 yr, the lack of observed outflow tracers at greater distances along the SE–NW implies that accretion onto this source only began in the last 100 yr. Given its high luminosity (L ≳ 10⁴ L_⊙), such an age is implausible.

Zoom In Zoom Out Reset image size

The dramatic change of direction of the W51north outflow suggests a sudden major event that redirected the driving jet over the course of just 150–200 yr. This, in turn, indicates a substantial change in the orientation of the outflow-launching region, i.e., the accretion disk, over that same timescale. The large change in angle (≳60°) suggests that small perturbations or instabilities in the disk are inadequate. Given the presence of a substantial reservoir of surrounding material distributed asymmetrically around the source, we argue that accretion onto the disk is a likely cause for the change of outflow direction. Under the assumption that the present (young) outflow is the result of a single major accretion event and that the infall rate remained constant over the transition period (175 yr), from the infall rate estimated in Section 4.3, we infer that a mass of about 0.25–0.5 M_⊙ was dumped onto the protostar/disk system in such accretion event.

Is this mass enough to flip the disk? We assume that the disk has a radius of 350 au and a mass of 1 M_⊙ (i.e., a quarter of the compact core mass is in the disk; see Table 1) and is in Keplerian rotation about a 20 M_⊙ central star. The total angular momentum in such a disk would be L_disk ∼ 3 × 10⁵⁴ g cm² s⁻¹ (see calculation in Appendix D). If a 0.25–0.5 M_⊙ condensation impacts the disk at a radius of 350 au with a velocity of 10 km s⁻¹ (this corresponds to the infall velocity at 350 au toward a 20 M_⊙ star), the angular momentum of the impactor is ${L}_{\mathrm{dump}}\sim (2.6\mbox{--}5.3)\times {10}^{54}\ \sin \alpha$ g cm² s⁻¹, where α is the angle between the impactor trajectory and the disk plane. For α = 90°, the disk angular momentum vector is reoriented by $\theta ={\tan }^{-1}({L}_{\mathrm{dump}}/{L}_{\mathrm{disk}})=41^\circ \mbox{--}61$° (for a smaller α, the flip would be smaller by a factor $\sin \alpha$). We explicitly notice that θ is independent of the protostellar mass, since both L_dump and ${L}_{\mathrm{disk}}\propto \sqrt{{M}_{* }}$. We conclude that the accretion-driven disk reorientation model is plausible.

4.4.1. Alternative Scenarios

A prominent bipolar outflow, which was first detected in the CO (3–2) line with the SMA by Shi et al. (2010b) and then imaged with ALMA in the CO (2–1) line at 02 resolution by Ginsburg et al. (2017; Figure 5, left panel), is driven by W51e2e. This outflow has a high relative velocity, v ± 100 km s⁻¹. On the largest scales traced by ¹²CO, both ends (red- and blueshifted) of the outflow are sharply truncated at ∼25 (0.07 pc or 15,000 au). To the southeast (SE), the high-velocity flow lies along a line that is consistent with the extrapolation from the NW flow (i.e., directly bipolar); however, at lower velocities (10 km s⁻¹ < V_LSR < 45 km s⁻¹), the direction appears more to the south before being abruptly truncated. For the NW redshifted part of this outflow (70 km s⁻¹ < V_LSR < 120 km s⁻¹), there is apparently a clear bend at an intersection with a blueshifted flow (22 km s⁻¹ < V_LSR < 45 km s⁻¹) from another source, W51e2nw, that is also believed to be an intermediate-to-high-mass protostar (Goddi et al. 2016). The significant change in direction suggests that these two outflows interact. Although such a scenario seems implausible given the outflows' small volume filling factor, both the morphology and the distinct velocity range of the two outflows (which make it easy to distinguish the two) point to a genuine encounter. Therefore, at least for the NW flow, this encounter is responsible for the truncation of the redshifted lobe.

The ALMA long-baseline images unveil the smallest scales of this outflow (Figure 5, right panel). Although the SiO structure is broadly consistent with the lower-resolution ¹²CO structure at the base of the outflow, the P.A. of the line connecting the NW and SE flows is −56° (whereas the ¹²CO outflow is closer to P.A. = −37°). The blueshifted flow extends only across ∼045 or 2430 au. The bulk emission of the redshifted lobe in the NW has a similarly compact size (∼06 or 3240 au), but some weaker emission extends across ∼10,000 au and is cospatial with the ¹²CO outflow (this more extended, weaker structure is recovered in lower-resolution images of SiO tapered to 01, not displayed here). Interestingly, the redshifted lobe shows a two-stream limb-brightened morphology with a central cavity, possibly due to the presence of a dusty filament at the center (see left panel of Figure 1). The dynamical age of the SiO outflow is ∼116 yr (see Appendix C.1), shorter than the 600 yr estimated by Ginsburg et al. (2017) at the peak observed velocity of ¹²CO.

Sato et al. (2010) used the VLBA to measure proper motions of water masers, which also show a picture consistent with the thermal ¹²CO and SiO emission. Interestingly, the blueshifted and redshifted water masers closer to the protostar imply a flow P.A. of −56°, perfectly consistent with thermal SiO, while the redshifted masers along the outflow at larger distances have a P.A. of −43°, closer to that of the ¹²CO outflow axis spatially distant from the protostar. Considered together, these findings are consistent with an indication of a change in the P.A. of the outflow with increasing distance from the protostar.

The outflow from W51north is remarkably extended and complex. Figure 6 offers a full view of this outflow on scales from 25,000 au down to hundreds of au. The first notable feature is that the redshifted and blueshifted lobes on the large scales as traced by ¹²CO are significantly asymmetric. A collimated high-velocity component (covered velocities of ±100 km s⁻¹) extends directly north across ∼4'' (or 20,000 au), while the blueshifted component points to the SE and is sharply truncated, extending only across ∼06 (or 3200 au; Figure 6, left panel; Ginsburg et al. 2017). At the smallest scales, i.e., within a few thousand au from the driving source (Figure 6, right panel), further surprises are revealed. At the base of the ¹²CO outflow, SiO reveals a very fast and well-collimated bipolar outflow, with the blue lobe extending across 034 (or 1600 au) and having a maximum speed of 96 km s⁻¹. The red lobe is slightly more compact, with a radius of 024 or 1300 au and a maximum speed of 66 km s⁻¹. Both SiO lobes are consistently oriented at P.A. ∼ −70°, clearly rotated from the large-scale ¹²CO, which has a P.A. of −13° (Figure 6, left panel).

Close to the central protostar, the redshifted outflow lobe shows a remarkable structure. It has an orientation consistent with that of the blueshifted lobe at P.A. ∼ −70°, up to a radius of 024 (1300 au). Starting at about 015 (800 au) from the center, it suddenly turns by 80° (P.A. ∼ 10°) and curves until it settles at a P.A. of 0° along the base of the redshifted lobe of the ¹²CO outflow. The high-velocity portions of the flow (up to ∼70 km s⁻¹ from the systemic velocity) are found closer to the protostar, while the material in the "transition" region has a lower velocity (65 km s⁻¹ < V_LSR < 90 km s⁻¹, i.e., within 30 km s⁻¹ from the stellar systemic velocity). Some high-velocity material (up to 50 km s⁻¹ from the systemic velocity) is also observed at the northern end of the SiO emission (see Appendix B.3 and Figure 9 for a description of the velocity structure of the outflows). The large-scale red lobe shows contiguous structure tracing onto the small scale, with the change in direction starting at around 035 (1900 au) and completing by about 015 (800 au). Assuming a constant velocity of 30 km s⁻¹, the outflow turn occurred over just 175 yr. The general structure, including the sharp turn seen in the redshifted lobe, is reflected in the blueshifted lobe, which shows a southward "protrusion" starting at about 022 (1200 au) from the protostar. Although much less prominent, this southern component in the blueshifted SiO lobe provides a (weaker and more compact) counterpart to the redshifted northern component. In both SiO and CO emission, the blueshifted lobe is barely detected beyond about 030 or 1600 au (although the weakest CO emission extends up to 3200 au); one possibility is that it breaks out of the cloud (assuming the source is located at the front of it).

Zoom In Zoom Out Reset image size

The peculiar structure of the outflow revealed by ALMA images of SiO emission is also tracked by the H₂O maser 3D velocities (Imai et al. 2002), shown as colored arrows in the right panel of Figure 6. The H₂O masers reside in two complexes separated by ∼3000 au at the front edges of the compact SiO outflow lobes ¹³ and are moving away from each other in ballistic motions at about 200 km s⁻¹ and at a P.A. of −72°, consistent in both P.A. and velocity with the SiO compact outflow. In particular, the H₂O masers detected in the SE appear to trace a bow shock at the tip of the blueshifted SiO lobe, while the H₂O maser complex on the opposite side traces both the high-velocity NW redshifted lobe and the base of the northern outflow component. In the redshifted part of the outflow, the measured proper motions appear to smoothly rotate from NW to north moving away from the tip of the SiO lobe, closely following the SiO emission. In addition, the expansion velocity of the H₂O maser outflow decreases from 90 to 50 km s⁻¹ at radii 02–05 from the outflow origin, also consistent with the SiO redshifted emission. Finally, a least-squares fitting analysis to the measured maser 3D motions provides evidence of a triaxial symmetry, inconsistent with a simple bipolar symmetric outflow but consistent with multiple outflows or a single outflow with a more complex structure.

Figures 8 and 9 showcase the structure of the SiO protostellar outflows as a function of velocity within 3000 and 4000 au from the central protostars W51e2e and W51north, respectively. Note that in W51e2e, the high velocities preferentially trace the central parts (i.e., primary wind), and the low velocities preferentially trace the limb-brightened edges (i.e., entrained gas) of the outflow, a phenomenon often observed in low-mass protostellar outflows (Frank et al. 2014). In W51north, the high velocities trace material at larger radii from the protostar both in the compact NW–SE outflow and the larger-scale north–south outflow, reminiscent of the "Hubble flows" observed in low-mass protostellar outflows. It is also interesting to note that the bulk of the compact NW–SE outflow expands at high velocities, while the bulk of the larger-scale north–south outflow expands at low velocities, in agreement with the scenario where the former represents the latest outflow event and the latter is a fossil outflow.

Zoom In Zoom Out Reset image size

To estimate the dynamical age of the outflows, T_dyn, we take the ratio of the projected length (measured from the protostar), R_max, to the the maximum speed of the outflow, V_max, corrected for the inclination i: ${T}_{\mathrm{dyn}}({V}_{\max })=({R}_{\max }/{V}_{\max })/\tan (i)$. ¹⁴

Toward W51e2e, the outflow is clearly bipolar and fairly symmetric, though the redshifted side is partly obscured by the dust continuum; therefore, we use the blue lobe in our estimates. We observe a maximum velocity of 105 km s⁻¹ (v_LSR = −50 km s⁻¹) at a separation of about 2560 au, which gives a dynamical age ${T}_{\mathrm{dyn}}({V}_{\max })=116\times (\tan (45^\circ )/\tan (i))$ yr, where the inclination i is expressed in units of degrees.

In W51e8, the highest observed outflow velocity is 42 km s⁻¹ at a separation of about 390 au, which corresponds to ${T}_{\mathrm{dyn}}({V}_{\max })\,=44/\tan (i/45^\circ )$ yr. However, we warn the reader that the outflow appears to be nearly in the plane of the sky; therefore, the age estimate is not reliable.

In W51north, we first consider the blueshifted lobe, which displays a much simpler structure. We observe a maximum velocity of 95 km s⁻¹ (v_LSR = −35 km s⁻¹) at a separation of about 1540 au, which gives a dynamical age ${T}_{\mathrm{dyn}}({V}_{\max })=77\,\times (\tan (45^\circ )/\tan (i))$ yr.

The redshifted lobe from W51north breaks into three major components, representing three consecutive events; these are labeled as "young," "transition," and "old" in Figure 7. The "young" flow is the component along the NW–SE close to the protostar. For a maximum velocity of 66 km s⁻¹ at a separation of about 1300 au, the dynamical age is $\sim 100\times (\tan (45^\circ )/\tan (i))$ yr, higher than inferred from the blueshifted lobe. Adopting an average velocity of ∼30 km s⁻¹ at a separation of 800 au (where the second component starts), one gets 125 yr. The second component, which identifies the "transition" between the young NW–SE component and the old north–south component, extends from 800 to 1900 au, corresponding to 300–125 yr (or a duration of 175 yr) using an average velocity of 30 km s⁻¹. The third component is the larger-scale "old" north–south outflow and extends up to 08 or about 4300 au, yielding an age of 680 yr for a 30 km s⁻¹ velocity. An even older component of the outflow is traced by ¹²CO, up to a length of about 4'' or 22000 au, corresponding to an age of about 3500 yr (for an average velocity of 30 km s⁻¹).

In order to determine the momentum of the outflowing gas, we compute its total mass from the SiO emission integrated in the full velocity range. We follow a standard approach in two steps.

As a first step, we measure the number of SiO molecules per unit area along the line of sight. We first calculate the column density of SiO in the energy level corresponding to the observed transition J = 5–4 using the measured intensity of that transition (channel by channel) and the Boltzmann equation for statistical equilibrium coupled with the standard radiative transfer equation, e.g., the optically thin version of Equation (30) in Mangum & Shirley (2015). Then we relate the number of molecules in the given energy level to the total population of all energy levels in the molecule, assuming that they are populated according to a Boltzmann distribution; we use the rotational partition function, a quantity that represents a statistical sum over all rotational energy levels in the molecule, and we assume a constant temperature defined by the excitation temperature T_ex, e.g., Equation (31) in Mangum & Shirley (2015). We adopt T_ex =  250 K (note that the partition function is approximately linear with T_ex in this regime, so if the molecules are twice as hot, the column estimate should be halved).

In the second step, we use the SiO column density to compute the total mass of the molecular gas. To convert from the SiO column to the H₂ column, we assume a conservatively high abundance of SiO: N(SiO)/N(H₂) = 10⁻⁷. The total momentum is then obtained by integrating over the full velocity range of SiO emission: P_o  = ∑_v m_v  v.

In our analysis, we treated the blueshifted and redshifted lobes independently, but the latter provided mass values 4–10 times lower. In fact, in W51e2e, the redshifted lobe is partly obscured by the dust continuum, whereas in W51north, the redshifted side shows a more complex structure than the blueshifted side, yielding more uncertain estimates. To mitigate these observational biases, we used the blue flows alone in our estimates. The outflow masses and momenta derived with this analysis are reported in columns 6 and 7 of Table 2. In order to account for both outflow lobes, we thus need to multiply those numbers by a factor of 2, yielding total masses of 0.36, 0.72, and 0.22 M_⊙ and total momenta of 18, 50, and 3 M_⊙ km s⁻¹ for the outflows driven by W51north, W51e2e, and W51e8, respectively.