The Magnetic Field in the Milky Way Filamentary Bone G47

G47 was observed with SOFIA HAWC+ in Band E, which is centered at 214 μm and provides a resolution of 182 (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 Band E field of view is 42 × 62. 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. ¹⁹ The pixel sizes are 37 × 37, which oversamples the 182 beam. Errors along the bone for the Stokes I, Q, and U maps varied from about 0.5–0.8 mJy pixel⁻¹. From the Stokes parameters, the polarization angle, χ, at each pixel is calculated via

$$\begin{eqnarray}&&\chi =\displaystyle \frac{1}{2}\arctan 2\left(\displaystyle \frac{{\rm{U}}}{{\rm{Q}}}\right),\end{eqnarray} \tag{ 1 }$$

where the arctan2 is the four-quadrant arctangent. The positively biased polarization fraction, P_frac,b, at each pixel is calculated via

$$\begin{eqnarray}&&{P}_{\mathrm{frac},{\rm{b}}}=\displaystyle \frac{\sqrt{{Q}^{2}+{U}^{2}}}{I}.\end{eqnarray} \tag{ 2 }$$

Polarization maps have been de-biased in the pipeline via P_frac = $\sqrt{{P}_{\mathrm{frac},{\rm{b}}}-{\sigma }_{{P}_{\mathrm{frac}}}^{2}}$, where ${\sigma }_{{P}_{\mathrm{frac}}}$ is the error on P_frac,b (and P_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_frac values. However, large P_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 μm continuum data from Hi-GAL (Molinari et al. 2016), and H₂ column density maps (${N}_{{{\rm{H}}}_{2}}$) generated by Zucker et al. (2018) from the multi-wavelength Hi-GAL Herschel data. These were subsequently converted to N_gas (i.e., ${N}_{{{\rm{H}}}_{2}}$ + N_He) by multiplying by the ratio of the mean molecular weight per H₂ molecule (${\mu }_{{{\rm{H}}}_{2}}$ = 2.8) divided by the mean molecular weight per particle (μ_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'' (∼1.4 pc), and they have pixel sizes of 115 × 115.

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) line widths tend to be higher (factor of ∼2) as they are 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 5 × 115), and the PA of the rectangle is given by the instantaneous slope of the 24th ordered polynomial fit the spine, as discussed in Section 3. The sliding rectangular box has a width, w, of 20 HAWC+ pixels (74'') and a height, h, of 15 HAWC+ pixels (555), equivalent to a width and height of ∼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 the 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 }$, δ v_los, and δ θ for each rectangular box. In the column density cutout, we take the median value as our measure of the mean column density, ²⁰ ${\bar{N}}_{\mathrm{gas}}$, and subsequently convert it to a number density ${\bar{n}}_{\mathrm{gas}}$ and then $\bar{\rho }$ assuming a cylindrical bone (see Appendix A.2). δ v_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, δ θ_obs, and for its error cutout, we take the median value, which we call σ_θ . The estimated intrinsic angle dispersion, δ θ, can be corrected for observational errors such that $\delta \theta =\sqrt{\delta {\theta }_{\mathrm{obs}}^{2}-{\sigma }_{\theta }^{2}}$. From this we can estimate the plane-of-sky B-field strengths, B_pos. We do not calculate the B-fields at locations with δ θ > 25° (Ostriker et al. 2001) since the turbulence driving the angular dispersion would be super-Alfvénic.

The B_pos map is shown in the left panel of Figure 4, with fields ranging from ∼20 to 160 μG. We then solve for the critical ratio, λ, by taking the ratio of the observed and critical mass to magnetic flux ratios, i.e.,

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{\left(M/{\rm{\Phi }}\right)}_{\mathrm{observed}}}{{\left(M/{\rm{\Phi }}\right)}_{\mathrm{crit}}}\end{eqnarray} \tag{ 4 }$$

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

Zoom In Zoom Out Reset image size

The observed mass within the w × h box is M_observed = μ_p ${m}_{p}{\bar{N}}_{\mathrm{gas}}{wh}$. The magnetic flux is calculated via Appendix A.3, and we approximate ${\left(M/{\rm{\Phi }}\right)}_{\mathrm{crit}}$ as $1/(2\pi \sqrt{G})$ (McKee & Ostriker 2007). The λ 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_pos and λ 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, λ tends to be close to or larger than one, indicating that the areas are typically supercritical to collapse. However, there are several areas along the bone where λ < 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_pos and λ as a function of the average number density, ${\bar{n}}_{\mathrm{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}}_{\mathrm{gas}}$ ≲1700 cm⁻³ (${\bar{N}}_{\mathrm{gas}}\lesssim 1.6\times {10}^{22}$ cm⁻²), λ is typically less than 1 (subcritical), while for higher densities (locations of most YSOs), λ is typically higher and often supercritical. Overall, there is little change in B_pos as a function of ${\bar{n}}_{\mathrm{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 ±0.006 μG/cm⁻³. However, given the dispersion of points over a small range of densities, we cannot draw conclusions from this relation.

Zoom In Zoom Out Reset image size

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 λ is equal to or less than one, this indicates that either we are underestimating λ or that the high values of λ are more localized to YSOs and would necessitate higher resolution polarimetric observations for proper measurement of the higher λ values.

We note that for the B_pos vs. ${\bar{n}}_{\mathrm{gas}}$ panel, there are 3 points with higher field strengths (>100 μ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 appear to be potentially two velocity components that cause the velocity dispersion (and thus the B-field strength) to be overestimated by a factor of ∼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 nonhomologous 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, ²¹ 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.

Zoom In Zoom Out Reset image size

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 δ θ = 11° ²² , which is ∼75% of the value of δ θ 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 δ v_los = 0.8 km s⁻¹. From our model, we find an average number density of 4300 cm⁻³. ²³ From these values, we calculate a mean field in the plane of sky in the unmasked area of B_pos = 56 μG. The peak total field strength for the SOFIA beam is then B₀ = 108 μG, assuming ${\bar{B}}_{\mathrm{tot}}=(4/\pi ){B}_{\mathrm{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_⊙, resulting in a mass-to-flux parameter of λ = 1.7. The sliding box analysis along the spine of this area found comparable B-field strengths in this region of ∼30–75 μG with values of λ between 0.8 and 1.4. We note the sliding box is ∼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).

In Section 4.1, a rectangular box was moved across the spine of the bone, and we only considered pixels within this box for the DCFST technique. We want to determine whether or not a pixel at location x_c , z_c is within a particular rectangle. Consider a rectangle centered at location x_c0, z_c0 with height h, width w, and angle θ which is measured counterclockwise from up (north), as shown in Figure 7. We define one set of axes with respect to the rectangle, where x is in the direction of the width and z is in the direction of the height. We define a second set of axes in the coordinate system of the map (Galactic north–south, west–east for our case), with the axes centered on the pixel of interest relative to the center of the rectangle, i.e., x_c –x_c0 and z_c –z_c0. The transformation between the coordinate system is then

$$\begin{eqnarray}&&x=({z}_{c}-{z}_{c0})\sin \,\theta +({{\rm{x}}}_{c}-{{\rm{x}}}_{c0})\cos \,\theta \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&z=({z}_{c}-{z}_{c0})\cos \,\theta +({{\rm{x}}}_{c}-{{\rm{x}}}_{c0})\sin \,\theta \end{eqnarray} \tag{ A2 }$$

Zoom In Zoom Out Reset image size

Pixel x_c , z_c will be inside the rectangle if ∣x∣ ≤ w/2 and ∣z∣ ≤ h/2. All pixels meeting these criteria in a given SOFIA map are considered within a particular rectangular box. For our analysis, we used w = 20 pixels = 74'' and h = 15 pixels =555. In our particular case, the chosen pixel size oversamples the beam. However, oversampling does not change the true dispersion, mean, or median, and potentially can estimate these parameters more accurately since at least Nyquist sampling is needed to capture all features of a map.

We want to calculate the volume density for the sliding box to apply the DCFST method. Consider a cylinder (approximation for a bone or filament) with the x-axis along the length of the cylinder, the z-axis perpendicular to the x-axis and in the plane of sky, the y-axis along the line of sight, and a center at (x, y, z) = (0, 0, 0). This box has a height h (Figure 7), which we define to extend from − z₁ to z₁ so that z₁ = h/2. If the observed mean column density within the box, $\bar{N}({z}_{1},R)$, is calculated, then the mean number density within the box, $\bar{n}({z}_{1},R)$, is

$$\begin{eqnarray}&&\bar{n}({z}_{1},R)=\displaystyle \frac{\bar{N}({z}_{1},R)}{2\bar{Y}({z}_{1},R)}\end{eqnarray} \tag{ A3 }$$

where $2\bar{Y}({z}_{1},R)=(1/{z}_{1}){\int }_{0}^{{z}_{1}}{dz}\sqrt{{R}^{2}-{z}^{2}}$ is the path length averaged over the heights 0 to z₁. Substituting the path length in the above equation, we arrive at the final equation for $\bar{n}({z}_{1},R)$ of

$$\begin{eqnarray}&&\bar{n}({z}_{1},R)=\displaystyle \frac{\bar{N}({z}_{1},R)}{\sqrt{{R}^{2}-{z}_{1}^{2}}+\tfrac{{R}^{2}}{{z}_{1}}{\tan }^{-1}\tfrac{{z}_{1}}{\sqrt{{R}^{2}-{z}_{1}^{2}}}}.\end{eqnarray} \tag{ A4 }$$

For G47, we take R = 1.6 pc (Zucker et al. 2018) and h = 2z₁ = 15 pixels = 555. To convert to the mean mass volume density, $\bar{\rho }$, $\bar{n}$ should be multiplied by the mean molecular weight times the hydrogen mass, m_H . If $\bar{n}$ is the gas volume density (used in this study), the mean molecular weight per particle, μ_p = 2.37, should be used; if $\bar{n}$ is the H₂ volume density, then the mean molecular weight per H₂ molecule, ${\mu }_{{{\rm{H}}}_{2}}$ = 2.8, should be used (Kauffmann et al. 2008).

We draw the sliding rectangular box we use in Figure 8 from the plane-of-sky perspective (left panel) and a three-dimensional view (right panel). The x, y, and z axes are along the bone's short axis, the line of sight (LOS), and the bone's long axis, respectively. The magnetic flux is simply the sum of the flux in each dimension, i.e.,

$$\begin{eqnarray}&&{\rm{\Phi }}={A}_{x}{B}_{x}+{A}_{y}{B}_{y}+{A}_{z}{B}_{z}\end{eqnarray} \tag{ A5 }$$

where A_i and B_i are the areas and B-fields for each axis direction. We measure the B_pos from the DCFST method, and we have a typical B-field direction, θ_B , which we take to be the median angle in the sliding box analysis. The true B-field direction has an inclination i along the line of sight (B_pos = B for i = 0). For a rectangular box rotated so that the z-axis has a PA of θ, we define Δθ ≡ θ_B −θ. The solutions for the parameterization for each A_i and B_i are indicated in Figure 8. The formula for _Ax accounts for the fact that it is the crosssectional area of a circular cylinder with parallel planar sides. The inclination i is unknown. The true B-field is typically chosen based on a statistical average with inclination such that B_pos/B = π/4 (Crutcher et al. 2004). We choose i to be based on this average so that $i={i}_{\mathrm{avg}}={\cos }^{-1}(\pi /4)\approx 38\buildrel{\circ}\over{.} 2$. For a given bone radius R and w × h rectangular box, we can now calculate Φ from the measurements of B_pos and θ_B in each box.

Zoom In Zoom Out Reset image size

For a given field strength and orientation, the magnetic flux depends on the geometry of the smoothing box since the projected box area in the field direction varies with the relative box dimensions and with the box inclination. As such, the measurements should be placed in the context of the box dimensions and field direction. Nevertheless, for the adopted box dimensions, changing box orientation by 90° changes the derived flux by less than 25%. For a typical box orientation, reducing the adopted box width w by a factor of 2 reduces the derived flux by ∼35%.