The Magnetic Field in the Milky Way Filamentary Bone G47

G47 was observed with SOFIA HAWC+ in Band E, which is centered at 214 $\mu$m and provides a resolution of 18$\farcs$2 (Harper et al., 2018) or 0.58 pc resolution at a distance of 6.6 kpc. The observations were taken over multiple flights in September 2020 during the OC8E HAWC+ flight series as part of the FIELDMAPS legacy project. The polarimeter’s field of view is 4$\farcm$2 $\times$ 6$\farcm$2. The entire bone was mapped by mosaicking together four separate on-the-fly (SCANPOL) maps. The total time on source for the combined observations was 4070 s. We use the Level 4 delivered products from the SOFIA archive.¹¹1https://irsa.ipac.caltech.edu/applications/sofia/ The pixel sizes are $3\farcs 7\times 3\farcs 7$, which oversamples the 18$\farcs$2 beam. Errors along the bone for the Stokes $I$, $Q$, and $U$ maps varied from about 0.5 to 0.8 mJy pixel^-1. From the Stokes parameters, the polarization angle, $\chi$, at each pixel is calculated via

where the arctan2 is the four-quadrant arctangent. The positively biased polarization fraction, $P_{\text{frac,b}}$, at each pixel is calculated via

Polarization maps have been de-biased in the pipeline via $P_{\text{frac}}$ = $\sqrt{{P_{\text{frac,b}}-\sigma_{P_{\text{frac}}}^{2}}}$, where $\sigma_{P_{\text{frac}}}$ is the error on $P_{\text{frac,b}}$ (and $P_{\text{frac}}$).

SOFIA HAWC+ is not sensitive to the absolute Stokes parameters, so some amount of spatial filtering via on-the-fly maps affects the data. These can lead to artifacts that show up as unrealistically large $P_{\text{frac}}$ values. However, large $P_{\text{frac}}$ values ($>$20%) are outside of the areas of interest in this paper, and thus will not be used for any analysis.

Delivered data were in equatorial coordinates, and we rotated them to Galactic coordinates via the python package reproject (Robitaille et al., 2020) and properly rotating position angles (Appenzeller, 1968).

We use the $Herschel$ 250 $\mu$m continuum data from Hi-GAL (Molinari et al., 2016), and H₂ column density maps ($N_{\text{H$_{2}$}}$) generated by Zucker et al. (2018) from the multi-wavelength Hi-GAL $Herschel$ data. These were subsequently converted to $N_{\text{gas}}$ (i.e., $N_{\text{H$_{2}$}}$ + $N_{\text{He}}$) by multiplying by the ratio of the mean molecular weight per H₂ molecule ($\mu_{\text{H$_{2}$}}$ = 2.8) divided by the mean molecular weight per particle ($\mu_{p}$ = 2.37; Kauffmann et al. 2008). Zucker et al. (2018) also fit the “spine” of the bone – equivalent to a one-pixel wide representation of its plane-of-the-sky morphology – using the RADFIL algorithm (Zucker & Chen, 2018). The resolution of the column density and spine maps are 43$\arcsec$ ($\sim$1.4 pc), and they have pixel sizes of 11$\farcs$5 $\times$ 11$\farcs$5.

We also use ¹³CO(1–0) data from the Galactic Ring Survey (GRS; Jackson et al., 2006) and NH₃(1,1) data from the Radio Ammonia Mid-plane Survey (RAMPS; Hogge et al., 2018), each of which we convert to velocity dispersion maps via Gaussian fits following Hogge et al. (2018). While ¹³CO(1–0) is detected everywhere along the bone, NH₃(1,1) is only detected toward the densest parts. We make a “final velocity dispersion” map, which we will use to estimate B-field strengths, where we use the NH₃(1,1) velocity dispersion when it is available for a particular pixel and ¹³CO(1–0) otherwise. We combine the two together in which we use NH₃(1,1) velocity dispersion when it is available for a particular pixel; otherwise we use ¹³CO(1–0). In the dense regions where NH₃(1,1) is detected, ¹³CO(1–0) linewidths tend to be higher (factor of $\sim$2), as it is the combination of diffuse and compact emission. Locations where NH₃(1,1) is not detected are expected to be more diffuse, and thus ¹³CO(1–0) widths are mostly accurate in these areas.

Locations of the Class I and II YSOs were taken from Zhang et al. (2019), which were identified via $Spitzer$ observations. Zhang et al. (2019) estimated the survey completeness for Class I YSOs to be a few tenths of a solar mass and for Class II YSOs to be a few solar masses. Class I YSOs are likely at locations of the highest star formation activity along the bone. Xu et al. (2018) identified several more YSO candidates toward G47, but unlike Zhang et al. (2019), they did not use criteria to exclude contaminants such as AGB stars.

To estimate how the B-field changes across the bone, we apply the DCFST technique along the bone’s spine. We do this by “sliding” a rectangular box down the spine, allowing the box to rotate as the bone’s spine change directions in the sky. The center of the box changes one column density/spine pixel at a time (one spine pixel is 11$\farcs$5 $\times$ 11$\farcs$5), and the PA of the rectangle is given by the instantaneous slope of the 24${}^{\text{th}}$ ordered polynomial fit the spine, as discussed in Section 3. The sliding rectangular box has a width, $w$, of 20 HAWC+ pixels (74$\arcsec$) and a height, $h$, of 15 HAWC+ pixels (55$\farcs$5), equivalent to a width and height of $\sim$4 and 3 HAWC+ beams, respectively. These dimensions allow for just over 10 independent beams for each box, which is a sufficient amount of data points for calculating the angular dispersion for the DCFST technique. We do not use a larger box since our underlying assumption is a uniform field in each box, and the field becomes less uniform at larger scales, resulting in an overestimate of the angle dispersion.

We create four image cutouts for each rectangular box sliding down the spine: one for the B-field PA map and another for its error map, one for final velocity dispersion map, and one for the column density map (see Section 2.2 for discussion of the latter two). Appendix A.1 discusses how to determine whether a given pixel of a map is located within the sliding box. To apply the DCFST technique (Equation 3), we need to estimate $\bar{\rho}$, $\delta v_{\text{los}}$, and $\delta\theta$ for each rectangular box. In the column density cutout, we take the median value as our measure of the mean column density,²²2This removes potential outliers. The percent difference between the mean and median for each box is typically less than 5% and never more than $\sim$10% $\bar{N}_{\text{gas}}$, and subsequently convert it to a number density $\bar{n}_{\text{gas}}$ and then $\bar{\rho}$ assuming a cylindrical bone (see Appendix A.2). $\delta v_{\text{los}}$ was chosen to be the median value in the final velocity dispersion cutout. From the B-field PA cutout, we calculate the standard deviation of the cutout, $\delta\theta_{\text{obs}}$, and for its error cutout, we take the median value, which we call $\sigma_{\theta}$. The estimated intrinsic angle dispersion, $\delta\theta$, can be corrected for observational errors such that $\delta\theta=\sqrt{\delta\theta_{\text{obs}}^{2}-\sigma_{\theta}^{2}}$. From this we can estimate the plane-of-sky B-field strengths, $B_{\text{pos}}$. We do not calculate the B-fields at locations with $\delta\theta>25^{\circ}$ (Ostriker et al., 2001) since then the turbulence driving the angular dispersion would be super-Alfvénic.

The $B_{\text{pos}}$ map is shown in the left panel of Figure 4, with fields ranging from $\sim$20 to 160 $\mu$G. We then solve for the critical ratio, $\lambda$, by taking the ratio of the observed and critical mass to magnetic flux ratios, i.e.,

(Crutcher et al., 2004). $\lambda$ parameterizes the relative importance of gravity and magnetic fields. For $\lambda>1$, the gas is unstable to gravitational collapse, and when $\lambda<1$, fields can support the gas against collapse.

The observed mass within the $w\times h$ box is $M_{\text{observed}}$ = $\mu_{p}$$m_{p}\bar{N}_{\text{gas}}wh$. The magnetic flux is calculated via Appendix A.3, and we approximate $(M/\Phi)_{\text{crit}}$ as $1/(2\pi\sqrt{G})$ (McKee & Ostriker, 2007). The $\lambda$ map is shown in the right panel of Figure 4. Note that errors on these values are difficult to quantify given that some input parameters have non-Gaussian errors, and we make assumptions about the geometry of the bone. These values of $B_{{\text{pos}}}$ and $\lambda$ reflect our best guesses from the data, and we expect them to be correct within a factor of 2–3. However, since many of the uncertainties in our results are not dominated by random effects but are correlated, for example via the column density or geometrical assumptions, the relative change in these parameters along the bone are likely to be more accurately determined.

Along the spine of the bone, there are two main groups of YSOs: one toward the northeast and one toward the southwest. At both these locations, $\lambda$ tends to be close to or larger than one, indicating that the areas are typically supercritical to collapse. On the other hand, there are several areas along the bone where $\lambda<1$, indicating that B-fields are potentially strong enough to resist local gravitational contraction. For each position along the spine for which we have a measurement of the B-field, Figure 5 shows $B_{\text{pos}}$ and $\lambda$ as a function of the average number density, $\bar{n}_{\text{gas}}$. The figure also indicates whether most of the velocity dispersion pixels in the sliding box are based on ¹³CO or NH₃(1,1) since this tracer governs the median velocity dispersion. For densities $\bar{n}_{\text{gas}}$ $\lesssim$ 1700 cm^-3 ($\bar{N}_{\text{gas}}$ $\lesssim$ $1.6\times 10^{22}$ cm^-2), $\lambda$ is typically less than 1 (subcritical), while for higher densities (locations of most YSOs), $\lambda$ is typically higher and often supercritical. Overall, there is little change in $B_{\text{pos}}$ as a function of $\bar{n}_{\text{gas}}$. However, if we only consider the field strengths where NH₃ primarily traces the velocity dispersion, i.e., the black points in Figure 5, the field strength increases slightly as a function of density. Based on the linear regression fit to these points, the slope is 0.022 $\pm$ 0.006 $\mu$G/cm^-3. However, given the dispersion of points over a small range of densities, we cannot draw conclusions from this relation.

Together, these results indicate that the field in some parts of the bone can support against collapse, while in other areas, it is insufficiently strong and thus the gas collapses to form stars. Since YSOs are forming in areas where $\lambda$ is equal to or less than one, this indicates that either we are underestimating $\lambda$ or that the high values of $\lambda$ are more localized to YSOs and would necessitate higher resolution polarimetric observations for proper measurement.

We note that for the $B_{\text{pos}}$ vs $\bar{n}_{\text{gas}}$ panel, there are 3 points with higher field strengths ($>$100 $\mu$G) than others. These points are sequentially located next to each other along the spine (see Figure 4). While at these scales G47 has mostly one main velocity component, at this location there appears to be potentially two velocity components that causes the velocity dispersion (and thus the B-field strength) to be overestimated by a factor of $\sim$2. Nevertheless, our fitting routine finds that a one component is slightly better than two, so we only consider it as one component.

The column density peaks toward the southwest show a pinched B-field morphology. Field lines that are frozen to the gas can create such pinched morphology during collapse. Mestel (1966) and Mestel & Strittmatter (1967) calculated the B-field distribution via non-homologous spherical collapse assuming flux-freezing. Myers et al. (2018, 2020) extended these calculations for a uniform field collapsing to Plummer spheroids. Since the southwest peaks of G47 have two column density peaks,³³3SOFIA 214 $\mu$m and $Herschel$ 160/250 $\mu$m maps resolve the top core into two more cores, but for simplicity and to perhaps better reflect the initial collapse, we consider them as one core. we apply this technique using two Plummer spheroids. We first rotate the delivered G47 data (i.e., in equatorial coordinates) clockwise by 12^∘ to align the column density peaks in the up-down direction. We then fit the column density maps with two $p=2$ Plummer spheroids. With the assumption of an initial uniform field and flux freezing, we can use the resulting Plummer spheroids to predict the field morphology by summing the contributions to horizontal and vertical B-field components from each spheroid (see Sections 2 and 3 of Myers et al. 2020). The resulting column density and plane-of-sky B-field lines for the model are shown in the top left panel of Figure 6.

We overlay the model on top of the HAWC+ polarimetric map (Figure 6, top right panel). We mask out angles where the uncertainty in the angle is $>$5^∘. We also mask out angles that differ from the model by more than 25^∘ since these angles are poorly described by the SFF model (i.e., they are outliers that make the distribution non-Gaussian), and these areas may harbor systematic gas flows not included by the model. We calculate the difference in angles between the unmasked inferred field directions and the model (bottom left panel). The dispersion of this distribution is $\delta\theta=11^{\circ}$.⁴⁴4Observational errors, i.e., $\sigma_{\theta}$, in this area are typically only $\sim$2 degrees and thus do not significantly affect $\delta\theta$., which is $\sim$75% of the value of $\delta\theta$ if we did not apply a model. From these data, we can again apply the DCFST technique on the unmasked area. The area is 5.4 pc² in size and the median velocity dispersion based on NH₃ data is $\delta v_{\text{los}}=0.8$ km s^-1. From our model, we find an average number density of 4300 cm^-3.⁵⁵5The number densities of the SFF analysis are slightly higher than that of the sliding box analysis because the line of sight path length for the SFF model is the width of masked region, which is smaller than the diameter of the bone. From these values, we calculate a mean field in the plane of sky in the unmasked area of $B_{\text{pos}}$ = 56 $\mu$G. The peak total field strength for the SOFIA beam is then $B_{0}=108\,\mu$G, assuming $\bar{B}_{tot}=(4/\pi)B_{{\text{pos}}}$ (Crutcher et al., 2004) and $(B_{0}/\bar{B})=(n_{0}/\bar{n})^{2/3}$. The total mass in this region is 1170 $M_{\odot}$, resulting in a mass-to-flux parameter of $\lambda=1.7$. The sliding box analysis along the spine of this area found comparable B-field strengths in this region of $\sim$30–75 $\mu$G with values of $\lambda$ between 0.8 and 1.4. We note the sliding box is $\sim$70% of the size of the unmasked region. These mass-to-flux ratios lie within the range of mass-to-flux ratios in low-mass star-forming cores, according to a recent study (Myers & Basu, 2021).