SOFIA and ALMA Investigate Magnetic Fields and Gas Structures in Massive Star Formation: The Case of the Masquerading Monster in BYF 73

We mapped BYF 73 on 2019 July 17 at 0832–0905 UT with HAWC+'s band D (λ154 μm) filter. ¹¹  Chopping and nodding were done asymmetrically due to the nearby FIR emission to the Galactic west and south. The total on-source integration time was 784.4 s. Pipeline processing with HAWC-DRP produced final Level 4 quality image products that were downloaded from the SOFIA archive. This processing produces data that has all known instrumental and atmospheric effects removed, giving an absolute Stokes I calibration uncertainty of 20%, a relative polarization uncertainty of 0.3% in flux and 3° in angle, and astrometry that should be accurate to better than 3'' (Harper et al. 2018). However, we found the HAWC+ L4 astrometry was still consistently offset ∼2'' to the Galactic south compared to the Gemini 10 μm, Herschel 70 μm, and ALMA and ATCA 3 mm maps, all of which are strongly and consistently peaked on the massive protostellar core MIR 2 (allowing for MIR 1's proximity to MIR 2 in the Gemini data), so we inserted this correction by hand into the HAWC+ data files.

At a distance of 2.50 ± 0.27 kpc (near NGC 3324; Barnes et al. 2010; Samson 2021), the scale for BYF 73 is 001 = 36'' = 0.44 pc, or 0.1 pc = 825 = 00023. Thus, HAWC+ band D gives us a useful SDR from 0.16–3.6 pc, a linear factor of 22 and almost 500 resolution elements in area. The resulting full-field images in both total intensity Stokes I and the debiased polarized flux $P^{\prime}$ = $\sqrt{{Q}^{2}+{U}^{2}-{n}_{\mathrm{rms}}^{2}}$, where n_rms is the combined instrumental and sky noise, are shown in Figure 1, overlaid also with the inferred B field polarization vectors.

Zoom In Zoom Out Reset image size

BYF 73 was observed with ALMA at 3 mm wavelength on 2020 January 1 in the C-43 array (baselines 15–314 m) and in two correlator setups and mapping modes; the total on-source integration time was ∼7500 s. The first mode mapped a standard, 13-pointing mosaic of size 28 centered on the peak molecular line emission as measured in the Mopra maps (Barnes et al. 2010), similar in extent to the ATCA mosaic reported by P18, in both the 3 mm cold dust continuum plus the J = 1 → 0 lines of ¹³CO and C¹⁸O and N = 1 → 0 line of CN. The second mode was a single-pointing, full-polarization, deeper integration at the peak emission position near MIR 2, to map the B field strength and structure with (1) the cold dust continuum, and potentially with (2) the Goldreich–Kylafis effect in the line wings of ¹²CO (Goldreich & Kylafis 1981), and/or (3) the Zeeman effect in the hyperfine structure of CN (e.g., Hakobian & Crutcher 2011). Mosaicking with ALMA was not available for polarization modes in Cycle 7.

Standard reduction pipelines were applied to the data, including bandpass, complex gain, flux, and polarization calibration on nearby quasars; images were formed by a joint deconvolution for the mosaics with cleaning and restoration as implemented in the CASA task tclean; and primary beam correction was made within a cutoff at 0.2. The resulting science-ready FITS files were either downloaded from the ALMA Science Archive for the pipeline-reduced data, or the ALMA-North America servers for manually reduced data not included in the automated pipeline processing.

The continuum mosaic is shown in Figure 2 (top panel), with a maximum recoverable scale (MRS) ∼55''. This is larger than the nominal single-field ALMA MRS due to the joint deconvolution for a mosaic, which recovers some of the larger spatial scales missed in a single pointing (see caption and Section 3.1). The mosaic also produced data cubes of the ¹³CO, C¹⁸O, and CN emission at 0.16 km s⁻¹ velocity resolution, but with slightly smaller beams and MRS compared to the continuum, due to the higher frequencies; see Section 4 for results and details.

Zoom In Zoom Out Reset image size

In the top panel of Figure 2 it is clear that, except for the extensions off-frame to the northwest and southwest (i.e., into the H ii region), the HAWC+ and ALMA continuum I maps seem to generally trace the same structures, including the weaker point sources to the east and along the ionization front to the north, despite the 20- and five-fold difference, respectively, in wavelength and resolution.

In the single polarization field, the MRS = 28'' for all polarization products, and the synthesized beams are also the same. An image of the debiased polarized continuum flux $P^{\prime}$ is shown in Figure 2 (bottom panel), but vignetted to exclude spurious features due to missing short spacings along the north and south field boundaries. The $P^{\prime}$ image is also overlaid with the inferred B field polarization vectors.

At 2.5 kpc, the ALMA mosaics of BYF 73 give a useful SDR from 0.030–0.67 pc (6,000–140,000 au): coincidentally with the HAWC+ maps, this is a linear factor of 22 as well.

It is instructive to compare the earlier 90 GHz ATCA data (P18), which only detected MIR 2 within the mapped mosaic of BYF 73, with the 50× more sensitive ALMA 102 GHz images. The inferred flux density of MIR 2 was 50% higher (34 mJy) in ATCA's ∼2× larger synthesized beam than in Figure 2, but on convolving the ALMA data to the ATCA resolution, we recover the identical flux density for MIR 2. Further, P18 found that MIR 2's 90 GHz continuum and 89 GHz HCO⁺ line flux were <30% of the values expected from the ∼arcminute-resolution spectral energy distribution (SED) fit to BYF 73's Herschel data (∼120 mJy; Pitts et al. 2019) and the Mopra single-dish line flux (Barnes et al. 2010), respectively, which P18 attributed to a smooth overall structure in BYF 73 that ATCA apparently resolved out. This turns out to be very close to half-true: we find the total flux density in our mosaic to be ∼70% of the projected SED value at 102 GHz, despite similar shortest baselines in both interferometers.

These results are explained by the millijansky-level structures around MIR 2, which were too weak to be separately detected in the ATCA map (noise σ_rms = 7 mJy beam⁻¹) but raised the measured flux density in the unresolved structure of MIR 2; also, these structures together contribute ∼half the additional flux expected from the SED fit to the ALMA map. With this insight, we see that the older ATCA and current ALMA data are completely consistent with each other, allowing for their respective sensitivities.

We also note that the deconvolved size of MIR 2 in the ALMA data is measured to be 32 × 28 in both the mosaic and deeper polarization field, which is only slightly smaller than the 42 × 30 derived from the ATCA map despite its ∼2× larger synthesized beam (∼4× in area). Therefore it seems MIR 2's protostellar structure is close to being resolved at this scale, 3'' = 7500 au, and future subarcsecond imaging may reveal useful information about its mass distribution.

The spatially resolved SED fit to Herschel data of Pitts et al. (2019) not only provides the missing short-spacing flux information as above, but also allows for the calculation of a merged single-dish and ALMA 3 mm continuum image. The SED fits were used to project how BYF 73 would look at 3 mm with Herschel's λ500 μm resolution of 36'', assuming that at 3 mm, MIR 2 has a negligible contribution from free–free emission in an unresolved UCH ii region, since that would tend to push MIR 2's flux density above the SED fit. The derived image was combined with the ALMA map via the Miriad task immerge (Sault et al. 1995) to recover the missing flux density in Figure 2 that resides in larger angular scales. The result changes the appearance of the image very little; it also does not change the brightness of the individual small-scale structures, except to fill in the broad (∼50'') but shallow (∼0.3 mJy beam⁻¹) negative bowl underlying MIR 2 and its environs. This shows that the missing 30% of the flux density is distributed very smoothly across BYF 73 after all.

To provide context, we show in Figure 3 composite mid-IR images (similar to that in P18), overlaid with selected contours of the HAWC+ and ALMA data from Figures 1 and 2. We have added five more mid-IR point-source designations to the eight identified by P18. MIR 9 and 10 are mid-IR bright in the IRAC images, especially band 2 (4.5 μm), and would have been easily detected with T-ReCS (8–18 μm) if the imaged area had been slightly larger. MIR 11–13 are really millimeter-continuum sources, but we use the mid-IR designations to avoid new nomenclature that might be confusing.

Zoom In Zoom Out Reset image size

At 7 mJy, MIR 11 is the second-brightest point source in the ALMA mosaic, and is also clearly detected by HAWC+ (154 μm). However, in IRAC bands 2–4 (4.5–8.0 μm), it is detected not as a point source, but as a biconical nebula around the 3 mm position (Figures 3(e) and (f)), with a bright northeast lobe and fainter southwest lobe; evidently the central object is deeply embedded and undetected shortward of 10 μm. The northeast lobe also shows redder mid-IR colors at locations closer to MIR 11 itself. These are all classic hallmarks of another (massive) protostar. MIR 13 is weakly detected at K and IRAC bands 2–4, while MIR 12 is not detected shortward of 10 μm, like MIR 11. Together, however, MIR 12 + 13 are marginally detected in the HAWC+ I image as distinct extensions to the FIR emission; in combination with their 3 mm continuum detections at a signal-to-noise ratio (S/N) ∼10, we consider them to also be probable (lower-mass) protostars. In the spectral-line data (Section 4), while MIR 11–13 are all outside the single ¹²CO field of view, the mosaics show interesting features near MIR 11 consistent with its protostar status (Section 4.5). In contrast, the mosaics near MIR 12 and 13 are completely unremarkable. Based on the spectroscopy, none of MIR 11–13 seem to have any impact on the wider cloud's evolution.

This extensive multiwavelength data, showing a relative paucity of millimeter-wave point sources and almost-as-scarce mid-IR (i.e., 8–18 μm) point sources, supports P18's inference that most of the plethora of near-IR (i.e., 1–5 μm) stars are likely to be in the foreground of the BYF 73 cloud. That is, while scores of stars within the T-ReCS field show embedded near-IR colors (Andersen et al. 2017), most of these cannot be deeply embedded, since P18 only directly detected eight of them with T-ReCS, suggesting a lack of embedding envelopes. Based on comparisons with their near-IR visibility, T-ReCS would likely only have detected two more sources outside the observed mosaic, MIR 9 and 10, at P18's sensitivity level.

Even among these 10 mid-IR bright sources, only MIR 2 is detectable at all in the 3 mm continuum; specifically, even the very mid-IR-bright stars MIR 1, 3, 9, and 10 are not detected with ALMA. By comparison with MIR 2, this suggests that these other four mid-IR bright stars have very minimal (if any) protostellar dust envelopes, of mass <3–4 M_⊙ (ALMA's 3σ detection limits in the two observing modes). Therefore, it is reasonable to conclude that, among all of these point sources, only MIR 11–13 have similarly "cold" SEDs to MIR 2, and are in a similarly early stage of protostellar evolution. Scaling MIR 11's 3 mm flux density to MIR 2, which is three times as bright, suggests that its mass may very approximately be 80 M_⊙, still large by protostellar standards. Similarly scaling MIR 12 + 13's 3 mm flux densities yields dust masses ∼7 M_⊙ and 10 M_⊙, respectively.

For the extended emission, both the polarized and unpolarized 154 μm structures simultaneously trace two very different dust populations, each in their own way: (1) the warm dust permeating and surrounding the H ii region, arcing out to the west and northwest from the molecular clump, and seen well in Herschel and Spitzer images at 70 μm and shorter wavelengths, and (2) the cold dust in the massive (2 × 10⁴ M_⊙) molecular clump to the east, traced well by the usual millimeter-wave molecular lines and longer wavelength (≥250 μm) continuum.

Similarly, the 3 mm emission mostly traces the cold molecular structures, but apparently also some high-density warm dust associated with the north–south ionization front (IF) between the molecular cloud and H ii region. Circumstantially, MIR 3 appears to be the main driver of the IF in Figures 3(a)–(d); MIR 4 also lies close to the southern end of the IF, but seems not to have as much impact on its surroundings. The main extended features in the 3 mm continuum are the rather striking arcs of emission running mainly east and west of MIR 2, which, for lack of a better term, we call the "Streamer." ¹²  There is also a notable ∼5'' × 10'' gap (the "Hole") in the 3 mm emission within the Streamer, immediately adjacent to MIR 2 on its eastern side. It is unclear from its continuum properties whether this is a true lack of emission due to an absence of material in the Streamer, or whether it is the shadow of an extremely cold, high-optical-depth component in the foreground of the Streamer, completely absorbing the 3 mm emission beyond it. The spectral-line maps, however, resolve this question; see Sections 4.1 and 4.4.

We begin our exploration of the molecular core's B field as revealed by HAWC+. Zooming in to the inner portion of the molecular cloud near the massive protostar MIR 2 (P18), a very striking feature of the polarized emission stands out immediately—see Figure 4. There is a strong, narrow null in $P^{\prime}$ curving around the western side of the molecular peak, between it and the H ii region, in particular the darker blue colors indicating very low $P^{\prime}$. The east–west width of this null is quite small, apparently only 1 pixel or ≲7000 au (i.e., much less than the angular resolution of HAWC+, but we show that this is not unphysical below). Further, this null can be traced most of the way around the molecular peak, although it broadens out somewhat to the north, south, and east, to an approximate width of 3–6 pixels, or 0.1–0.2 pc.

Zoom In Zoom Out Reset image size

The area enclosed by this boundary layer, signified by where $P^{\prime}$ rises above 0 on the inside, is an approximately elliptical zone of size 0.55 × 0.40 pc at PA ≈ 80° (labeled IBL for inner boundary layer in Figure 4). The outer boundary layer of the null is a vaguely elliptical polygon of approximate height 0.91 pc and width 0.60 pc, at PA ∼20° (labeled OBL in Figure 4). The space between these boundary layers has essentially zero (i.e., S/N < 3) FIR polarization and B_⊥.

To see why the narrow western null is real, we show in Figure 5 the Stokes Q and U maps in the same zoomed area as Figure 4. Carefully comparing the three images in Figures 4–5 on the western side of the MIR 2 core, one sees that where $P^{\prime} \approx 0$ (e.g., close to the I = 5.6 Jy pixel⁻¹ contour), Q drops from positive to negative values going into the center, and U behaves similarly (Q and U are both 0 where the images are orange). Therefore, the $P^{\prime}$ null is very narrow here because both Q and U are changing rapidly through 0 at the same locations (but in a manner consistent with HAWC+'s angular resolution), from larger positive to larger negative values as one traverses into the center of the core.

Zoom In Zoom Out Reset image size

In other words, each of these components of $P^{\prime}$ is strongly reversing sign on this boundary, corresponding to a null in $P^{\prime}$ and 90° change in direction (due to the definition of the Stokes parameters) for both the observed polarization angle and inferred B field direction ${\theta }_{{B}_{\perp }}$ across this boundary. West of MIR 2, this change is very sharp, much less than a beamwidth. Indeed, this change of direction is visible in the $p^{\prime}$ vectors themselves, as shown in Figure 4. The vectors just outside the null, i.e., in the H ii region and also east of the MIR 2 core, are all oriented roughly east–west, while the vectors inside the null, especially close to the MIR 2 core, are mostly oriented roughly north–south. The change in PA is sharpest where Q+U are reversing most sharply, just west of MIR 2. For the small-$P^{\prime}$ locations between the two boundary layers to the north, east, and south of the MIR 2 core, the changes in sign for Q and U are more gradual. But they do change sign in each case, when looking from the outer areas of the cloud to the center.

In the more gradually changing boundary north, east, and south of the core, the $P^{\prime}$ nulls seem to indicate an area where the 2D plane-of-sky B field component (B_⊥) is dropping to zero, since the $P^{\prime}$-null boundary's width is roughly beam-sized and there are few pixels with $P^{\prime}$>3σ. In contrast, the sharp $P^{\prime}$-null boundary west of MIR 2 is much thinner. Rather than merely dropping to zero, this seems to be due to a 90° change in direction alone, i.e., that B_⊥ is sharply changing in direction at the boundary layer around MIR 2, producing a null in $P^{\prime}$ without a corresponding null in B_⊥.

An alternative situation at the IBL is that the 3D B vector is aligned very close to our line of sight there, i.e., that B consists only of B_∣∣. In other words, looking progressively from the H ii region (west) through the IBL and toward MIR 2 and its immediate east, the 3D B field orientation points one way (∼EW) in the H ii region, twists through 90° at the null so that B points directly toward or away from us at the IBL, and then twists another 90° to point in another 3D direction (∼NS) nearer to MIR 2 within the IBL, with Q and U changing sign in the process. We consider this a less likely proposition, however, since the $P^{\prime}$ is so narrow, any B_∣∣ would put additional unresolved structure into the IBL, and it would require a rather large flux annulus to be pointed right at us, a fairly significant "finger of God" effect in our view. On the other hand, a pure 90° change to B_⊥ alone might require a discontinuity in B_⊥. To resolve this, higher-resolution polarization data could reveal more details to this narrow change in B_⊥ (see Section 3.4), while Zeeman measurements could measure B_∣∣ directly, potentially distinguishing between a single 90° B_⊥ twist or additional B_∣∣ in the IBL.

A third possibility for the IBL's apparent null is that the polarization signal comes from dust emission on the low-opacity (west, H ii region) side, as opposed to dust absorption on the high-opacity (east, molecular core) side, making the null more of a polarization radiative transfer effect, without necessarily implying anything significant for the cloud's inherent B field. This would be similar to the FIR polarization pattern in Sgr B2, although observed at a much coarser physical resolution of 1.5 pc (35'' beam) with the Kuiper Airborne Observatory (Novak et al. 1997). In order to be relevant, the dust opacity in the core would need to be ≳1; however, based on P18's SED fitting to Herschel and other data, we compute an effective peak τ ≈0.001 for BYF 73 at 154 μm and 37'' resolution. Allowing for HAWC+'s 7.4× smaller beam area, it is unlikely that the peak τ within the IBL is more than 7× higher than this. Even in ALMA's beam, ∼30× smaller in area than HAWC+'s, τ could rise to ∼0.2 at 154 μm if all of the flux were coming from MIR 2, but this is still less than 1, and we know MIR 2 contains less than half the flux there, Section 3.1. Therefore, we can discount this possibility. For now, we prefer the pure B_⊥-twist interpretation since it is the simplest.

Upon further inspection of the inner 0.5 pc ∼001, one can see another significant feature around MIR 2 in the polarization maps, even where the central I structure is very smooth. While the I and $P^{\prime}$ emission peak exactly on MIR 2's position, the I morphology is slightly more extended to the east, compared to the sharper decline toward the west/H ii region. This morphology is mimicked in the inner $P^{\prime}$ distribution, i.e., where $P^{\prime}$ > 0 inside the IBL, except that in $P^{\prime}$ the point-source nature of MIR 2 is much more distinct, while the extension to the east is revealed as a semicircular ring structure adjacent to MIR 2. We dub this the "eastern polarization lobe" (hereafter EPL). Morphologically, it seems unlikely that this lobe is associated with any of MIR 1–8, as can be seen in Figure 6, which overlays their positions on both the $P^{\prime}$ and $p^{\prime}$ images. This is underscored by indications from P18 that MIR 1 and 3–8 are possibly on the near side of the BYF 73 cloud, and not as deeply associated with the MIR 2 core.

Zoom In Zoom Out Reset image size

Intriguingly, the eastern polarization lobe (EPL) shows a similar polarization signature to the MIR 2 core proper (see top panel of Figure 5) but inverted in Q, going from negative values outside the EPL to positive values across it. This is equivalent to a sharp rotation of the inferred B field between each structure as we will see in the next section, further suggesting that the EPL is distinct from the MIR 2 core/protostar, each having its own physics, despite the much more amorphous appearance around MIR 2 in Stokes I (Figure 1). Identification of these features was based solely on the HAWC+ data, and before the ALMA maps (next) were in hand.

Turning now to the ALMA data, the only high-column-density dust structures seen in the 3 mm continuum (Figure 2) are (i) MIR 2; (ii) the 1–3 mJy beam⁻¹ structures east and west of it, which we call the "Streamer" and "Streamer-west," respectively; (iii) a ∼0.5 mJy beam⁻¹ linear feature aligned almost exactly north–south with the H ii region's IF; (iv) another ∼0.5 mJy beam⁻¹ patch northeast of MIR 2 aligned with the EPL; (v) some even weaker diffuse features ∼1' to the east of MIR 2; plus (vi) three other eastern point sources, which we have designated MIR 11–13 in Figure 3. The larger-scale features in Stokes I at the two different wavelengths have a very nice overall correspondence, despite the different resolutions: the east–west extension of the brightest core emission, the weaker emission extending north along the IF, and the new point sources MIR 11–13 and diffuse emission extending to the east all look mutually consistent. Even the EPL's structure, inferred from the HAWC+ data alone, is easily and gratifyingly verified in the ALMA I images. The detectable ALMA polarization, however, is limited to a subset of these features, namely MIR 2, both sides of the Streamer, the EPL, and possibly the southernmost parts of the IF (near MIR 4).

The ALMA $P^{\prime}$ map is therefore much more spatially compact than the HAWC+ $P^{\prime}$ map. However, the ALMA $p^{\prime}$ values in the molecular core are quite large, typically 5%–20% or more in high-S/N areas, as opposed to the more typical HAWC+ $p^{\prime}$ values around 1% within the IBL/molecular core (HAWC+ $p^{\prime}$ is 10% or more only in the H ii region, but does rise to ∼5% in the diffuse, eastern extremes of the molecular cloud). This lower percentage polarization at shorter wavelengths could be due to two effects:

(A) The polarization signal is being diluted in the larger HAWC+ beam due to its origin in small structures, such as those found in the ALMA maps, but which to some extent cancel each other out in the HAWC+ beam. For example in the MIR 2 core, the correspondence between HAWC+ and ALMA vectors is modest, and the ALMA vector PAs vary more strongly than the HAWC+ PAs. However, due to the correspondence between both HAWC+ and ALMA inferred B field morphologies described below, we discount this effect.

(B) More probably, in the denser parts of the cloud, the 3 mm emission is more efficiently polarized by the cold dust than the 154 μm emission: the "polarization spectral index" (PSI) is >1. This would run counter to the situation in the ρ Oph cloud, where radiative torques from external illumination are thought to more efficiently align grains in the less dense parts of that cloud, giving a PSI < 1 (Santos et al. 2019). Here, we argue that MIR 2's radiation could be aligning grains more efficiently in the cloud core, if radiative torques from internal illumination are the cause (Lazarian 2007).

We overlay both instruments' polarization maps of the IBL, the peak column density area in all maps, in Figure 7. Within the cold, high column density dust preferentially traced by the 3 mm maps, we note two distinct magnetic domains comprising five substructures, each with its own orientation and character.

Zoom In Zoom Out Reset image size

(1a) Close in to MIR 2, the field is oriented mostly north–south, which is very similar to that inferred from the HAWC+ data, but with amplitudes $p^{\prime}$∼1% for HAWC+ and $p^{\prime}$ ≳ 3% for ALMA. We call this the "MIR 2 core."

(1b) Just to the southeast of MIR 2, at the western end of the main Streamer, there is a small patch of polarization with a similar north–south orientation, which we call the "MIR 2 extension."

These two structures comprise the predominantly north–south magnetic domain inside the IBL. The following three structures comprise a different magnetic domain, oriented mostly east–west or somewhat northeast–southwest.

(2a) Across the EPL, the uniformity is almost as good as in the MIR 2 core, with most HAWC+ vectors $p^{\prime}$∼1% @ N60°E, while the ALMA vectors range over $p^{\prime}$ = 10%–20% and run mostly east–west, although some vectors turn toward ∼N60°E at the more distant fringe from MIR 2.

(2b) Along the main Streamer east of MIR 2, ALMA vectors are $p^{\prime}$∼5%–15% while running mostly east–west nearer to MIR 2, but again turning more toward ∼N60°E as they move away from MIR 2. HAWC+ does not detect high-S/N polarized emission from the Streamer; thus, its vectors are somewhat jumbled in orientation there, but their alignment with the ALMA vectors is still reasonably good.

(3) In the Streamer-west, while the HAWC+ vectors continue to align north–south, the ALMA vectors turn east–west, but this is beyond the reliably calibrated polarization radius in the ALMA field, so the divergence may not be significant. It is also possible that, because of the larger beam, the HAWC+ data are dominated by the bright emission from MIR 2 (polarized north–south) farther into the Streamer-west, IF, and H ii region than are the ALMA data, before HAWC+ finally picks up the east–west B field orientation in the H ii region itself. If true, this would make the polarization signal from both instruments more consistent with each other here too, as per effect B above.

In terms of the $P^{\prime}$ null and sharp 90° B_⊥ twist seen in the HAWC+ data, Figure 7 seems to suggest that it is mostly an artifact of resolution and sensitivity, as per the pure B_⊥-twist explanation (Section 3.3). In other words, we can see the inferred B_⊥ direction change quickly between the MIR 2 core and the Streamer-west, right under the edge of the north–south B_⊥ vectors in Figure 7, even if the ALMA Streamer-west vectors are less reliable there.

In summary, the significant B field structures in the molecular core of BYF 73 are MIR 2, the EPL, and the Streamer, where both the HAWC+ and ALMA inferred B fields are broadly consistent. The B field structures seen by both facilities in the H ii region may also be consistent with each other. We reserve discussion of the B field structures in the H ii region for Section 5.1.

The continuum structures seen in the molecular core at both the SOFIA/HAWC+ and ALMA wavelengths are intriguing, both in total intensity and polarized emission. However, while apparently related to the dominance of MIR 2, from their structure alone their physical significance is not entirely obvious, nor is how they are connected to BYF 73's star-forming activity.

Not surprisingly, the ALMA spectroscopy provides important insights; perhaps more surprising is that it is the ¹²CO data that provide the key. Although the ¹²CO spectral polarization cubes only cover a single ALMA field as in the bottom panel of Figure 2, the information they reveal about the nature of the continuum emission in BYF 73 is pivotal. First, the brightest ¹²CO emission by far lies in the highly Doppler-shifted line wings, extending up to ±35–40 km s⁻¹ from the cloud's systemic V_LSR, as illustrated in Figure 8. Spatially, these line wings delineate a massive bipolar outflow clearly emanating from MIR 2. The opening angle appears small near MIR 2, θ_open≲10°, and the outflow appears to impact the Streamer: the flow directions are apparently strongly affected by the large inertia of ambient cloud material in the Streamer, and deviate from their initial vectors. Based on the small θ_open and lack of overlap between the red and blue wings, we estimate the outflow's inclination to the line of sight lies in 40° ≲ θ_incl ≲ 80°.

Zoom In Zoom Out Reset image size

From detailed inspection of the ¹²CO Stokes I cube, the intrinsic outflow direction from MIR 2 is along the magenta line in Figure 8 at a Galactic PA = 120°, but this terminates at the magenta boxes at each end of that line. The redshifted outflow then deviates around both sides of the Streamer-west along the paths indicated by the gold arrows in Figure 8. The blueshifted outflow is apparently deflected by the highest-density portion of the Streamer into the direction shown by the cyan arrow of Figure 8. Spectrally, the highest relative speeds appear to lie near MIR 2 in the red wing, but are displaced from MIR 2 by ∼14'' = 0.17 pc downstream in the blue wing. Apart from this offset, the outflow speeds are generally more modest as one looks farther downstream.

In the case of the blue wing, the outflow direction along Galactic PA = 120° close to MIR 2 is seemingly deflected by a clear 47° "bounce" into a single new direction along PA = 73°. The flow then continues to at least the eastern edge of the single polarization field, 0.6 pc away from MIR 2. This deflection is clearly seen in the individual ¹²CO Stokes I cube's channel maps, and is not an artifact, for example, of opacity-masking of a more southerly portion of the blue wing by the Streamer, hiding a continued blue outflow along PA = 120°. This is because (1) the Streamer and outflow are well separated in V_LSR, so there is no opportunity for the Streamer to mask some parts of the blue wing (see also below); and (2) the ¹³CO line, which is much lower opacity than ¹²CO, shows the same structural features coincident with the deflection 9'' east of MIR 2.

In the case of the red wing, the initial outflow from MIR 2 along PA = –60° appears to be somewhat blocked by the Streamer-west, such that the flow deviates to either side of this obstruction, before continuing to flow at the same PA to the western edge of the field, 0.3 pc away from MIR 2. These deviations are similarly easy to see in the channel maps.

Note that the outflow widths, at ∼10''–15'', are well under the ALMA MRS in the single ¹²CO field, so we believe we recover essentially all of the outflow structure in the line wings. It is difficult to say, however, if the outflow continues beyond these boundaries (e.g., into the H ii region), since the ¹²CO data are limited by the field of view. But it is fairly obvious that the outflow and Streamer interact strongly, one sculpting the other, including the appearance of the EPL and Hole.

The second noticeable feature of the ¹²CO data, apart from where the outflow can be specifically traced in to MIR 2, is that the rest of the cloud is much fainter in the line core between roughly −22 and −17 km s⁻¹, with any non-outflow features being <10% as bright. This confirms that the cloud is extremely optically thick everywhere, and that except for the large outflow-driven Doppler shifts, virtually no ¹²CO emission can escape from the cloud's interior.

Third, and even more interestingly, we detect the linearly polarized Goldreich–Kylafis effect almost everywhere in the outflow, and at high S/N in P_d in almost all channel that trace the outflow in I. This is presented and analyzed in Section 5.3.

The ALMA ¹³CO, C¹⁸O, and CN mosaics provide further details for analysis of the BYF 73 cloud emission, but we focus here on kinematic features associated with the Streamer and MIR 2, in order to shed further light on the B field structures described above, and their dynamics.

In contrast to the very compact 3 mm continuum emission (compared to the mosaic size, Figure 2), the ¹³CO, C¹⁸O, and CN mosaics illustrated in Figure 9 all show much more extended structure, although this emission is brightest near the continuum features, and for ¹³CO, across much of the IF visible in the Spitzer images as well (Figure 3). The ¹³CO and C¹⁸O emission fills most of the mosaic, seemingly even extending beyond it, in all directions for ¹³CO, and to the north and east for C¹⁸O. Even the CN is somewhat extended, although less so than the isoCO lines. The ¹³CO+C¹⁸O extents include parts of the H ii region (the area west of the IF), presumably due to residual molecular gas on its near and far sides that has not yet been ionized by the UV field or swept up in the general H ii region expansion. Some of this effect is more easily visible in the LV diagram of Figure 10, where the ¹³CO is brightest at velocities slightly redward and blueward of the C¹⁸O across the cold cloud.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In fact, in the data cubes this is very widespread: the effect can be seen at positions and velocities of nearly all structures, even deep within the molecular cloud. There are many extended, often filamentary features with shallow velocity gradients and a distinct ¹³CO layer lying just westward of, and slightly red- or blueshifted from, each bright C¹⁸O structure. Evidently, the cloud is actually somewhat porous to the UV field emanating from the H ii region, despite the cloud's more opaque appearance in the near-IR. The C¹⁸O structures then seem to delineate the colder, more shielded parts of the cloud's interior, while the cocooning ¹³CO around each feature may define its more excited side, facing the H ii region. Curiously, the brightest CN emission seems to track better with bright ¹³CO in position and velocity rather than with C¹⁸O, although widespread fainter CN does lie across the mosaic and various C¹⁸O features. The variation of line ratios with position and velocity is difficult to portray here, but Figures 9 and 10 give some idea of the complexity.

The most noticeable kinematic feature in the spectral-line cubes is an east–west velocity gradient across the Streamer, one remarkable in several respects. Near the bright continuum emission (and thus the brightest parts of each emission line), the gradient is consistent across all three species, reaching a maximum blueshift of −22.0 ± 0.1 km s⁻¹ about 20'' = 0.25 pc east of MIR 2, and a maximum redshift of −17.0 ± 0.2 km s⁻¹ around 11'' ± 2'' = 0.13 ± 0.02 pc northwest of MIR 2 (see Figure 11). The gradient is also at its sharpest exactly across the middle of MIR 2 itself, ∼2–3 km s⁻¹ across only 1 ALMA beam (∼6000 au = 0.03 pc), or ∼75 km s⁻¹ pc⁻¹. The steepest part of the gradient, defining a symmetry axis, lies on a nearly north–south curve, just like the B field orientation in Figure 7. Moreover, this symmetry axis across the middle of MIR 2 (straddling the width of the Streamer) looks identical in all three lines, underscoring its dynamical importance and strongly suggesting a flat north–south feature within MIR 2 as the origin for the outflow (more on this below). Meanwhile, the full extent of the east–west blue-to-red velocity gradient lies along a curved line across l ∼ 28623–28620, roughly 1.3 pc.

Zoom In Zoom Out Reset image size

These details suggest the possibility that the main and western extensions of the Streamer form part of a rather large structure (perhaps including a disk), in which inward flow to MIR 2 must occur and there generate the outflow. While the case is somewhat circumstantial so far, the evidence becomes much stronger upon closer inspection.

For BYF 73 as a whole, we estimate that its systemic velocity is V_sys = –19.6 ± 0.2 km s⁻¹ on the LSR scale, based on inspection of the C¹⁸O data cube. As seen in Figures 9 and 10, nearly all of the mosaics' line emission lies between −24 and −16 km s⁻¹. However, there is clear evidence in the ¹³CO cube of high-velocity line wings close to the position of MIR 2: up to ΔV_blue = −11.5 km s⁻¹ and Δ V_red = +16 km s⁻¹, for a total V_LSR range of 27.5 km s⁻¹ (i.e., from −31 to −3.5 km s⁻¹). This emission is quite small in spatial extent, with length ∼20'' = 0.25 pc and width 10''–15'' = 0.12–0.18 pc for each lobe (see Figure 11).

This compact configuration is completely different to the massive, more extended bipolar outflow clearly visible in ¹²CO (Section 4.1), which is indicated in the bottom panel of Figure 12 to help distinguish the LV patterns of the two species. In particular, the ¹²CO outflow emission that extends beyond an area ±10'' around MIR 2 must be relatively low opacity τ and high excitation T_ex, since it is much brighter than the superimposed ¹³CO emission, where the latter is even detectable at the same (l,V) coordinates. In contrast, the high-velocity emission coincident with MIR 2 has ¹³CO almost as bright as ¹²CO, suggesting a much higher τ and lower T_ex. The respective excitation conditions are consistent with gas being entrained by a powerful mechanical outflow, and gas responding to the local gravitational potential.

Zoom In Zoom Out Reset image size

Finally, the extended, low-velocity ¹³CO wings are not visible in ¹²CO due to the latter's high τ, although much of the ¹³CO line wing emission (Figure 11) is spatially oriented similarly to the inner ¹²CO outflow, including the 43° bend in the blue wing and the deviations around the Streamer-west in the red wing (Figure 8). Line wings similar to either the ¹²CO or ¹³CO high-velocity patterns are not detectable in either the C¹⁸O or CN cubes.

Thus, while the brightest ¹³CO emission is probably also tracing the outflow, the small extent of the high-ΔV emission is more peculiar. It becomes progressively tinier as the velocity channel being viewed moves farther from the cloud's systemic value, contracting to within a beamwidth of MIR 2 at the highest velocities. This is the opposite of what is typically seen in protostellar jets, where the highest velocities are usually at the most distal parts of the outflow (Lee et al. 2000).

Instead, the LV-moment diagrams in Figure 12 suggest this pattern might arise from a Keplerian disk. Sample rotation curves  V_rot = $\sqrt{{GM}/R}$ are overlaid for three different central masses in Figure 12 as well. The only free parameter in fitting such curves is the central mass: ¹³  the position of MIR 2 is well constrained, as is the velocity extent of the emission. Only the highest mass curve of the three examples fits the high-ΔV¹³CO emission envelope adequately. The lower-mass curves and P18's mass estimates are all much too small, and are strongly ruled out under this interpretation.

Indeed, a 1350 ± 50 M_⊙ curve fits the data remarkably well over a longitude extent of 0008 = 0.35 pc, or a radius of 36,000 au, which is about half the longitude extent (286203–216) of the IBL as measured in the HAWC+ data (Figure 6). Such a disk would also be ∼half the size of the Keplerian disk seen in another massive cloud, K3-50 (Howard et al. 1997), and about half or less of K3-50's disk mass (Barnes et al. 2015), so these numbers are not unheard of in high-mass SF. But this means that the central mass of MIR 2 (presumably comprising a massive protostar and its envelope) is five to six times larger than P18's preferred estimate, and that it dominates the dynamics of the gas over that span.

Alternatively, the kinematics can also be explained by gas in free-fall toward a 950 ± 35 M_⊙ object, since V_ff = $\sqrt{2}{V}_{\mathrm{rot}}$ decreases the central mass required to produce the curves seen in Figure 12 by $\sqrt{2}$. Reference to either scenario hereafter is meant to include both as feasible physical settings near MIR 2.

In our view, such rotation/free-fall curves are so distinctive that there is no other reasonable explanation for the motions of the gas, at least within the IBL (i.e., excluding a much smaller amount of apparently counterrotating gas in the same window). Outside the IBL, the deviations from Keplerian rotation are stronger, and may be due to more typical, modestly turbulent cloud motions and/or internal structures.

Even if we have constructed a plausible picture of BYF 73's internal structure and mass flows, the EPL is a distinctive feature in both the SOFIA and ALMA polarization maps (Figure 7) with an as-yet-undetermined role or import. It lies north of the plane of the Streamer/disk around MIR 2, yet the inferred B field directions through it still point toward/away from MIR 2. It is bright in line emission as well (Figure 9), particularly C¹⁸O, but also ¹³CO and CN in turn around its arc. The velocity fields (Figure 10) reveal little except that it is close to systemic (−20 km s⁻¹) in C¹⁸O at its northern apex, and slightly blueshifted (−21.5 km s⁻¹) in both CN and C¹⁸O at its base near the Streamer, but blending smoothly with the Streamer's rotational pattern. Scanning through the channels in the ¹³CO cube, the changing pattern gives the impression of a splash effect or prominence-like tendril, driven by the outflow off ambient Streamer material downstream of the 47° bend in the blue wing, with a relatively gentle relative blueshift of 1.5 km s⁻¹.

If true, this would be quite interesting: does it signify the expulsion of surface material from the Streamer? Or is it a wisp from the wider cloud, falling in at a slightly higher speed, unconstrained by the Streamer due to its northerly approach? While at a lower surface density than the disk, perhaps 10²⁷ m⁻² or a tenth of the Streamer (Section 5.2), the continuum data show that its mass is not insignificant, perhaps 15 M_⊙ in total. Higher-resolution polarization data would be desirable to determine exactly what the EPL signals.

Around MIR 11 there appears to be a second strong velocity gradient, similar to that straddling MIR 2 (Section 4.1). This gradient is somewhat vague in ¹³CO, clearer in C¹⁸O, but most obvious in CN (see Figure 13), and is oriented with the strongest ΔV along a similar axis (∼N40°E) as the biconical mid-IR nebula (Figures 3(e) and (f)): blueshifted emission on the more IR-visible side to the northeast, and redshifted emission to the southwest. The various line profiles show only a narrow velocity range, however, ±2–5 km s⁻¹, rather than a strong outflow signature, and for ¹³CO and C¹⁸O, the red- and blueshifted emission overlaps somewhat on the sky. Both lines have strong self-absorption around the line center of at –19.8 km s⁻¹. The CN is better separated on the sky into blue and red lobes, with no self-absorption. These biconical characteristics suggest the possibility of MIR 11 being in an even earlier, pre-outflow stage of evolution than MIR 2.

Zoom In Zoom Out Reset image size

5.1.1. Preamble

5.1.2. HAWC+ Data Analysis

In the HAWC+ data, the measured ${\theta }_{{B}_{\perp }}$ in any given region with S/N ≳ 5 will have an uncertainty in orientation dominated by instrumental noise, ${\rm{\Delta }}{\theta }_{{B}_{\perp }}$ ≲ 3°. This occurs at a level slightly less than the $P^{\prime}$ = 50 mJy pixel⁻¹ contour in Figures 1 and 3, and we show the corresponding ${\theta }_{{B}_{\perp }}$ maps in Figure 14. There are three regions of interest (hereafter ROI) satisfying these criteria: two arcs of polarized emission in the H ii region (labeled north and south), and the area within the IBL containing the MIR 2 core (but excluding the EPL due to the paucity of statistics, ∼20 pixels). To begin the DCF analysis, we show in Figure 15 the ${\theta }_{{B}_{\perp }}$ distribution in all pixels within each ROI. None of these is really Gaussian, but we overlay such fits in order to compute effective dispersions in ${\theta }_{{B}_{\perp }}$ for reasons that will become clear shortly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next construct histograms for subsets (comprising square boxes of area A) of each ROI, computing a dispersion s(A) in ${\theta }_{{B}_{\perp }}$ for each subset. The smaller the boxes, the more choice we have of where to fit them inside each ROI. We then compute a mean dispersion <${s}_{{\theta }_{{B}_{\perp }}}(A)$> (± a standard deviation) in the field orientation for all small boxes of a given area A, no matter where they are placed within the ROI. Finally, in Figure 16 we plot all such results as a function of box size A^1/2, ranging from a minimal useful size of A = 3 × 3 pixels (roughly one Nyquist sample given the HAWC+ band D beam) to the full size of each ROI.

Zoom In Zoom Out Reset image size

In each ROI, the mean dispersion <s> within all boxes of area A rises as the box size increases, meaning that ${\theta }_{{B}_{\perp }}$ is more correlated on small scales (e.g., s < 5° within the H ii region over spans < 0.5 pc), and becomes less correlated over longer distances (s ∼ 10° across ≳1 pc within the H ii region). The scale at which s seems to plateau is where we identify the correlation length as per Myers & Goodman (1991). Thus, the HAWC+ polarization data suggest that, as far as the B field is concerned, the H ii region consists of one or perhaps two coherent structures, with a correlation length ∼0.5 pc as evidenced by the plateauing of s at 6°–7° in the smaller H ii North (Figure 16), and the slow rise of s across H ii South as the B field orientation slowly changes across the 1 pc arc of the H ii region.

In contrast, the interior of the IBL (already much smaller than the H ii region) shows two closely spaced and distinct B field configurations, namely the EPL and MIR 2 core. While no useful statistics could be compiled for the EPL, even the MIR 2 core has insufficient resolution to discern more than a strongly rising s at all measured scales, and no correlation length can be defined beyond the ∼0.2 pc extent of the core itself, as in Figure 7.

5.1.3. ALMA Data Analysis

The same approach can be used with the ALMA polarization data, which has high enough resolution to resolve the DCF analysis for the MIR 2 core and its environment, unlike its minimal representation in Figures 14–16. We therefore define the ALMA ROIs in Figure 17 where ${\theta }_{{B}_{\perp }}$ has S/N > 2.5 but is typically 5–10, as in Figures 2 and 7 and conforming to the description in Section 3.4. Then in Figure 18 we show the ALMA ${\theta }_{{B}_{\perp }}$ distributions, compiling the ALMA DCF statistics in Figure 19.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We first notice that, in the overlapping range of scales (0.1–0.3 pc), values of the dispersion s in the IBL from both ALMA and HAWC+ data are in agreement, rising approximately from 10° to 20° when averaging broadly over the five structures in Figure 19. Focusing next on the smallest ALMA scales (0.02 pc), s in the IBL also starts out at a few degrees, then gradually rises. In each of the five ROIs, Figure 19 hints at an s plateau for each structure, before rising further as other uncorrelated structures are included. In Table 1 we compile the (${A}_{\mathrm{corr}}^{1/2}$, s_corr) pairs that can be read off the trends in Figure 19, and for completeness also include the values for the H ii region ROIs from Figure 16. The EPL seems to have the fastest-rising s and the least well-defined plateau in Figure 19, suggesting that it has not mapped out a single correlation length in its structure.

Table 1. Davis-Chandrasekhar-Fermi Correlation Statistics for BYF 73 Polarization Structures from SOFIA and ALMA Data

Structure	Correlation	B Field Angular	σB /${\bar{B}}_{\perp }$
	Scale ${A}_{\mathrm{corr}}^{1/2}$	Dispersion scorr
H ii North	0.5 pc	6°	0.1
H ii South	1.5 pc	12°	0.3
Streamer-west	0.14 pc	18°	0.3
MIR 2 core	0.08 pc	16°	0.4
EPL	0.12 pc	22°	0.6
Streamer	0.14 pc	13°	0.3
MIR 2 extn.	0.05 pc	6°	0.1

Download table as:  ASCIITypeset image

Therefore, the measurable correlation lengths in the IBL are perhaps only a tenth of those in the H ii region, 0.05–0.15 pc compared to ∼1 pc. The dispersions in the B field direction at those lengths are, respectively, ∼6°–20°, compared to 7°–12° for the H ii region. Such lengths and dispersions are the scales above which B field directions are not correlated with each other between neighboring areas, and are therefore the scales we should explore for other physical thresholds, and particularly for any constraints on B field strength.

5.1.4. Numerical Results

As described by Ostriker et al. (2001), Crutcher et al. (2004), Barnes et al. (2015), and many other workers, the standard DCF analysis (Davis 1951; Chandrasekhar & Fermi 1953) directly links the dispersions in polarization angle δ θ = s (tracing variations in the B field orientation) to three other physical parameters that simply and naturally describe the propagation of a transverse MHD wave in a turbulent plasma: the gas density n, the line-of-sight velocity dispersion δV, and the plane-of-sky magnetic field strength B_⊥. That is, s will increase as (1) B decreases, since then the magnetic restoring forces are reduced; (2) n increases, since then the medium's inertia to the MHD wave is greater; or (3) δV increases, since that describes the strength of the MHD wave.

According to Crutcher et al. (2004), with appropriate SI unit conversions (1 nT = 10 μG), the projected B field strength

$$\begin{eqnarray}&&{B}_{\perp ,\mathrm{DCF}}=0.543\,\mathrm{pT}\,\sqrt{\mu n}\,({\rm{\Delta }}V/s)\,,\end{eqnarray} \tag{ 1 }$$

where n (m⁻³) is the gas density with mean molecular weight μ = 2.35, ΔV = $\sqrt{8\mathrm{ln}2}\,\delta V$ is the velocity FWHM (kilometers per second) in the cloud, and s is measured in degrees. Included in the constant is a numerical factor Q = 0.5 from Crutcher et al. (2004) to correct for various smoothing effects (e.g., see Ostriker et al. 2001).

For an illustrative example, consider the region of densest gas inside the IBL. From modeling of HCO⁺ emission by Barnes et al. (2010) at the 40'' Mopra resolution (roughly 3× the HAWC+ beam), we take δ V ∼ 1 km s⁻¹ as a median intrinsic value and the estimated peak n ∼ 5 × 10¹¹ m⁻³ ostensibly near MIR 2, to connect s_corr in our polarization maps to the field strength B_⊥. Equation (1) then becomes

$$\begin{eqnarray}&&{B}_{\perp ,\mathrm{DCF}}(\mathrm{MIR}\,2,\mathrm{Mopra})\sim 92\,\mathrm{nT}\,{\left(s/15^\circ \right)}^{-1}\end{eqnarray} \tag{ 2 }$$

at these values of n and δ V around MIR 2 (or 0.9 mG in cgs units). Indeed, the value for n may well be higher in the smaller HAWC+ beam, and is certainly much higher in the 0.03 pc structures revealed by ALMA (see Section 5.2), peaking at 3.6 × 10¹³ m⁻³. Then,

$$\begin{eqnarray}&&{B}_{\perp ,\mathrm{DCF}}(\mathrm{MIR}\,2,\mathrm{ALMA})\sim 1.18\,\mu T\,{\left(s/10^\circ \right)}^{-1}\end{eqnarray} \tag{ 3 }$$

(12 mG), a value which has not been observed in any star-forming region outside of maser spots. Even the smaller value in Equation (2) is among the highest nonmaser B field strengths in similar massive star-forming clouds, according to the compilation of Crutcher (2012, his Figure 7).

Despite the possibly record-setting value for B_⊥ near MIR 2, it is commensurate with MIR 2's high gas density, i.e., together they indicate a mass-to-flux ratio that is very close to critical (see below). In other words, any field strength much less than this (or density much greater) would probably not provide sufficient support against gravitational collapse (allowing, of course, for the likelihood of ∣B∣ > B_⊥). However, this density is derived from SED fitting, which, as we have already noted, may be significantly underestimating the gravitational potential near MIR 2, based on the apparently Keplerian or infalling motions seen in the ¹³CO data (Section 4.3). In that case, even this high B field cannot avoid criticality.

As a contrasting example, we also consider the H ii region ROIs. In such regions, bulk expansion speeds are typically two to three times the velocity dispersion (=sound speed) in the roughly 8000 K ionized gas, ∼12 km s⁻¹ (Habing & Israel 1979; Franco et al. 1990). Such flows are thought to dominate the energetics in the gas. For BYF 73, the H ii region has a total flux density at 843 MHz of only 60 mJy (Green et al.1999). This corresponds to a small emission measure, EM = n² D = 9.5 × 10¹⁵ pc m⁻⁶ (Mezger & Henderson 1967). With an apparent diameter 2R = D ∼ 0.5 pc, this yields a much lower electron density n_e ≈1.4 × 10⁸ m⁻³ than in the molecular cloud, ¹⁴  but we also have a somewhat larger dust-based estimate of n_d ≈ 7 × 10⁹ m⁻³ from SED fitting (Section 5.2), which may average in material from outside the H ii region. We bracket this uncertainty by combining these values in Equation (1) with two estimates (e.g., using a lower μ = 1.28 in the ionized gas),

$$\begin{eqnarray}\begin{array}{rcl}{B}_{\perp ,\mathrm{DCF}}({\rm{H}}\,{\rm\small{II}},\mathrm{ions}) & \sim & 21\,\mathrm{nT}\,{\left(s/5^\circ \right)}^{-1}\,\,\,\mathrm{and}\\ {B}_{\perp ,\mathrm{DCF}}({\rm{H}}\,{\rm\small{II}},\mathrm{dust}) & \sim & 195\,\mathrm{nT}\,{\left(s/5^\circ \right)}^{-1}\,\,,\end{array}\end{eqnarray} \tag{ 4 }$$

for the H ii region ROIs, depending on which parts of the line of sight through the H ii region are being sampled by the HAWC+ polarization data.

Even the lower (pure-H ii) value seems somewhat high compared to a more typical 1 nT in other H ii region studies (Crutcher 2012; Barnes et al. 2015); whether this value is reasonable is unknown, but B field measurements could also be obtained, for example, via high-resolution HI Zeeman observations. The higher dust-based estimate would apply if the polarization contribution is predominantly from outside the H ii region, but then the B field configuration still suggests a connection to the H ii expansion. This could be reconciled with an origin in a sheathing, higher-density PDR layer rather than the ionized cavity.

5.1.5. Energetic Considerations

For our purposes, though, the point is whether the energy density ${\mathfrak{M}}$ in these somewhat strong B fields exceeds or is less than the kinetic energy density ${\mathfrak{K}}$ in the ionized outflow. Borrowing from Equation (15) in Section 5.3, we can write this ratio as $\tfrac{{\mathfrak{M}}}{{\mathfrak{K}}}$ = 5.2% ${\left(\tfrac{{\rm{\Delta }}V/{V}_{\mathrm{rel}}}{s/6^\circ }\right)}^{2}$. From the flow in the H ii region, we have that ΔV/V_rel ∼ 1, while from Figure 16 we have s = 6°, 12° in the H ii region-north and -south, respectively. Thus, the B field energy density in the H ii region is probably still small compared to the kinetic energy in the ionized flow.

This approach is only valid, however, if the situation of the DCF analysis holds, namely the presence of an MHD wave with turbulent motions. If other processes operate, then the B field strength may be indeterminate without direct measurements, either smaller or larger than the DCF value. For example, if other motions enhance variations in ${\theta }_{{B}_{\perp }}$, s may be larger than the DCF-only value, underestimating B_⊥.

In expanding H ii regions, the kinetic energy density of the expansion ${\mathfrak{K}}$ often exceeds the thermal energy density ${\mathfrak{T}}$ by a wide factor (i.e., in addition to exceeding ${\mathfrak{M}}$): ${\mathfrak{K}}$/${\mathfrak{T}}$ = $\tfrac{2}{3}$ M ² (Barnes et al. 2015), where M (=2–3 in the example above) is the Mach number of the flow. For SF in cold molecular gas, the most interesting question is the relationship between B fields and gravity. If ${\mathfrak{M}}$ ≪ ${\mathfrak{G}}$, gravity dominates and the gas is considered "supercritical;" if ${\mathfrak{M}}$ ≫ ${\mathfrak{G}}$, it is "subcritical." We discuss this question further in Sections 5.2, 5.3, and 6.3.

However, there is an additional factor in criticality: we need to understand not only the value of the dispersion s, but also its behavior on different length scales. As explained by Myers & Goodman (1991), where the DCF method applies, these dispersions are related to the ratio of the disordered versus ordered B field strengths via

$$\begin{eqnarray}&&{s}_{\mathrm{corr}}=\displaystyle \frac{{\sigma }_{B}}{{L}^{1/2}{\bar{B}}_{\perp }}\,.\end{eqnarray} \tag{ 5 }$$

Here s_corr is measured in radians, L is the number of B field correlation lengths (assumed ≈A_corr) in the line of sight, ${\bar{B}}_{\perp }$ is the strength of the ordered component of the projected B field, and σ_B is the dispersion in the strength of the random component of the projected B field. For now, we estimate L from the behavior of s as seen in the DCF structure functions (Figures 16, 19), i.e., where s plateaus in each structure at some size scale A, and so constrain somewhat the ratio σ_B /${\bar{B}}_{\perp }$.

In the H ii region, L = 1–2, s_corr ≈0.1, so we estimate (σ_B /${\bar{B}}_{\perp }$)_HII ≈ 0.1–0.3. This is another way of describing the highly ordered appearance of the B field vectors over large scales (Figure 1). Likewise, for the five structures in the IBL, effectively L = (1–3) × ${A}_{\mathrm{corr}}^{1/2}$ by construction, and we estimate the random:ordered B field strength ratios for all seven structures described here in the same way, and list them in Table 1.

The histogram of relative orientations (HRO) method of analyzing B field orientations in star-forming gas is by now a fairly standard technique (e.g., Soler et al. 2017, and references therein). In the Vela C molecular complex, for example, Soler et al. (2017) used BLASTPol data with a resolution 30 = 0.6 pc at Vela to examine how the B field orientation changes with column density. They found that at lower molecular gas column densities ∼10²⁶ m⁻², the B field direction tends to be mostly parallel to or not show any preferred direction relative to gas structures, whereas the B field is mostly perpendicular to higher column density structures ∼10²⁷ m⁻². This generally confirms a series of results by the lower-resolution (∼10') Planck Collaboration (e.g., Planck Collaboration 2016) over a wider range of molecular gas column densities.

In the various structural components of Vela C, the transition from mostly parallel to mostly perpendicular can be rather sharp at a certain column density for each structure, but this transition density is different in each structure. This is widely attributed to a transition from subcritical gas at lower densities, where the flow is at least guided to some extent by the B field, to near-critical or slightly supercritical gas at higher densities, where gravity is capable of overwhelming the magnetic pressure, allowing stars to form.

5.2.1. HAWC+Data

With our substantially higher-resolution ALMA and SOFIA/HAWC+ data, it should be instructive to conduct a similar HRO analysis from several parsecs down to <0.1 pc scales in the massive SF environment of BYF 73. We show first in Figure 20 the relative alignment of B field vectors with the SED-fit column density N map (Pitts et al. 2019) as a proxy for "structure" in the molecular gas, in all of the HAWC+ data as shown in Figure 1. That is, where the rotated polarization vectors ${\theta }_{{B}_{\perp }}$ are aligned with the tangent to the isoN contours, the relative angle is close to 0°, and the field is considered to be "parallel" to the gas structures. Where ${\theta }_{{B}_{\perp }}$ is perpendicular to the contours and aligned with the column density gradient ∇N, the relative angle is close to 90°, and the field is considered to be "perpendicular" to the gas structures.

Zoom In Zoom Out Reset image size

This approach has the advantage of not imposing any preconceived interpretation of whether the gas structures represent "clumps," "cores," "filaments," or any other potentially subjective term (see, e.g., Planck Collaboration 2016; Soler et al. 2017). The distribution is quantified by computing histograms on each N-bin separately as in Figure 21, including the HRO shape parameter −1 < ξ < 1 computed on each N bin's HRO, as described by Soler et al. (2017). This parameter objectively indicates whether there is a preponderance of parallel (ξ > 0) or perpendicular (ξ < 0) alignments in the data, and can be plotted as a function of N (Figure 22) to reveal any trends via linear regression,

$$\begin{eqnarray}&&\xi \,=\,{C}_{\mathrm{HRO}}\,(\mathrm{log}N-{X}_{\mathrm{HRO}})\,.\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Already in Figure 20 we can see that the distribution of relative alignments has definite patterns in various column density ranges. These observations are reflected numerically in Figures 21 and 22. Thus, in the lower-N ranges, there is an overabundance of parallel alignments (relative PA < 20°) between the inferred B field orientation ${\theta }_{{B}_{\perp }}$ and the isoN contours, and ξ > 0 at high significance, 3σ–13σ each across eight N-bins. In the top two N ranges, however, there is a sudden transition to clearly more perpendicular alignments, PA ∼40°–70° in the penultimate N bin (ξ ≈ 0), and 60°–90° in the top bin (ξ < 0 at 6σ). Indeed, compared to Vela C (Soler et al. 2017), the slope C_HRO is substantially closer to –1 in the HAWC+ data for BYF 73, indicating an even stronger alignment trend with increasing N and a sharper transition (X_HRO intercept) than in the Vela cloud, at N_crit = (3.9${\pm }_{0.8}^{1.0}$) × 10²⁶ m⁻². Interestingly, the steepness of C_HRO as seen in Planck+BLASTPol large-scale maps may be correlated with the inclination angle of the mean B field (e.g., Sullivan et al. 2021). One sees a shallower slope in clouds where the polarization fraction levels indicate that the B field is inclined closer to the line of sight. Thus, the steeper slope in BYF 73 may be related to its B field lying close to the plane of the sky (see Sections 4.1, 5.4, 6).

The distribution of points in Figure 20 can be more intuitively understood in Figure 23, which overlays the HAWC+ $p^{\prime}$ vectors and N contours from the Herschel-based SED fits (Pitts et al. 2019). This map shows that the large number of points with N ≲ 10²⁶ m⁻² and preferentially parallel alignments arises in the H ii region, while the other large concentration of points with N ≈ 2 × 10²⁶ m⁻² arises mostly from the extended IF to the north and the similar arc bounding the H ii region to the southwest of MIR 2. The transition between these two column density levels contains relatively few points due to the sharp density gradient across the IF. For N ≥3 × 10²⁶ m⁻² and progressively closer to MIR 2, the alignments become preferentially more perpendicular.

Zoom In Zoom Out Reset image size

5.2.2. ALMA Data

We can repeat the HRO analysis on the smaller scale of the ALMA field. Figures 24–26 similarly show the overall relative alignment distribution as a function of N, the HROs in separate N-bins, and the ξ versus N plot for all ALMA data. The relative alignment distribution shows similarly striking changes with N as in the HAWC+ data. In the three lowest-N bins, the B field shows no particular preference for parallel or perpendicular alignments in the ALMA maps (ξ ≈ 0 within the uncertainties). In the middle six N bins, though, the distribution changes to show a substantial preference for perpendicular structures (ξ < 0 at an S/N of 2.5σ–5σ in each bin). These nine bins together behave similarly to the BLASTPol results in Vela C, again including the existence of a sharp transition from positive to negative ξ, but now at a four-times-higher N ≈1.6 × 10²⁷ m⁻² than in the HAWC+ data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, in the highest-N bin, the distribution changes back to a very strong (7σ) parallel signature, ξ = 0.60 ± 0.08. This produces a very atypical nonnegative result in Figure 26 for the fitted HRO parameter C_HRO (black labels and dashed line). In Vela C and elsewhere, a negative C_HRO means that ξ changes systematically from weakly positive to definitely negative values as N rises, meaning a transition from parallel or random B field alignments to perpendicular ones, often with a sharp transition across ξ = 0 at a particular N. In this context, the last data point in Figure 26 may be anomalous, but as it turns out, this may not be that significant.

To see this, consider the ${\theta }_{{B}_{\perp }}$ and N maps together (Figure 27), where one can see where each of the 10 N bins are located. The three lowest-N bins with ξ ≈ 0 arise in the weaker emission features of the Streamer to the west and farthest east, and southern parts of the IF. The middle six N bins with ξ < 0 arise in the brighter emission of the EPL, MIR 2-ext, and the main part of the Streamer. The highest-N bin arises exclusively from the brightest parts of the MIR 2 peak, where the structure is actually not well resolved in the 26 ALMA beam. This would not only preclude accurate ${\theta }_{{B}_{\perp }}$ measurements at MIR 2, but also might include Q and U cancellation within the ALMA beam, underestimating $P^{\prime}$. Resolution alone would make any alignment inferences questionable, but in addition we recall that MIR 2 is near the limit of the reliably calibrated window of the ALMA polarization field. Therefore, we cannot accurately quantify the alignment measurements or their uncertainties right at the MIR 2 peak, and conclude that the ξ value in this N bin should be discounted.

Zoom In Zoom Out Reset image size

As an exercise, therefore, we also computed the regression parameters C and X for the nine lower-N bins in Figure 26, and show these as red labels and a dashed line. In this case, C is definitely negative (3σ) and more in line with the Vela C results, while the X_HRO intercept gives N_crit = (9.5${\pm }_{3.4}^{5.3}$)× 10²⁶ m⁻². Based on this alone, it seems desirable to obtain an even higher-resolution polarization map of MIR 2 and its immediate surroundings. Such a map would allow us to explore the massive protostellar core's B field in much finer detail and track the ξ trend to even higher N, not to mention better resolving the core itself (e.g., 250 au at 01).

5.2.3. Combined Data

The ξ–N trends in Figures 22 and 26 overlap nicely in column density, and we present a combined plot in Figure 28. There, the slope C is slightly shallower compared to either the HAWC+ or ALMA-only results, but the overall trend is firmer (the uncertainties in C and X are smaller) due to the wider range of N being sampled. In combination, the data suggest that the transition to perpendicularity in BYF 73's Streamer and MIR 2 core occurs (from the red intercept X in Figure 28) at N_crit = (6.6${\pm }_{0.9}^{1.2}$) × 10²⁶ m⁻², near the geometric mean of the transitions from the individual instruments. This can also be converted into an equivalent critical gas density if we assume a line-of-sight depth to the Streamer approximately equal to its projected width, n ∼ N/D, where D ≈ 0.087 pc. Then, n_crit = (2.0${\pm }_{0.4}^{0.5}$) × 10¹¹ m⁻³.

Zoom In Zoom Out Reset image size

This is actually a rather suggestive threshold: in Figure 27, it is the column density of the second-lowest contour, and includes the MIR 2 core, ∼all of the main and Streamers-west, and much of the IF. It suggests that the Streamer's width may be related, locally at least, to the transition between MHD forces governing the gas dynamics and self-gravity, as seen in Section 4.

It is instructive to compare this result with other HRO studies. For example, Planck Collaboration (2016) used 10' resolution Planck data to study B field orientations in 10 nearby (150–450 pc) Gould Belt clouds, with a finest physical resolution similar to our HAWC+ data and ranging up: ∼0.4–40 pc. At this scale the median gas densities are n ≈5 × 10⁸ m⁻³, substantially less than is typical in the Streamer as estimated above. The threshold column densities in these 10 clouds are also lower than that for BYF 73, by a factor of 10 on average, X ≈ 6 × 10²⁵ m⁻². Similarly in Vela C, a massive but relatively unevolved cloud at 700–900 pc (Soler et al. 2017; Zucker et al. 2021) used Herschel and BLASTpol data at 3' resolution for their HRO analysis (i.e., with a similar physical resolution to Planck Collaboration 2016), and found a typical X ≈ 3 × 10²⁶ m⁻³ or about half BYF 73's value. Finally, in two portions of L1688 in Ophiuchus at 140 pc, Lee et al. (2021) combined HAWC+ and Planck data (giving similar physical resolutions to our ALMA data) to confirm a column density threshold similar to Planck Collaboration's (2016) 10 clouds, and a volume density threshold n ≈10¹⁰ m⁻³.

We can relate this column density threshold to the equivalent B field threshold if gravity and the magnetic pressure were critically balanced. Using the mass-to-flux ratio approach of Crutcher et al. (2004) as adapted by Barnes et al. (2015), we have

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{(M/{\rm{\Phi }})}_{\mathrm{obs}}}{{(M/{\rm{\Phi }})}_{\mathrm{crit}}}=0.064\,\displaystyle \frac{{N}_{{{\rm{H}}}_{2}}/{10}^{24}{{\rm{m}}}^{-2}}{{B}_{\mathrm{TOT}}/\mathrm{nT}}\,\,,\end{eqnarray} \tag{ 7 }$$

or

$$\begin{eqnarray}&&{B}_{\mathrm{crit}}\,(\mathrm{nT})={N}_{{{\rm{H}}}_{2}}/(1.57\times {10}^{25}{{\rm{m}}}^{-2})\,\,,\end{eqnarray} \tag{ 8 }$$

when λ = 1. With the above threshold, we obtain B_crit = 25 ± 6 nT for the HAWC+ measurements, B_crit = 61 ± 27 nT for the ALMA data, and B_crit = 42 ± 7 nT averaged over the mapped HAWC+ and ALMA emission in BYF 73, based on the combined HRO analysis.

Again, we can compare this result to equivalent B_crit for the nearby clouds of ∼4 nT (Planck Collaboration 2016), ∼20 nT for Vela C (Soler et al. 2017), and ∼4 nT for L1688 (Lee et al. 2021), placing BYF 73 at a higher-N_crit and -B_crit level than these other clouds. Our B_crit result for BYF 73 generally is also about half the DCF value at the peak of MIR 2 with Mopra's resolution (Equation (2)), which is not unreasonable given the respective density levels. Thus, both the DCF and HRO analyses give us mutually consistent clues about the B field strengths in BYF 73, which seem to be significantly stronger than in other clouds.

5.3.1. Widespread Polarization in the Outflow

The GK effect can arise when B fields (even weak ones) and velocity and excitation gradients in molecular gas combine to produce imbalances from thermal equilibrium in populations of magnetic sublevels M of spectral lines with opacity ∼1 (Goldreich & Kylafis 1981; Girart et al. 2004; Crutcher 2012). This can produce linearly polarized spectral-line emission that is either aligned with (π transitions) or perpendicular to (σ transitions) the local B field, depending on the radiative transfer circumstances, namely the unknown angles between the radiation anisotropy, the line of sight, and the B field direction. In addition, the classical Zeeman effect can give rise to circularly polarized σ transitions parallel to B, observable for that component of B oriented along the line of sight.

We report here the widespread detection of strong, linearly polarized emission in the ¹²CO line wings from BYF 73 (i.e., its bipolar outflow) that is consistent with the GK effect. The native Stokes data have high S/N, up to 21 for $P^{\prime}$ and $p^{\prime}$, in individual 0.16 km s⁻¹-wide channels and across much of the outflow visible in both line wings. However, where the polarization signal weakens, the $p^{\prime}$ values tend to rise, and, with the larger uncertainties, the resulting vector maps become somewhat confusing to look at. Therefore, we averaged (binned) the native Stokes data into ∼3 km s⁻¹-wide channels for display purposes only, and formed the polarization products on the binned cubes with proper noise weighting: the significant polarization features are then more easily visualized in the binned data. Figure 29 shows overlays of the blue- and redshifted I and $P^{\prime}$ emission in these binned data, together with all observed polarization vectors above 4σ_rms (and up to 72σ) across the full velocity range of the line wings.

Zoom In Zoom Out Reset image size

Interpreting GK polarization vectors requires some resolution of the 90° ambiguity described above, depending on whether σ transitions from the M = ±1 magnetic sublevels (polarized perpendicularly to the B field) will be stronger or weaker than the π transitions from M = 0 sublevels (parallel to B). In general, whether outflows are driven by magnetocentrifugal forces anchored in protostellar disks, or by collimated protostellar jets carrying their own B fields, the nominal expectation is that inferred B fields should be aligned along outflows. In Figure 29 we have rotated the vectors by 90° from those observed, and it is this orientation which, remarkably clearly, shows an overwhelming orientation along the outflow direction in each wing, especially for the higher-S/N pixels. Equivalent plots of the observed vectors show a near-universal circumferential alignment around MIR 2, which would seem to be unphysical based on the above understanding.

Indeed, the high-S/N vectors track both the bend in the blue wing and the fork in the red wing inferred solely from the I emission pattern (Figure 8). There are some low-S/N vectors that do not align with this general pattern, however, typically near the $p^{\prime}$ threshold. This is most notable in the −24 km s⁻¹ panel, both north and south of the outflow itself. In the northern portion of this polarized emission, the alignment is instead approximately across the EPL as mapped by HAWC+ (Figures 6–8). South of the outflow, the emission appears to be an artifact of missing short spacings in the field of view, so we discount it.

An arrangement with outflows oriented along the B field direction is typical of Crutcher's (2012) summary of outflow studies via GK mapping. The data for BYF 73 show that the observed polarization is preferentially oriented 90° from the presumed B field direction down the outflow axes, and so supports an excitation condition in which the σM = ± 1 transitions are robustly overpopulated in the outflow relative to the πM = 0 transitions.

5.3.2. Simplified DCF Analysis

Quantifying this description, however, is challenging due to the sheer volume and effectively 4D nature of the data. As a first attempt, we perform a simplified DCF analysis per unbinned 0.16 km s⁻¹ -wide channel in the data. We argue that this is reasonable, even though the original DCF method was not developed for outflows. Indeed, we believe that DCF analysis of spectral-line linear polarization will give a better result than for dust polarization, in the following sense.

A GK-imbued spectral line directly samples the turbulence, density, and polarization dispersion in the same region. A problem with application of DCF to dust polarization is that one needs to estimate the density sampled by the dust polarization, plus a turbulent linewidth. Although we have dust-based column density maps of BYF 73 (Figures 23, 27) from which density estimates can be simply inferred, densities and linewidths are more generally inferred from observations of spectral lines: excitation analysis for density, and directly measured linewidth. However, different spectral lines sample different density regimes, and different lines have different linewidths, so just what is appropriate for the dust polarization analysis is never clear. One often ends up with some sort of ill-defined average along the line of sight. On the other hand, for spectral-line polarization, things are self-consistent. One can infer the polarization dispersion, density, and measure the turbulence directly in the same parcel of gas, namely that sampled by the spectral line being observed. DCF then gives an estimate of the B field strength in that spatial and density parcel, not for some ill-defined and possibly different regions along the line of sight.

If GK polarization can be detected in multiple spectral lines that sample different density regimes, one can in principle build up a 3D picture of the B field. So far, however, detections of the GK effect in species other than CO have been rare, but presumably that will improve as time goes on.

For this channel-DCF analysis of BYF 73, we do not include sub-ROIs of each channel at smaller scales, as in Figures 16 and 19, and assume instead for simplicity that the whole-channel-ROI gives an approximate measure of the ${\theta }_{{B}_{\perp }}$ correlation length for that channel, since the polarized emission is dominated by the outflow structure, as seen in Figure 29. The results are shown in Figure 30 for both ¹²CO line wings.

Zoom In Zoom Out Reset image size

Several features are immediately evident. The most significant are the clear trends in ${\theta }_{{B}_{\perp }}$, for the blue wing from V_LSR = −22 to −36 km s⁻¹, and for the red wing from V_LSR = −17 to −2 km s⁻¹, showing a B field orientation that changes gradually, in both cases, from east–west to more along the outflow axis and then back to east–west, as one looks from the lower to higher outflow speeds. This internal consistency is not so surprising since the statistics in these channels (a few hundred pixels each) are quite robust. Observationally, however, there is no reason to expect the polarization to line up so reliably, channel by channel, unless the polarization signal in all channels is strongly governed by the intrinsic physics of the outflow. Thus, over these velocity ranges, the dispersion in ${\theta }_{{B}_{\perp }}$ for each channel is quite small, s = 11° ± 4°, even where a few pixels appear as outliers in the ${\theta }_{{B}_{\perp }}$ distribution of some channels. This is not much larger than the average noise-derived Δθ_rms ≈ 7°, giving an intrinsic mean dispersion s ≈ ± 86, or as little as ±7° in some places. Overall, the polarized emission at these velocities appears very well organized.

At velocities from V_LSR = −36 to −55 km s⁻¹ and >–2 km s⁻¹, however, the mean ${\theta }_{{B}_{\perp }}$ direction in each channel becomes more erratic as the outflow speed increases, on average still lying near the outflow directions but with a dispersion among the channels s ≈ ± 40° for the blue wing. This wider variation probably reflects poorer statistics, with typically only a few, or a few dozen, pixels per channel.

At velocities closer to the line core, V_LSR = −22 to −17 km s⁻¹, aliasing of extended line emission throughout the field of view introduces many polarization features with probably unreliable ${\theta }_{{B}_{\perp }}$, indicated in Figure 30 by both the increasing density of dots at all ${\theta }_{{B}_{\perp }}$ and higher χ² from the non-Gaussian ${\theta }_{{B}_{\perp }}$ distribution, all becoming more noticeable as the velocity approaches V_sys = –19.6 km s⁻¹.

5.3.3. Magnetic Field Strength Calculations

As with the continuum data (Section 5.1), we can use this basic DCF information on the dispersion in ${\theta }_{{B}_{\perp }}$ per channel from the GK effect, to make estimates of the B field strength in the outflow. For Equation (1) from Section 5.1, we first estimate the ¹²CO column density in the line wing emission ${I}_{{}^{12}\mathrm{CO}}$ via the velocity-resolved, opacity-corrected conversion law from Barnes et al. (2018, their Figure 5(b), not their integrated law in Figure 9(b)). Next, we convert that to an H₂ column density with Pitts & Barnes' (2021) dust temperature-dependent abundance law. Finally, we turn this into a volume density assuming a line-of-sight depth D ≈ 0.087 pc through the outflow (and correlation length, later) equal to the outflow's average projected width:

$$\begin{eqnarray}&&{n}_{\mathrm{ch}}=\displaystyle \frac{{N}_{0}}{D}\,\displaystyle \frac{{({I}_{{}^{12}\mathrm{CO}}{\rm{d}}V/{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1})}^{p}}{{10}^{[-10\,{\mathrm{log}}^{2}({T}_{d}/{T}_{0})+\mathrm{log}{X}_{0}]}}\,.\end{eqnarray} \tag{ 9 }$$

From Barnes et al. (2018), N₀ = 1.27 × 10²⁰ m⁻² and p = 1.92; the ALMA channel width dV = 0.159 km s⁻¹ converts I to the proper units, and from Pitts & Barnes (2021), T₀ = 20 K is the dust temperature at which the gas-phase CO abundance relative to H₂ peaks, at a value X₀ = 0.74 × 10⁻⁴.

Equation (9) thus converts the ${I}_{{}^{12}\mathrm{CO}}$ data cube into a cube of H₂ density per channel at the observed velocity. To combine this with Equation (1), we need a turbulent velocity dispersion. We can choose a velocity FWHM ΔV in the gas corresponding to 1 ALMA channel to be consistent with the above formulation, but in reality it may be several times larger, since the outflow is likely to be turbulent at some level related to the ±25 km s⁻¹ range of flow speeds. In that case, the true velocity FWHM in the gas would re-scale the single-channel B_⊥,DCF estimated via Equation (1) by ψ = (ΔV/dV)^1+p/2, since we would need to evaluate I in Equation (9) over the same ΔV-wide bins (the column density inferred from I, and hence the density, is additive across channels). With p as above, ψ ∝ (ΔV/dV)² approximately, or as the ∼square of the number of channels in a ΔV bin.

So for molecular gas with μ as before and s = 7° in the inner part of the flow, where minimal gas-phase densities inferred from Equation (9) are n_ch ∼ 10⁹ m⁻³ ch⁻¹, we obtain

$$\begin{eqnarray}&&{B}_{\perp ,\mathrm{DCF},\mathrm{outflow}}=0.61\,\mathrm{nT}\sqrt{\displaystyle \frac{{n}_{\mathrm{ch}}}{{10}^{9}{{\rm{m}}}^{-3}}}\left(\displaystyle \frac{\psi }{s/7^\circ }\right)\,\,\mathrm{or}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\quad \approx 4.5\,\mathrm{nT}{\left(\displaystyle \frac{{I}_{{}^{12}\mathrm{CO}}{\rm{d}}V}{10\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{p/2}\left(\displaystyle \frac{\psi }{s/7^\circ }\right)\,\,\end{eqnarray} \tag{ 11 }$$

as a minimum for B in the molecular outflow, when approximating the dust temperature at 20 K throughout to compute a minimal H₂ density at the peak CO abundance.

As described above, one channel probably slices the turbulent structure in the outflow rather finely: for ΔV = 3 km s⁻¹, for example, the scaling would increase the single-channel B_⊥,DCF coefficient by ψ = 320×, e.g., to 1.4 μT for the same value of I in Equation (11). Indeed, the brightness of the mapped ¹²CO outflow emission in Figure 29 ranges up to 90 K km s⁻¹ in 3 km s⁻¹-wide bins, suggesting that B fields might be stronger still in some locations. Of course, this scaling may not actually be valid: while simulations of turbulent plasmas suggest that DCF estimates are reasonable up to Mach numbers M ∼ 5–9 (Ostriker et al. 2001), such values may be far exceeded ( M > 100) in the outflow.

5.3.4. Energy Densities

Despite these somewhat large uncertainties, we at least have rough estimates for B_⊥,DCF field strengths in the outflow. As a final exercise, we compare the energy density ${\mathfrak{M}}$ that would exist in such B fields with the kinetic energy density ${\mathfrak{K}}$ of the outflowing gas. ${\mathfrak{M}}$ follows directly from Equation (10),

$$\begin{eqnarray}&&{\mathfrak{M}}={B}^{2}/2{\mu }_{0}=1.49\times {10}^{-13}{\mathrm{Jm}}^{-3}\,{n}_{9}{\left(\displaystyle \frac{\psi }{s/7^\circ }\right)}^{2},\end{eqnarray} \tag{ 12 }$$

where μ₀ is the permeability of free space (or the magnetic constant in preferred SI usage) and n₉ = n_ch/10⁹m⁻³. Thus, smaller values for ${\mathfrak{M}}$ will be obtained at lower gas densities (and approximately I, through Equation (9)), either at the edges of the outflow or at higher velocities, whereas larger ${\mathfrak{M}}$ will lie closer to the outflow origin where the density or I is higher.

For ${\mathfrak{K}}$ we can assume a cylindrical outflow geometry (diameter D, length L) to approximately compute

$$\begin{eqnarray}\begin{array}{rcl}{\mathfrak{K}} & = & \displaystyle \frac{\tfrac{1}{2}{{MV}}_{\mathrm{rel}}^{2}}{\tfrac{1}{4}\pi {D}^{2}L}=\displaystyle \frac{1}{2}\rho {V}_{\mathrm{rel}}^{2}\\ & = & 3.9\times {10}^{-12}\,{\mathrm{Jm}}^{-3}\,{n}_{9}\,{\left({V}_{\mathrm{rel}}/\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}\end{array}\end{eqnarray} \tag{ 13 }$$

where L ≈ 0.6 pc is the physical length of the 50''-long blue lobe of the outflow. The second expression is much simpler to use, since the mass density ρ = 2μ m_H n_ch, with the number density n_ch from Equation (9). This can actually be done separately for each channel if we use its relative outflow speed V_rel as measured from V_sys = −19.6 km s⁻¹. Thus, the value of ${\mathfrak{K}}$ will be larger at higher ρ ∝ I (approximately) but especially at higher V_rel, or smaller at lower I or especially lower V_rel.

We can now compute a datacube of the ${\mathfrak{M}}$/${\mathfrak{K}}$ ratio,

$$\begin{eqnarray}&&\displaystyle \frac{{\mathfrak{M}}}{{\mathfrak{K}}}=3.8 \% \,\displaystyle \frac{{n}_{9}{\left(\tfrac{\psi }{s/7^\circ }\right)}^{2}}{{n}_{9}{({V}_{\mathrm{rel}}/\mathrm{km}\,{{\rm{s}}}^{-1})}^{2}}\,,\end{eqnarray} \tag{ 14 }$$

which turns out to be independent of the density as long as we measure both over the same channels or velocity bins. Since ψ is part of the density scaling, this simplifies to

$$\begin{eqnarray}&&\displaystyle \frac{{\mathfrak{M}}}{{\mathfrak{K}}}=3.8 \% \,{\left(\displaystyle \frac{{\rm{\Delta }}V/{V}_{\mathrm{rel}}}{s/7^\circ }\right)}^{2}\,,\end{eqnarray} \tag{ 15 }$$

which we can evaluate per native channel of width dV (or any other binning), even while using a larger ΔV to represent the turbulence in the flow. In other words, the ${\mathfrak{M}}$/${\mathfrak{K}}$ ratio can reasonably be estimated with some knowledge of only the gas turbulence ΔV and polarization dispersion s in the B field directions in each channel at V_rel, which is ultimately just a restatement of the DCF method, per unit volume.

For purposes of illustration, we take ΔV = 3 km s⁻¹ and a more conservative s = 13°, and present the ${\mathfrak{M}}$/${\mathfrak{K}}$ ratio results in Figure 31 at some representative channels V_rel. Different ΔV or s values would obviously scale the ratios as (ΔV/s)². Despite the larger s and smaller ${\mathfrak{M}}$ in Figure 31 than discussed above, ${\mathfrak{M}}$/${\mathfrak{K}}$ peaks at 37, i.e., ≫1. We discuss this further in Section 6.

Zoom In Zoom Out Reset image size

As the only observational technique capable of directly measuring B field strengths, the Zeeman effect has been widely utilized over five decades (Crutcher 2012). However, successful detections are notoriously difficult: for extended thermal emission from molecular clouds, only HI, OH, and CN have yielded B field detections, and among these, only CN can provide information on field strengths in dense (≳10¹¹m⁻³) gas. Despite considerable effort, there still exist only 14 individual CN-Zeeman measurements from a heterogeneous sample of clouds (Falgarone et al. 2008). But with the advent of full-polarization capability in Cycle 7, anticipation has been high that ALMA might fundamentally change the state of play in this field.

Unfortunately, despite the very high S/N (∼200) in the ALMA Stokes I data for BYF 73, covering eight of the nine hyperfine transitions of the CN J = 1 → 0 line, the V cube shows nothing discernible above the noise. Computing the ratio of Stokes V to dI/dV and scaling this to the Zeeman splitting coefficient of any of the brightest hyperfine transitions (as in Table 1 of Falgarone et al. 2008) yields only 3σ limits ∼1 μT (10 mG), as seen in Figure 32. This is near the upper end of the range of field strengths seen before in dense gas (Crutcher 2012), but the noise level would need to be at least halved to obtain reliable measurements even at those levels. A further issue was the Cycle 7 limit for accurately calibrated V data lying within the inner 10% of the primary beam.

Zoom In Zoom Out Reset image size

This also means we cannot use Zeeman data to distinguish between the scenarios (a pure B_⊥-twist or an additional B_∣∣ component) put forward to explain the HAWC+ $P^{\prime}$ null on the western edge of the IBL (Section 3.3).

While somewhat discouraging, the nondetection may partly be due to the cloud's orientation. That is, the Zeeman effect can only measure the line-of-sight component B_∣∣. The fact that the outflow is viewed close to side-on (Section 4.1) suggests that most organized structures in the molecular cloud, such as an accretion disk around MIR 2, would also probably be presented edge-on to us, as might any structures being accelerated away from it, thus possibly maximizing B_⊥ and minimizing B_∣∣.

Based on 40'' resolution Mopra HCO⁺ maps, Barnes et al. (2010) first described a massive infall of dense, cold material within the wider BYF 73 cloud, without any evidence of an outflow characteristic of lower-mass YSOs. This suggested an extremely early evolutionary state for a very massive protostar, which seemed to be confirmed by the mid-IR and FIR data of P18. With the higher-resolution SOFIA and ALMA data presented here, particularly the strong bipolar ¹²CO outflow, we see that the original appearance of outflow-free, extremely young massive SF may have been something of a masquerade. ¹⁵  Nevertheless, through the ALMA ¹³CO data, we are able to discern more specific clues to the configuration of the inflow originally seen in HCO⁺.

However, the 3D relationship between the Streamer, outflow, and disk or infall as described in Section 4 remains puzzling. The disk can be traced from an outer radius of 0.18 pc = 36,000 au to an inner radius no larger than the limit of the ALMA resolution, 18 = 4500 au. Apparently, this disk is close to edge-on based on the sharp velocity gradient across MIR 2, so the filamentary impression of the Streamer may be an illusion. For example, we see that both the outflow and rotational/infall patterns are oriented east–west, but this arrangement would seem physically counterintuitive. Undoubtedly, there is some depth to these features in the line of sight, and it is possible that any inflow to the disk might be approaching MIR 2 from behind its eastern side, even while the blue jet is receding from MIR 2 on the same side. Likewise, inflow overlying the western Streamer may be from the front, while the red jet recedes from MIR 2 as it encounters the H ii region. Separating these features in the line of sight requires only a few ×10⁴ au ∼0.1 pc, so we consider this scenario reasonable. Moreover, the Streamer/disk is clearly not flat; there is a measurable width and curvature to its structure.

What are the dimensions of the disk/infall zone centered on MIR 2? The modal value of the latitude across all of the disk emission, from the LV-1st moment map in Figure 12 (middle panel), is b = 016997 ± 000005, only 08 = 0.3 ALMA beamwidths north of MIR 2 itself. Three-quarters of this emission lies within 0169 < b < 0172, a span of only 11'' or ∼4 ALMA beams, strongly suggesting a somewhat narrow structure for the high-ΔV material. We also computed an LV-2nd moment map (Figure 12, bottom panel) to examine the latitude width, confirming that it is indeed thin, from only ∼3 ALMA beams = 7'' thick to <1 beam. From an inspection of all three LV-moment maps, we see that the eastern side of the disk at the higher velocities (V_LSR < −25 km s⁻¹) seems to lie mostly at one latitude close to that of MIR 2, b = 01698 ± 00005, and so is indeed quite flat in the east–west plane to within 1 ALMA beamwidth. This is across an extent of 001 = 36'', for an aspect ratio of 15–20:1 oriented east–west.

Given this, it is hard to imagine a disk oriented in the same direction as the outflow it is supposed to be driving. This favors the ¹³CO data tracing free-falling material onto a 950 M_⊙ MIR 2 within a 36'' × 2'' structure, rather than a Keplerian disk, since such a disk ought to be oriented close to north–south, parallel to the sharpest velocity gradient across MIR 2. If this is indicative, the gradient suggests a disk thickness perhaps ≲2'' or 5000 au, but possibly even narrower. On the other hand, even if we separate the infall from the outflow along our line of sight, it is equally hard to see how a predominantly east–west infall (i.e., along a polar direction) produces a disk oriented north–south. So the puzzle persists.

Additionally, if both the MIR absorption and FIR emission mass estimates at MIR 2's peak position are 5×–12× too small (P18), this is possible evidence for significant grain growth in MIR 2's protostellar envelope; any free–free emission from MIR 2 (Section 3.1) would make this discrepancy worse. P18's gravitational energy release luminosity also scales with the mass, raising it to perhaps 20%–33% of MIR 2's total luminosity. Indeed, if future higher-resolution observations revealed an impact radius for the inflow only five times smaller at 1000 au, this could not only account completely for MIR 2's brightness via gravitational energy release, but also possibly reveal it to be the first example of a massive "first hydrostatic core."

At velocities closer to systemic, the brightest LV emission in the eastern disk, corresponding to the Streamer at 28622 > l > 28621, also lies very close to this east–west plane, although slightly north of it. However, it is not distributed along any of the Kerplerian curves: instead, its velocity drops linearly from −23 km s⁻¹ to systemic over its 36'' length, with kinematics mimicking that of solid-body rotation. This portion of the eastern disk is thicker, 5''–7'', than the high-velocity emission there, <3'', giving it an aspect ratio around 6:1. Continuing the non-Keplerian behavior, east of l = 28622 or outside the Streamer's distance from MIR 2, there appears to be some material in "superrotation" in the bottom-left quadrant of the curves, i.e., with V_LSR exceeding the rotational curve for 1350 M_⊙. This lies at b = 0172 (red in the middle panel of Figure 12), or 7'' north of MIR 2, but is again about as thin (1 ALMA beam) as the high-velocity disk material. However, the envelope of this material's superrotation is moving at close to $\sqrt{2}\times$ this curve, suggesting either free-fall of ambient material toward the Streamer/disk from the rear of the cloud, or that the enclosed mass at this radius has ∼doubled.

In contrast, the western side of the disk seems to curve somewhat north of the east–west plane of the eastern disk, to a latitude as far as 12'' north of MIR 2 at 0173, and at a moderately high velocity (V_LSR = −10 km s⁻¹) from systemic. The rest of the western disk ranges in latitude from MIR 2's value up to this limit, and the line of maximum velocity in the redshifted wing map (Figure 11, right panel) is clearly curved to the northwest from MIR 2. The western disk's thickness (Figure 12) is also broader than for the eastern disk overall, up to 7'', but is also thin (<3'') in many places, even where it curves to the northwest. It is worth noticing that the solid-body portion of the Streamer/eastern disk seems to continue partway (0006 = 22'') into the western disk in both the 1st- and 2nd-LV-moment maps, but then seems to reverse bluntly back to MIR 2's longitude at V_LSR = −16 km s⁻¹, while thickening to a width of almost 10'' just west of MIR 2, almost as if the Streamer's infall (if that's what it is) were being deflected from the east–west plane by some obstacle west of MIR 2.

What of the counterrotating parts of the LV diagrams (i.e., the "empty" top-left and bottom-right quadrants of the rotation curves)? Much of this emission, especially the brighter portions thereof, lie north (b > 0172, magenta) or south (b < 0168, black) of the disk, and outside the longitude range of the IBL (286216 > l > 286203): they appear to be associated with other internal structures of the cloud, supporting the rotational interpretation for the inner parts of the Streamer.

While MIR 2's mass seems dynamically dominant within the IBL, the mass in the Streamer/disk must nevertheless also be significant. In Section 5.2 we found a mean column density 3 × 10²⁷ m⁻² ≈ 6000 M_⊙ pc⁻² along the Streamer, or roughly 15 M_⊙ per 4'' box, assuming there is not the same mass deficit/degree of grain growth as in the MIR 2 core. A rough total along the full 72'' length and 8'' width of the streamer is then perhaps 500 M_⊙. This would explain why MIR 2's gravitational influence seems to drop beyond the outer Keplerian radius determined above, since there, the gas mass of the broader cloud starts rivaling MIR 2's effects.

In such a disk, the 0.034 M_⊙ yr⁻¹ mass accretion rate determined by Barnes et al. (2010), still a record as far as we know, can be supported by a merely 0.01%/yr "leakage" of mass through the disk onto MIR 2's core, or alternatively, that the Streamer can supply this accretion rate for another 10⁴ yr. On the other hand, the time required to build up the more massive MIR 2 core at this accretion rate is closer to 40,000 yr instead of the 7000 yr estimated by P18, assuming that Barnes et al.'s (2010) accretion rate is correct. Without detailed modeling, we cannot refine the accretion rate value here beyond an approximate calculation below, but the true rate seems unlikely to be too much less than this, with such a massive reservoir available for accumulation.

As one example, consider the velocity difference between the brightest disk material within the bottom-left quadrant of the rotation curves, and the 1350 M_⊙ curve itself, as seen in Figure 12 between 286217 ≳ l ≳ 286208 and at latitudes ∼4'' north of MIR 2. This difference runs from 0 km s⁻¹ at the eastern end of this window to ∼+10 km s⁻¹ at MIR 2, in the sense of being "subrotational" in our line of sight. If we suppose that the rotational motion here is being translated into proper motions inward to MIR 2, the effective accretion speed can be taken (very approximately!) as V_accr ∼5 km s⁻¹ = 5 pc Myr⁻¹. The emission along this feature (which covers perhaps half of the main part of the 8''-wide Streamer) averages ∼10 K/ch in brightness, or conservatively, ∼20 K km s⁻¹ integrated.

Using a simple conversion factor X(¹³CO) = 10²⁶ m⁻² (K km s⁻¹)⁻¹ (probably an underestimate; Barnes et al. 2018), the column density in this feature alone runs around N_accr ∼ 2 × 10²⁷ m⁻² = 4000 M_⊙ pc⁻², or perhaps $\tfrac{2}{3}$ of the Streamer's total column density as seen in Figure 27. This translates to a linear density Λ_accr ∼ 400 M_⊙/pc within the 0.1 pc width of the presumed accretion stream. Therefore we have a mass flux in this accretion stream of ${\dot{M}}_{\mathrm{accr}}$ = Λ_accr V_accr ∼ 2 × 10⁻³ M_⊙ yr⁻¹.

This very rough estimate is still on the large side compared to other massive protostars (Rygl et al. 2013), but smaller than Barnes et al.'s (2010) rate of 0.034 M_⊙ yr⁻¹. The true value is likely a multiple of the above example, however, due to several factors: our conservative starting column density, other accretion flows such as the western side of the Streamer, the superrotating material, higher-density streams traced better by C¹⁸O or CN, two-times-faster flow closer to MIR 2, and a probably five-times-larger effective X(¹³CO). Thus, the Barnes et al. (2010) rate may still be a reasonable global estimate.

In summary, the complex yet potentially understandable structure of the Streamer near MIR 2, and Keplerian disk/free-falling infall zone around it, may have much to tell us about heavy mass accretion onto a massive protostar. Clearly, MIR 2 is an exceptional and exciting object that demands further study.

The ¹²CO outflow from MIR 2 is fairly massive: with a typical line brightness ${I}_{{}^{12}\mathrm{CO}}$ ∼ 10–30 K per 0.16 km s⁻¹ ALMA channel or average integrated intensity ∼200 K km s⁻¹, we can use Equation (9) or its ilk (Barnes et al. 2018) to estimate gas column densities of about 8 × 10²⁵ m⁻² or 160 M_⊙ pc⁻² in the outflow. Inside the flow dimensions of roughly 1 × 0.1 pc, this gives a total outflow mass of perhaps 16 M_⊙. This mass is being driven to speeds of tens of kilometers per second, so the kinetic energy of the flow is similarly large, about 1.6 × 10³⁹ J. If this emerges over timescales of tens of thousands of years, then the mass outflow rate is ${\dot{M}}_{\mathrm{out}}$ = Λ_out V_out ∼ 5 × 10⁻⁴ M_⊙ yr⁻¹ (or ∼10% of the infall rate as estimated above, similar to other outflows; Pudritz & Ray 2019) and the mechanical luminosity of the outflow L_out ∼4 L_⊙. Could this mechanical power be imparted by B fields?

From Equation (15) and Figure 31, we see that the energy density in the B field is typically well below the kinetic energy density in the higher-velocity and lower-density gas: the B field is therefore likely a passenger in the flow at these points. On the other hand, ${\mathfrak{M}}$ may rival or even exceed ${\mathfrak{K}}$ where V_rel is small. This result, however, should be taken as merely suggestive, since ${\mathfrak{M}}$/${\mathfrak{K}}$ > 1 only in the 20 lowest-V_rel channels, where the ¹²CO opacity is still high and we may not be mapping much of the outflow via ¹²CO. But B fields do seem to be detected throughout the outflow, at a few × 10 nT. It seems reasonable to suppose that similar B fields (at least!) should exist close to BYF 73's V_sys, and specifically close to the base of the flow at MIR 2. If this were true, the B field would at least have the potential of being energetically important.

As such, this is circumstantial yet valuable evidence that the B field may be intimately involved in driving, or at least shaping, the outflow. Pudritz & Ray (2019) showcased some other recent observational results making this B field connection to the outflow, typically on ∼100 au (i.e., disk) scales. As far as we are aware, this is the first instance where the structure of the whole molecular outflow might at least partially be attributable to the B field at its origin. The connection is not always clear, however: a rare case where the B field seems to play an important role in massive SF is the compact H ii region K3-50, where a strong ionized outflow emanates from a high-mass protostellar object, surrounded by a Keplerian disk extending over radii 0.1–0.7 pc (Howard et al. 1997). There, the B field inferred at the inner edge of the disk seems to be strong enough to provide support against gravity; however, even in K3-50, the B field does not seem strong enough to influence the outflow (Barnes et al. 2015).

Our results for BYF 73 seem to provide additional observational support for the picture of magnetocentrifugally powered protostellar outflows (Shu et al. 1994; Ouyed & Pudritz 1997; Tomisaka 2011), as opposed to the main competing model of turbulent entrainment of gas from a bipolar jet (e.g., Raga et al. 1993). Recent numerical work on solar coronal mass ejections (Jiang et al. 2021) may even supply a specific mechanism for the high speeds in such outflows, namely sudden magnetic reconnection in bipolar loops, presumably anchored in the inner parts of a protostellar accretion disk. It is tempting to suppose such reconnecting loops drive the vigorous outflows widely seen in other star-forming clouds, as may be happening here with MIR 2/BYF 73. Future work in this area, supported by ALMA+SOFIA observations, should be very interesting.

We briefly also note that the critical N, n, and B_⊥ values derived from HRO analysis of the SOFIA+ALMA data (Section 5.2) place BYF 73 just below the maximum B_∣∣ trend line in Crutcher's (2012) n versus B summary plot (his Figure 6), nicely among other dense clouds' CN-Zeeman measurements. In contrast, BYF 73 is slightly above the line of criticality in Crutcher's (2012) N versus B_∣∣ plot (his Figure 7), on the side of being subcritical and a little above its supercritical counterparts in this regard. Given the uncertainties in our results, however, this may not be terribly significant. Further, BYF 73 is not the most extreme strong-B outlier compared to the line of criticality, but it may be the highest column density subcritical cloud. As Crutcher (2012) points out, this would be unusual in the sense of at least supporting the possibility of ambipolar diffusion playing an important role in cloud stability (Mouschovias & Ciolek 1999). However, our HRO result is not the same as a direct Zeeman strength measurement, and it should probably not yet be overinterpreted without some confirming evidence.

Nevertheless, it is tempting to wonder what higher-resolution and -sensitivity observations might reveal about the material closer in to MIR 2, where the infall might be more clearly imaged, the outflow may be driven, and its originating disk might be discerned.