Polarized Light from Massive Protoclusters (POLIMAP). I. Dissecting the Role of Magnetic Fields in the Massive Infrared Dark Cloud G28.37+0.07

Magnetic fields may play a crucial role in setting the initial conditions of massive star and star cluster formation. To investigate this, we report SOFIA-HAWC+ 214 μm observations of polarized thermal dust emission and high-resolution GBT-Argus C18O(1-0) observations toward the massive Infrared Dark Cloud (IRDC) G28.37+0.07. Considering the local dispersion of B-field orientations, we produce a map of the B-field strength of the IRDC, which exhibits values between ∼0.03 and 1 mG based on a refined Davis–Chandrasekhar–Fermi method proposed by Skalidis & Tassis. Comparing to a map of inferred density, the IRDC exhibits a B–n relation with a power-law index of 0.51 ± 0.02, which is consistent with a scenario of magnetically regulated anisotropic collapse. Consideration of the mass-to-flux ratio map indicates that magnetic fields are dynamically important in most regions of the IRDC. A virial analysis of a sample of massive, dense cores in the IRDC, including evaluation of magnetic and kinetic internal and surface terms, indicates consistency with virial equilibrium, sub-Alfvénic conditions, and a dominant role for B-fields in regulating collapse. A clear alignment of magnetic field morphology with the direction of the steepest column density gradient is also detected. However, there is no preferred orientation of protostellar outflow directions with the B-field. Overall, these results indicate that magnetic fields play a crucial role in regulating massive star and star cluster formation, and therefore they need to be accounted for in theoretical models of these processes.

1. INTRODUCTION Magnetic (B) fields have been suggested to play a significant role in the birth of massive stars and star clusters (e.g., see reviews by Tan et al. 2014;Pattle et al. cylaw.astro@gmail.com2022).However, many aspects remain uncertain, including the relative importance of B-fields compared to turbulence, converging flows, feedback and self-gravity in setting the initial molecular cloud conditions and fragmentation properties that ultimately lead to massive star and star cluster formation.Depending on the morphology and strength of B-fields, various models of magnetically regulated star formation have been proposed, including: fragmentation of magnetized filamentary clouds (e.g., Li et al. 2017); collisions of magnetised clouds (e.g., Wu et al. 2015); and collapse after compression in feedback bubbles (e.g., Inutsuka et al. 2015).For massive star formation/accretion, two main scenarios are (i) Core Accretion (e.g., the Turbulent Core Accretion (TCA) model of McKee & Tan 2003) and (ii) Competitive Accretion (e.g., Bonnell et al. 2001;Vázquez-Semadeni et al. 2019;Padoan et al. 2020;Grudić et al. 2022).The former is a scaled-up version of the standard model of low-mass star formation (Shu et al. 1987), although with the internal pressure of the massive prestellar core dominated by turbulence and/or magnetic fields rather than thermal pressure.The fiducial TCA model involves near equipartition between turbulence and magnetic fields, resulting in typical field strengths in massive pre-stellar cores of ∼ 1 mG.In Competitive Accretion, stars chaotically gain their mass in the crowded centres of protoclusters, typically via the global collapse of a cluster-forming clump, without passing through the massive prestellar core phase.This mode has generally been seen in simulations that are either unmagnetized or have relatively weak B-fields.
To uncover the connections between larger-scale environments and smaller-scale star formation properties and to test theoretical star formation models, it is crucial to obtain observational constraints on the magnetic field morphologies, strengths, and energy balances of molecular clouds that are progenitors of massive stars and star clusters.Infrared dark clouds (IRDCs), owing to their dense (n H ≥ 10 5 cm −3 ) and cold (T ≤ 20 K) conditions, are thought to be examples of such clouds (e.g., Tan et al. 2014).However, measuring the magnetic field strength and morphology in IRDCs is challenging due to their relatively far distances (d ≳ 3 kpc).Current and future advances in the polarimetric capabilities of various single dish telescopes (e.g., SOFIA-HAWC+, JCMT-POL2, LMT-TolTEC, and IRAM-30m-NIKA2) and the full polarization capabilities of interferometric facilities (e.g., ALMA, SMA, and VLA) are enabling studies of B-fields in IRDCs (e.g., Pillai et al. 2015;Xu et al. 2023a).Such studies typically rely on mapping polarized thermal dust emission to infer the plane of sky component magnetic field morphology and utilizing the Davis (1951) and Chandrasekhar & Fermi (1953) (DCF) method to infer the magnetic field strength.Other techniques to infer magnetic fields include via the Zeeman effect (e.g., Crutcher & Kemball 2019;Ching et al. 2022), employing velocity anisotropy observed in the spectral line cubes (Xu et al. 2023a), or via the Goldreich-Kylafis effect (e.g., Goldreich & Kylafis 1981, 1982;Forbrich et al. 2008).However, these non-DCF methods have so far been challenging to apply to IRDCs.
Observations of the magnetic field properties toward IRDCs have revealed ordered magnetic field morphology and relatively strong B-fields of up to ∼mG strength, which would be strong enough to influence the cloud dynamics, including fragmentation.A correlation between magnetic field morphology and filament density structure indicating perpendicular alignment with filament orientations has been reported (e.g., Liu et al. 2018;Soam et al. 2019).However, these studies have concerned only a handful of IRDCs, e.g., G11.11-0.12(the "Snake"), G34.34+0.24,G35.39-0.33,and the "Brick" (Pillai et al. 2015;Hoq et al. 2017;Liu et al. 2018;Soam et al. 2019;Tang et al. 2019;Vahdanian & Nejad-Asghar 2022;Chen et al. 2022;Bich Ngoc et al. 2023;Xu et al. 2023a).Hence, a systematic study of IRDCs that span a range of masses, sizes, and environments is an important next step further to understanding the links between the magnetic field and IRDC properties and, thus, better constrain the initial conditions of massive star and star cluster formation.
To more systemically map and constrain the magnetic fields properties in IRDCs, we have initiated the Polarized Light from Massive Protoclusters (POLIMAP) (PI: J. C. Tan) survey with an original goal to map the polarized dust continuum emission of a sample of 10 IRDCs (Clouds A to J) from the sample of Butler & Tan (2009, 2012).POLIMAP utilizes SOFIA-HAWC+ observations at 214 µm (Project ID: 09 0104)1 , complemented by joint Green Bank Telescope (GBT)-Argus 13 CO(1-0) and C 18 O(1-0) observations to probe molecular gas kinematics at high angular resolution of ∼ 7 ′′ .Before the shutdown of SOFIA, POLIMAP HAWC+ observations were completed for IRDCs B, C, F, H, I, J, while GBT-Argus observations have been completed for IRDCs B, C, F, G, H, I, J.In this paper we present the first results from POLIMAP for the massive IRDC G28.37+0.07(Cloud C).
Located at a distance of d = 5.0 kpc from the Sun (Simon et al. 2006), G28.37+0.07 is one of the most massive IRDCs known, with a mass estimated to be about 70, 000 M ⊙ (Butler et al. 2014;Lim & Tan 2014) and hosts multiple massive dense cores (Tan et al. 2013; Barnes et al. 2023).Figure 1 shows a Spitzer mid-infrared (MIR) view of the IRDC, along with MIR extinction (MIREX) contours that help illustrate the cloud's global structure (Butler & Tan 2012;Kainu-lainen & Tan 2013).Overall physical conditions of the cloud, including the mass surface density distributions, gas kinematics, excitation conditions, core mass function, high-mass protostellar populations and their outflows, have been explored in a variety of studies (e.g., Butler & Tan 2009, 2012;Kainulainen & Tan 2013;Hernandez & Tan 2015;Cosentino et al. 2018;Kong et al. 2019;Moser et al. 2020;Entekhabi et al. 2022).However, there has only been limited study of magnetic fields in this IRDC (e.g., Liu et al. 2020), especially on global scales.
This paper is structured as follows.In §2 we describe the observations and data reduction procedures used for the SOFIA and GBT observations, as well as various ancillary datasets.In §3 we present our main observational results for the 214 µm polarimetric imaging (tracing Bfield orientation), C 18 O(1-0) emission, derived B-field strength map, and associated B − N and B − n relations.In §4 we discuss the dynamical implications of the magnetic field, including presentation of the massto-flux ratio map and a virial analysis of dense cores in the IRDC.In §5 we explore the connections between magnetic field morphology and gradients in the mass surface density map and orientations of protostellar outflows.We present a discussion and summary of our results in §6.

SOFIA-HAWC+
We observed the POLIMAP IRDC sample with the Stratospheric Observatory for Infrared Astronomy (SOFIA) using the HAWC+ instrument (Harper et al. 2018) observing in band E (214 µm).HAWC+ is a Half-Wave Plate (HWP) and wire-grid polarimeter operating from ∼ 50 to 214 µm.Rotating the HWP allows the four required angles to be probed at θ = (0, 22.5, 45, and 67.5 degrees).Note, HAWC+ is not sensitive to circular polarization and thus only measures the linearly polarized emissions.
The SOFIA-HAWC+ observations of G28.37+0.07 were carried out on 21, 27 & 28 Sept. 2022 (Proposal ID: 09 0104) as part of SOFIA Cycle 9.The full width at half maximum (FWHM) of the beam is 18. ′′ 2, and the surveyed map size is about 12. ′ 0 × 11. ′ 4 (see Fig. 1).The observations were performed using the Nod-Match chop mode with a Lissajous scan pattern to have fast enough measurements such that the variations in the background do not overwhelm the polarization signal.
The raw data were processed by the HAWC+ data reduction pipeline version 3.2.02 .This pipeline includes different data processing steps, such as corrections for dead pixels and the intrinsic polarization of the instrument and telescope.The final data output of the pipeline is 'Level 4', flux-calibrated, and the resulting pixel size of these data products are 4. ′′ 55.

GBT-Argus
We mapped 13 CO(1-0) and C 18 O(1-0) emissions of the IRDC sample with Argus instrument of the GBT between Jan to April 2022.Argus (Sieth et al. 2014) is a 16-element heterodyne 'radio camera' operating in the 85-116 GHz range.Argus has 16 heterodyne pixel receivers mounted in a 4×4 layout and with a beam separation of 30.′′ 4 in both the elevation and cross-elevation directions.
The observations utilized a 'fast-mapping' method in which the heterodyne pixel receivers are scanned across the sky in rows of Galactic longitude, separated by 5. ′′ 58 (i.e., 0.8× the beam size).This might be expected to undersample the plane of the sky in comparison with Nyquist sampling, in which a separation of less than 0.5 × of the beam size is required.However, the effects of sky rotation in relation to the scan direction and the fact that multiple beams observe any given point in the sky allow for a relaxation of this requirement without a significant effect on the data quality.The map coverage toward G28.37+0.07 has an area of 16. ′ 2 × 13. ′ 2 (see Fig. 1).
The data were reduced using the GBTIDL package, with calibration achieved through observations of the Argus 'Vane', essentially a warm load placed in the beam of the receiver.Atmospheric opacity estimates were taken from the forecast values obtained from the CLEO3 application, which are the average of the values from three surrounding local sites.After baseline fitting and removal were performed on the spectra of each Argus feed, the data were then gridded using the GBO tool gbtgridder with a Gaussian beam assumption and a pixel size of 3. ′′ 0. The resulting cube was then subjected to further iterative baseline fitting and removal to equalize background noise levels.Both the final 13 CO(1-0) and C 18 O(1-0) position-position-velocity (PPV) cubes have a pixel scale of 2. ′′ 57 and a spectral resolution of 0.182 km s −1 .The root-mean-square (RMS) noise level is estimated by taking the average value of RMS noise levels measured from spectra obtained within beamsized regions over the velocity range of 50-65 km s −1 , i.e., away from the emission of the IRDC or other clouds.We find a noise level of 0.55 K per channel (of 0.182 km s −1 velocity range).

Ancillary Data
We utilize previously published maps of the mass surface density, Σ, of the IRDC.First, Butler & Tan (2012) constructed a mid-infrared (MIR) extinction (MIREX) map based on 8 µm Spitzer-IRAC imaging data.A nearinfrared (NIR) extinction corrected version of this map was published by Kainulainen & Tan (2013), which improves the accuracy of the large-scale, lower column density base level of the map.These MIREX maps have a relatively high angular resolution of 2 ′′ , set by the Spitzer-IRAC resolution.They also have the advantage of not requiring knowledge of the temperature of the IRDC material, with the main uncertainty arising from the choice of the opacity per unit mass at MIR wavelengths, estimated to be ∼ 30% from comparison of different dust models (Butler & Tan 2012).However, the MIREX maps have a disadvantage in that they cannot measure Σ in MIR-bright regions, e.g., MIR-bright sources, which appear as "holes" in the map.In addition, the MIREX method has a maximum level of Σ that it is able to probe, which for IRDC C is ∼ 0.5 g cm −2 , so in high-Σ regions it will tend to underestimate the true mass surface density.
An alternative method of estimating Σ is from modelling the far-infrared (FIR) emission.Lim et al. (2016) derived a Σ map of IRDC C by fitting a graybody emission model to Herschel-PACS and SPIRE fluxes from 160 to 500 µm.This analysis also yields a map of the dust temperature of the IRDC.Since the IRDC is close to the Galactic plane, there is significant FIR emission from material in the diffuse ISM along the line of sight.Lim et al. (2016) explored various methods for correcting for this background and foreground emission.Our analysis uses the map derived from their Galactic Gaussian (GG) foreground-background subtraction method.The FIR emission derived Σ map has a resolution of 18 ′′ (i.e., the "hi-res" version), while the temperature map has a resolution of 36 ′′ (see Lim et al. 2016, for details).

214 µm Polarimetric Imaging
The 214 µm intensity (Stokes I) map of the IRDC is shown in Figure 2.This map also shows 90-degreerotated polarization vectors, which are expected to trace the material's emission-weighted plane-of-sky magnetic field orientations along the line of sight.The spacing of these vectors is set to one SOFIA-HAWC+ beam at 214 µm, i.e., 18. ′′ 2. The inferred magnetic field position angle orientations follow the IAU convention, measured counterclockwise from Galactic north.In Figure 3 we also show the expected B-field orientation in the form of a "drapery" pattern overlaid on the Spitzer IRAC 8 µm map.Appendix A presents maps of Stokes Q, U , and the debiased polarization fraction (P ′ ), and signal-to-noise ratio (SNR) of Stokes I. Since SOFIA-HAWC+ does not measure circular polarization, the Stokes parameter V is assumed to be 0.
The Stokes I map shows strong, concentrated emissions at several locations.In particular, the brightest region in Figure 1 corresponds to a MIR-bright region at the eastern side of the IRDC, with I 214µm ≃ 10 4 MJy sr −1 .The main spine of the IRDC extends SW and then W from this location and is seen as an emission feature at 214 µm with intensity values of ∼ 10 3 MJy sr −1 .Two other FIR-bright regions with I 214µm ≃ 10 4 MJy sr −1 are present to the north of the main IRDC spine, with corresponding MIR-bright sources visible in Fig. 1.
The typical uncertainties of Stokes I in the main part of the map are at a near constant level of ∼ 9 MJy sr −1 , but with local peaks of ∼ 16MJy sr −1 at the locations of the brightest FIR regions.Thus the SNR map of Stokes I (see Appendix A) shows typical values of ∼ 150 in the main IRDC spine, rising to ∼ 500 in the FIR-bright sources.
In this paper, we focus on analysis of the orientation of the polarization vectors.We defer analysis of the degree of polarization, i.e., polarization fraction, and the implications for dust grain alignment mechanisms and size distributions for a future paper in this series.Nevertheless, we note that the map of P ′ (see Appendix A) shows typical values from a few to ∼ 10% in the main spine of the IRDC, but with the values dropping to ∼ 0.5% toward the brightest MIR region.Furthermore, in Table 1 we list the values of the derived Stokes parameters I, Q, U , magnetic field position angles (θ), and degree of polarization (P and P ′ ) for selected regions in the IRDC, i.e., the massive protostar (Cp23 Moser et al. 2020), which corresponds to the brightest FIR source, and 16 massive starless/early-stage core/clumps from the sample of Butler et al. (2014) (the locations of these sources are shown in Fig. 3).For each of these sources, the quantities have been evaluated within a circular aperture of radius of 9 ′′ , i.e., at the resolution of the 214 µm HAWC+ image.The polarization percentages shown in Table 1 follow what was noted above, i.e., the bright protostellar source has P ′ ≃ 0.5%, while the starless/early-stage cores/clumps generally have larger values, typically from a few to ∼ 10%.We return to the discussion of the dynamics of these sources in §4.2.Note- † Moser et al. (2020).* If the error in the polarization fraction is larger than the polarization fraction, then the debiasing fail, and the debiased polarization will be flagged as 0%.
When considering the derived orientations of the polarization vectors, it is important to note the uncertainties in θ. Figure 4 presents a map of these uncertainties.These are also presented in Table 1 for the example sources, which range from a few to ∼ 30 • .Such values are also typical in the main spine of the IRDC, but larger uncertainties are found in regions away from the cloud that are relatively faint in their FIR emission.

GBT C 18 O(1-0) spectrum and moment maps
In Figure 5, we present the average spectrum of C 18 O(1-0) emission observed by the GBT extracted from the IRDC ellipse region (see Fig. 1).This spectrum shows a single main emission feature extending from about 70 to 86 km/s and with a central peak near 79 km/s, which is the well-established value for dense gas in G28.37+0.07(e.g., Tan et al. 2013;Hernandez & Tan 2015).A simple Gaussian fit to this feature results in an average velocity centroid of 78.9 km s −1 and a total 1-D velocity dispersion of σ v = 2.58 km s −1 .
In Figure 6 we present the zero, first and second moment maps of C 18 O(1-0) emission of G28.37+0.07,evaluated over the velocity range from 70 to 86 km/s.Peak values of integrated intensity reach values of about 12 K km/s.However, there is significant emission at levels of ∼ 6 K km/s over most of the mapped region.Over most of the field of view, the first moment map shows quite constant values between about 78 and 80 km/s.However, some localized, higher-velocity features are seen in the SW region of the map.The second moment map shows values of just over 3 km/s in the dense regions of the IRDC, with moderately larger values in the surroundings.Finally, Figure 6d presents a map of the uncertainty of the velocity dispersion, estimated utilizing the bettermoments python package (Teague 2019).We see it takes typical values of about 0.5 km/s, i.e., about a 15% uncertainty.

Mapping Magnetic Field Strength
We construct maps of magnetic field strength of G28.37+0.07 using the data presented in the previous sections.We first use the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951;Davis & Greenstein 1951;Chandrasekhar & Fermi 1953) to estimate the plane of sky component of the magnetic field strength.The DCF method estimates magnetic field strength via: where Q is a factor related to the geometry (Heitsch et al. 2001;Padoan et al. 2001;Ostriker et al. 2001, usually set equal to 0.5), ρ is the gas density, σ v is the 1-D velocity dispersion, and σ θ is the dispersion in polarization position angles in radians.
A modified version of the DCF method has been proposed by Skalidis & Tassis (2021) (ST), which we will refer to as "refined" DCF (r-DCF).Here the magnetic field strength is estimated via: In general the r-DCF method is expected to be more accurate in regions sub-to trans-Alfvénic turbulence, where it derives lower magnetic field strengths than the DCF method.We also note that from the above equations, the two methods will give the same estimate of field strength (for the same input physical conditions and with Q = 0.5) when σ θ = 0.5 radians, i.e., 28.6 • .The procedure to construct these magnetic field strength maps is as follows.We first re-grid the original Stokes parameter maps into a coarser pixel scale with a size that corresponds to the SOFIA band-E beam FWHM = 18.2 ′′ so that the measurements in each pixel are spatially independent of each other.We then compute the polarization angle quantities, including the po- larization position angles and the rotated angles (i.e., magnetic field directions) from these re-gridded maps.
We then measure σ θ from these maps via where the summation is done over a set of measured polarization angles θ i and ⟨θ⟩ is the mean polarization orientation within the region of consideration.The expression inside involves the absolute values of angle difference to ensure the angles remain in the range -90 • to +90 • (Pelgrims 2019).
To construct the highest possible resolution B-field strength map we work with a 2 × 2 pixel grid (i.e., 36.′′ 4 × 36.′′ 4) as an input window in which to evaluate σ θ .This window is then translated one pixel at a time for the next evaluation.Thus most pixels in the image, except for those at the edge, receive four measurements of σ θ , with the final reported value being the average of these measurements.We also carry out a version of our analysis with a 3 × 3 pixel grid window, with these results presented in Appendix B.
The map of the mean polarization angle, ⟨θ⟩, is presented in Figure 7a and the map of polarization angle dispersion in Figure 7b.We notice a typical angle dispersion toward the main spine of the IRDC of a few degrees to ∼ 15 • .We also note an increase in the angle dispersion in the midpoint of the main spine, which is where a global change of field orientation from mainly perpendicular to mainly parallel to the filament is seen.A few thin "lines" of higher dispersion of ∼ 25 • are consistent with boundaries between regions that have particular orientations of field directions.
To measure σ v we utilize the C 18 O(1-0) data shown in Figure 6c, but with a small correction to subtract off, in quadrature, the thermal broadening of the line.For this we utilize the dust temperature map of Lim et al. (2016), assuming gas and dust are at similar temperatures.However, it should be noted that this correction is a very minor effect, given the observed values of velocity dispersion of ≳ 3 km/s and the cold temperatures (≲ 15 K) of the IRDC.
To estimate the density of the region we utilize the mass-weighted density map of (Xu et al. 2023b), which is derived from the mass surface density map via machine learning methods, with the algorithm trained on a set of numerical simulations of magnetized GMCs.The density map is shown in Figure 7d.
Combining the above information, we then create two versions of the magnetic field strength map, i.e., one based on the DCF (Eq. 1) and one based on the r-DCF (Eq.2).Logarithmically and linearly scaled versions of these maps are shown in Figure 8.We note that the resolution of the maps is 36.4′′ , i.e., 0.85 pc, set by the 2 x 2 pixel grid size used for the σ θ measurement.As expected, the DCF method returns higher B-field strength estimates than the r-DCF method.In the DCF B-field strengths map we see B-field strengths up to about 3 mG in localized regions, especially at the western end of the main filament of the IRDC.However, we regard the r-DCF map as being more accurate.Here the peak B-field strengths are about 1.5 mG at the western end of the IRDC.Going along the main spine to the east, these values decrease to be ∼ 0.5 mG.
The uncertainties in magnetic field strength are computed given estimated errors for the density, velocity dispersion, and polarization angle dispersion.Uncertainties in velocity dispersion over a 2 x 2 pixel grid are relatively small, i.e., about 5-10%, and so are negligible in the final uncertainty of field strength.For the uncertainty in density, based on the results of Xu et al. (2023b), we adopt ∆n H /n H = 0.5.We estimate the uncertainty in σ θ by sampling the individual polarization orientations within their estimated errors, running our algorithms on these maps, and then examining the dis-persion in the resulting derived values of σ θ .A map of the uncertainty in σ θ , i.e., ∆σ θ , is shown in Figure 9. Maps of the overall uncertainties of B-field strength via the DCF and r-DCF methods are shown in Figure 10.These have typical values of 20% and 15% in the IRDC for the DCF and r-DCF maps, respectively.
We also note that observational uncertainties in polarization position angle tend to inflate the measured dispersion in these angles, thus leading to a systematic bias of underestimation of B-field strengths.These uncertainties then also set a minimum value for σ θ , and thus a maximum value for B-field strength that can be measured in a given region (see also Pattle et al. 2017;Skalidis et al. 2023).The dependencies of magnetic field strength on mass surface density, Σ (or equivalently N H ), and volume den-sity, ρ (or equivalently n H ), including both mean and dispersion, are constraints on theoretical models and simulations of molecular clouds as they undergo collapse to form stars.The magnetic field strength estimated by the DCF and r-DCF methods is the plane of sky component, i.e., B 2D .To make comparison with Zeeman-based relations that measure 1D line of sight magnetic field components (e.g., Crutcher et al. 2004), we estimate the 1D field strength to be B 1D = B 2D / √ 2. Figure 11 presents the DCF and r-DCF B 1D versus N H relations for IRDC G28.37+0.07.Note that here we have measured N H from the Herschel, i.e., FIR, derived Σ map.In addition to the full data set shown by the black dots, we also display the median values of binned data, using bin widths of 0.2 dex in N H unless such a bin would contain fewer than five data points (relevant only for the lowest column density bin).We fit a power law relation of the form where N H,22 ≡ N H /10 22 cm −2 .For the DCF method we derive α B−N = 0.384 ± 0.028 and B N 22 = 0.269 ± 0.151 mG.For the r-DCF method we find α B−N = 0.384±0.017and B N 22 = 0.149±0.089mG.We note that the values of α B−N derived from the DCF and r-DCF methods are nearly identical.
Figure 12 presents the DCF and r-DCF B 1D versus n H relations for IRDC G28.37+0.07.Note that here we have estimated n H as the machine learning inferred value derived from the Herschel, i.e., FIR-based, Σ map.We fit a power law to the binned median values of the form: where n H,4 ≡ n H /10 4 cm −3 .For the DCF method we derive α B−n = 0.511 ± 0.045 and B n4 = 0.352 ± 0.182 mG.
For the r-DCF method we find α B−n = 0.512 ± 0.024 and B n4 = 0.193 ± 0.092 mG.Again, the DCF and r-DCF methods return very similar values of the power law index, with α B−n ≃ 0.5.The indices of the B − N and B − n relations provide constraints on the relative importance of magnetic fields in cloud contraction and dense core collapse (Crutcher et al. 2010;Li et al. 2015;Skalidis et al. 2022).In the paper by Crutcher et al. (2010), the B − n relation was characterized by a piece-wise function in which the magnetic field remains uniform until a threshold density of n H ∼ 300 cm −3 ).At higher densities, the field strength increases with an index α B−n ≃ 2/3.Such an index is expected for the case of isotropic collapse where the magnetic field has only a minor role in the cloud dynamics (Mestel 1966;Crutcher et al. 2010).On the other hand, for magnetic field regulated contraction (e.g., via ambipolar diffusion; Mouschovias 1976a,b) a scaling index shallower than 0.67 is expected, as suggested by multiple studies (e.g., Tritsis et al. 2015;Li et al. 2015;Cao & Li 2023).Thus, our results appear to favor a field regulated contraction and collapse.

DYNAMICAL IMPLICATIONS OF THE MAGNETIC FIELD
We have constructed magnetic field strength maps of IRDC G28.37+0.07 based the SOFIA-HAWC+ polarized dust emission data via the DCF and r-DCF (ST) methods combining the GBT-Argus 1-D velocity dispersion map of C 18 O(1-0) and the number density map derived from a machine learning method that is inferred from a Herschel FIR derived Σ map (Xu et al. 2023b).In this section we characterize the dynamical state of IRDC G28.37+0.07,its dense starless/early-stage core/clump population (Butler et al. 2014) and the massive protostellar core Cp23 residing at the MIR/FIR peak position (Moser et al. 2020).

Mass-to-flux ratio map
We first consider the mass-to-flux ratio (M/Φ) map of IRDC G28.37+0.07 to evaluate the relative importance of the magnetic field to gravity.The equation that computes the corresponding mass-to-flux ratio maps in units of the critical value is (Nakano & Nakamura 1978): where G is the gravitational constant, Σ is the mass surface density and B 1D is the 1-D magnetic field strength (obtained by scaling the plane of sky magnetic field by 1/ √ 2).(M/Φ) crit is the critical mass-to-flux ratio, where defined as If µ Φ > 1, magnetic fields are not able to support the cloud against gravity.Conversely, if µ Φ < 1 then magnetic fields are strong enough to counter gravitational collapse.
To evaluate the mass-to-flux ratio (M/Φ) conditions in the IRDC, we use our derived magnetic field strength maps (i.e., DCF and r-DCF cases, scaled by 1/ √ 2) and the FIR-derived mass surface density map.These maps are shown in Figure 13.As expected, the DCF-based map shows smaller values of µ Φ than the r-DCF-based map.Considering the cloud defined by the ellipse from Simon et al. (2006), we find an average mass-to-flux ratio of ∼ 0.38 ± 0.42 for DCF and 0.54 ± 0.36 for r-DCF.Considering the main spine of the IRDC, we notice that whether or not this structure is magnetically supercritical (µ Φ > 1) or not, depends on whether the DCF or r-DCF method is used.To the north of the main IRDC spine, we see two localized spots and a connecting arc region where there are high, supercritical mass-to-flux ratios.

Energetics
Here, we study the energy budget of the global IRDC and some of its dense substructures to characterise the relative importance of gravity, turbulence, and magnetic fields.
The dynamical state of a spherical cloud is described by the virial theorem.Following Bertoldi & McKee (1992) (hereafter BM92), a cloud in virial equilibrium satisfies the following condition: where E K , E B , and E G are the total kinetic, magnetic, and gravitational energies, respectively.The gravitational energy of an ellipsoidal cloud is given by where G is the gravitational constant, M is the cloud mass, R is the cloud radius, a 1 is a parameter that describes the effect of the internal density distribution of the cloud, and a 2 is a parameter that describes the effect of cloud ellipticity.For a power-law density distribution, ρ ∝ r −kρ , a 1 = (1 − k ρ /3) / (1 − 2k ρ /5).Following the study of Butler & Tan (2012), we adopt k ρ = 1 so that a 1 = 10/9.We will only consider the effects of ellipticity for the case of the global IRDC, which has a ratio of major to minor axes of 1.25.For this case, the value of a 2 = 1.05 (BM92).Following BM92, the kinetic energy is defined by where σ v is the total 1-D velocity dispersion, V is the volume of the cloud, and subscript "0" refers to ambient conditions around the cloud.Finally, the magnetic energy is where B 3D = √ 3B 1D is the total magnetic field in the cloud.
With our mass surface density, velocity dispersion and magnetic field strength maps we are able to estimate E K , E B , and E G , and thus test equation 8.For a cloud defined by a given aperture on the sky, the mass is most simply estimated from the total mass surface density times the projected area of the aperture.However, we expect that some of the material along the line of sight should be considered to be part of the surrounding ambient cloud, so for the dense core sample within the IRDC we also evaluate "background-subtracted" masses in which the background Σ is evaluated in an annulus from R to 2R around a given aperture.
Our estimates of the IRDC and core masses are listed in Table 2.The core masses span about an order of magnitude, ranging from ∼ 200 M ⊙ to ∼ 2000 M ⊙ .We also notice that the FIR-based masses are generally larger than the MIREX-based masses by about a factor of 1.5.This level of difference is within the systematic uncertainties of the two methods, e.g., involving assumptions about dust opacities.We also notice that in some cases (i.e., C3, C9-MIR, C12, C14) the background subtracted mass becomes very small, i.e., < 20% of the non-background-subtracted mass.Especially in these instances, we consider that the background subtracted mass is not reliably measured and we will set the level of the background as a measure of the uncertainty.Table 2 also reports an estimate of the volume density (in terms of n H ) of the sources, with this based on a simple assumption of spherical geometry of the structure.The density of the background, which is used for the surface terms around the core, is estimated in a similar way.
Table 2 next reports the 1D velocity dispersion of the sources, based on the C 18 O(1-0) emission.The velocity dispersion of the background region is also estimated from the spectrum of its C 18 O(1-0) emission.Note that we do not attempt to estimate a background subtracted velocity dispersion.The final primary measurements reported in Table 2 are the total magnetic field strengths, B 3D , as measured from the r-DCF (ST) map (i.e., scaling B 2D by 3/2) of the source and background.We note that the magnetic field in the annular region is often at a similar level to that in the core.For this reason, we also do not attempt to estimate a background subtracted magnetic field strength for the core.With the above measurements of primary quantities, we then estimate E G , E K , and E B .Note that E K and E B include a subtraction of surface terms and in certain cases where the background is relatively strong, these can then take negative values.To test the condition of virial equilibrium, in Figure 14 we plot 2E K +E B versus |E G |.Here we examine four cases resulting from whether we measure mass and density via the FIR or MIREX Σ maps and whether or not we carry out background subtraction to measure the core properties of mass and density.
From Figure 14 we see that for the FIR-derived core properties and the case before background subtraction, the core population scatters quite closely about the condition of virial equilibrium.The most extreme outlier to the gravity dominated regime is the massive protostellar source Cp23.The most extreme outlier to the 2E K +E B dominated regime is the C1 core/clump that is a prime candidate for hosting massive pre-stellar and earl-stage cores (Tan et al. 2013(Tan et al. , 2016;;Kong et al. 2018).Indeed, here we see that E B dominates the energetics.The effect of using MIREX based masses and densities is a modest shift toward supervirial conditions, expected given the systematically lower masses found by the MIREX method.However, as discussed by Butler & Tan (2012), the MIREX masses may be somewhat underestimated in the densest regions given the effects of "saturation" limiting the maximum Σ that can be probed by the method.
The effect of background subtraction on the mass and density estimates is more significant.As noted above, there are cases where the mass of a core is reduced by a large factor (> 5) after attempting background sub-traction.Such estimates thus become highly uncertain and care should be taken with the interpretation of the results.
To gauge the effects of the surface terms, Figure 15 shows the same quantities as Figure 14, but now with E K and E B evaluated without the influence of the surface terms.We see that the surface terms have a large impact the position of cores in this energy ratio diagram.Thus it is important to attempt to include them, even if their evaluation is somewhat uncertain.
Traditional studies of cloud dynamics have focused mostly on measuring the virial parameter, which describes the ratio between kinetic and gravitational energies (BM92): To place our sources in the context of such previous studies, here we also evaluate their virial parameters, which are shown in Fig. 16 and listed in Table 3.
For the case of no background subtraction in estimating core mass, we find average and dispersion of α vir = 1.11 ± 0.587 for FIR derived masses and 1.65 ± 0.666 for MIREX derived masses.If we consider background subtracted masses, the corresponding average virial parameter values are 9.22 ± 25.8 and 2.89 ± 7.89 for FIR and MIREX traced cores, respectively.
Next, we evaluate the Alfvén Mach number, m A , of the cores, given by m ).These values are also shown in Table 3.For the case of no background subtraction in estimating core mass, we find average and dispersion of m A = 0.438 ± 0.130 for FIR derived masses and 0.344 ± 0.104 for MIREX derived masses.If we consider background subtracted masses, the corresponding values are 0.277 ± 0.100 and 0.218 ± 0.0756 for FIR and MIREX traced cores, respectively.
Similarly, the mass-to-flux ratios are also presented in Table 3.For the case of no background subtraction in estimating core mass, we find average and dispersion of µ = 0.878 ± 0.284 for FIR derived masses and 0.539 ± 0.155 for MIREX derived masses.If we consider background subtracted masses, the corresponding mass to flux ratios are 0.380 ± 0.231 and 0.219 ± 0.0820 for FIR and MIREX traced cores, respectively.
Finally, we compare the observed quantities to the prediction of the Turbulent Core Model (TCM) of (Mc-Kee & Tan 2003).The TCM describes the properties of virial equilibrium pressure-truncated singular polytropic spheres that are partially supported by large-scale B−fields and partially by turbulence.The expected 1D velocity dispersion of such a core in the fiducial case of Alfven Mach number of unity is given by (Tan et al. 2013): (13) In Figure 17, we compare the observed velocity dispersion of the cores to the expected velocity dispersion from a fiducial virialized core of given mass, M c , and background clump mass surface density, Σ cl (see also Table 3).For the FIR-based masses, we find that most observed velocity dispersions are smaller than predicted by about a factor of 0.8.For masses based on the MIREX Σ map, the observed velocity dispersion is closer to that predicted by the TCM.
We also estimate the magnetic field at which a core is virialized as predicted by the TCM, i.e., assuming an Alfvén Mach of unity.We compare in Figure 18 the predicted magnetic field strengths and the observed field strengths.Typically, the observed field strengths are stronger than predicted by factors of up to a few.For the predicted TCM field strengths to be larger would imply sub-Alfvénic conditions in the cores.Tan et al. (2013) used the ratio of the observed velocity dispersion to the predicted velocity dispersion to estimate what magnetic field strength is needed to support a core in virial equilibrium.In Figure 19, we plot the corresponding corrected magnetic field strengths to the observed field strengths and find that with these corrections, the predicted field strengths are now in better agreement with those that are observed in the cores.This indicates that the observed magnetic field strengths are at a level to control support of the cores under sub-Alfvénic conditions.

Magnetic field orientation and mass surface density gradient
Here, we examine potential connections between magnetic field morphology and the formation of dense gas structures in the IRDC.First we examine the distribution of angles, Φ, between direction of plane of sky magnetic field orientation and gradient in the mass surface density.We carry out a histogram of relative orientations (HRO) analysis following the method described in (Soler et al. 2013).We examine cases with the Herschel (FIR)-derived Σ map and the MIREX Σ map, each regridded to the 18. ′′ 2 pixel scale of the HAWC+ image.We compute the gradient by convolving with a Gaussian kernel with a FWHM of 3 pixels (i.e., 55. ′′ 6).The HROs are presented in Figure 20.We set the ranges of Σ so that each HRO set contains 25% of the pixels..325(0.166) † An effective radius of 7.06 pc, i.e., 5.37 ′ , is used for the global IRDC, while a radius of 0.399 pc, i.e., 0.303 ′ , is used for the dense cores.For each source, the top row reports mass surface densities and masses based on the Herschel (FIR) mass surface density map, while the bottom row uses the MIREX map.Values in "()" are background subtracted estimates.When a background-subtracted mass or density is found to be negative, it is not evaluated and the symbol "-" is used.The values of σv and Btot in "[]" in the second rows are background values, i.e., σv,0 and Btot,0.An effective radius of 7.06 pc, i.e., 5.37 ′ , is used for the global IRDC, while a radius of 0.399 pc, i.e., 0.303 ′ , is used for the dense cores.For each source, the table reports: virial parameter (αvir); Alfvén Mach number (mA); mass-to-flux ratio (µ); fiducial TCM predicted virialized 1D velocity dispersion (σvir); fiducial TCM predicted core B−field strength (Bvir); and TCM B−field strength needed to achieve virial equilibrium give observed velocity dispersion (Bvir,corr).The top row for each source reports quantities based on masses and densities measured from the Herschel (FIR) mass surface density map, while the bottom row uses the MIREX map.Values in "()" are background subtracted estimates.When a background-subtracted mass or density is found to be negative, it is not evaluated and the symbol "-" is used.Background subtraction is not evaluated for the global properties of Cloud C. The results of Figure 20 generally show a greater number of pixels having smaller values of Φ, i.e., indicating that magnetic field direction tends to align with gradient in Σ (and thus also volume density).However, the distributions are quite broad, with significant numbers of pixels at large values of Φ.There is no clear trend of the HROs varying systematically with Σ.
We evaluate the significance of any overall parallel/perpendicular alignment by evaluating the projected Rayleigh statistic (PRS) (e.g., Jow et al. 2018;Fissel et al. 2019).The PRS is used to test for the existence of a preference for perpendicular or parallel alignment and has the following equation where N is the number of relative orientations between the two sets of angles.The corresponding uncertainty in the PRS is given by (Jow et al. 2018): We apply the PRS and present the Z x values in Table 4.A value of Z x > 0 indicates a preferred parallel alignment, while Z x < 0 indicates a preferred perpendicular alignment.While the values are all > 0, given the uncertainties, there is no strong evidence for preferred parallel alignment in any of the individual quartiles.Table 4 also reports the value for the PRS of the whole pixel sample.
The PRS can infer if there is a uni-modal distribution but cannot reflect on the existence of any multimodal, non-uniform distributions, e.g., a bimodal distribution.Thus, next we perform a Kuiper test, which is closely related to the Kolmogorov-Smirnov (KS) test.The Kuiper test compares the distribution to the cumulative density function of a certain distribution, which in this case is a uniform distribution with the same number of measurements within 0 to 90 degrees, i.e., no preferred direction.We use the python function kuiper two in astropy to carry out the Kuiper two samples test.The output p-value reflects whether the sample distribution differs from that of a uniform distribution.We also summarize the Kuiper test results in Table 4.We see that certain quartiles, corresponding to the ones with the largest PRS values, have small (≲ 0.05) values of p Kuiper , indicating distributions that are inconsistent with a uniform, i.e., random, distribution.In addition, the result for the entire pixel set shows a significant result that the distribution is inconsistent with that of a uniform distribution.

Magnetic fields and protostellar outflows
Here, we investigate if there is any correlation between magnetic field orientation and protostellar outflow axis orientation.We use the sample of 64 averaged blue-and red-shifted CO(2-1) plane-of-sky outflow axis position angles, θ outflow , presented by Kong et al. (2019).These authors already noted that their sample tends to align perpendicularly to the main IRDC filament.
Figure 21 shows a scatter plot of θ outflow versus θ B .Note, that multiple outflows can be located within a SOFIA beam, i.e., sharing the same 18.′′ 2-size pixel in the B−field map.In this case, we assume that all the outflows share that same field orientation.Indeed, the distribution of outflows is clustered so that 49 of the outflows are found to be launched from within just 7 of the pixels of the B−field map.Examining Figure 21, we notice that there seems to be a tentative trend that the relative orientations between outflows and local magnetic fields are either parallel or perpendicular.However, there are also significant numbers of sources that have no preferential orientations.We also color code the outflows with masses that are measured base on the ALMA continuum image (Kong et al. 2019).
In Figure 22 we present a histogram of the relative orientations between the average outflow orientations and the magnetic field orientations.The apparent bimodality is also apparent in this histogram.However, when we apply the PRS and Kuiper tests to the distribution of relative orientations, the results are Z x = 1.81 ± 2.59 and p = 0.65.Thus the data do not indicate any significant evidence for a preferential alignment direction or a distribution that is distinct from a random, uniform distribution.In this paper we have presented the first global, relatively high-resolution (18 ′′ ) magnetic field strength map of the massive IRDC G28.37+0.07utilizing SOFIA-HAWC+ 214µm polarized dust emission data, combined with a velocity dispersion map obtained with GBT-Argus C 18 O(1-0) observations and mass surface density and associated inferred volume density maps derived from Herschel FIR emission and Spitzer MIR extinction maps.We have assessed the role of magnetic fields on the star formation in of G28.37+0.07 based on an evaluation of the B-Σ and B-ρ relations in the maps and more directly by measuring the energy balances between magnetic fields, turbulence and gravity in the global cloud and a sample of dense cores.We have also examined the morphological connection of the magnetic field orientation with density (i.e., Σ) gradients and protostellar outflow axis orientations.Here, we discuss the implications of these results for a more general understanding of the role of magnetic fields in massive star formation.
The relative importance of magnetic fields compared to gravity may influence the observed B-Σ and B-ρ relations of the cloud.Gravity-dominated isotropic contraction is expected to lead to a B-ρ index of about 2/3, while our results imply a shallower index of about 0.5 ± 0.05.Such a shallow value may be consistent with dynamically important magnetic field regulated collapse via ambipolar diffusion (e.g., Cao & Li 2023).
More directly, the mass to flux ratio maps indicate that magnetic fields are dominant (or at least significant) in large regions of the IRDC.This is especially true for values of B-fields estimated by the DCF method.However, even for the lower B−field values of the r-DCF method of Skalidis & Tassis (2021), we find that major parts of the IRDC, including the eastern end of its dense filamentary spine are magnetically sub-critical (µ < 1).An assessment of the global energetics finds that magnetic energy is much larger than gravitational and turbulent energies.We note here that a caveat of the DCF and r-DCF methods is that they ignore the influence of self-gravity in bending B−field orientations.However, such effects would lead us to systematically underestimate true B−field strengths, so our conclusions about the dynamical importance of the fields are quite conservative.
Zooming into sub-regions of the cloud, while there are some variations in the mass to flux ratio, overall, we see that for most cores, magnetic fields play a dynamically important role compared to turbulence.A detailed study of virial balance indicates that cores may be close to virial equilibrium but with relatively strong B−fields that imply sub-Alfvénic conditions.Such Bfields with ∼ 1 mG strengths are also likely to be im-portant for setting the mass scale of the cores at a level of ∼ 100 M ⊙ and by preventing fragmentation to much smaller masses, such B−field strengths may be an essential requirement for the formation of massive stars (e.g., Butler & Tan 2012;Tan et al. 2013).
From the morphological point of view, we find strong evidence for a preferential alignment of B−field orientation with density gradient across the full range of Σ.This is consistent with previous results that found dense filaments have orientations perpendicular to B−field orientation (e.g., Soler et al. 2013).It appears that most of the range of Σ that is probed by our maps is in this regime.Again, the implication is that magnetic fields are helping to shape the collapse of gas into structures that then tend to align perpendicularly to the field.
Once protostellar sources form and drive outflows, one may also ask if the direction of the outflows is influenced by the large scale B−field.Here we do not find strong evidence for such an effect, even though a preferred orientation of the outflow axes has been claimed (Kong et al. 2019) and such effects have been noted for mainly low-mass protostellar outflows (Xu et al. 2022).
A summary of our main conclusions is as follows: • We derived the B−field strengths maps for G28.37+0.07 with both the DCF and r-DCF (Skalidis & Tassis 2021) methods, with the latter considered to be most accurate, given the subto trans-Alfvénic conditions in the IRDC.The B−field strengths in the r-DCF map range from about 0.03 to 1 mG.
• From the r-DCF map we derive a B-N H relation of the form B 1D = 0.15N 0.38 H,22 mG and a B-n H relation of the form B 1D = 0.19n 0.51 H,4 mG.The latter index is relatively shallow compared to some previous estimates, e.g., the classic relation of Crutcher (2012) with a power law index of 0.65.This difference may reflect the particular environment of IRDC G28.37+0.07 or be the result of different methodologies in the estimation of B-fields and density.
• Magnetic fields appear to be playing a dynamically important role on both the global and local, dense core scales of the IRDC.Conditions in the dense cores are consistent with virial equilibrium, although there are significant uncertainties, especially associated with estimation of mass and the effect of clump envelope background subtraction.The most massive protostellar source in the IRDC (Cp23) is seen to be the most gravitydominated source, while the highly deuterated massive pre-stellar core candidate (C1) is the most magnetically-dominated source.
• A morphological study of magnetic field orientation with gradients in the mass surface density map finds strong evidence for a preferentially alignment of the B−fields with the density gradients, again indicating the fields play a role in regulating the collapse.However, the is no evidence that the field is able to control directions of protostellar outflows.
The overall findings of the study demonstrate that magnetic fields play an important role in regulating the formation processes of massive stars by setting up the initial conditions in the massive IRDC G28.37+0.07.Such strong magnetic field may leave traces in anisotropic gas kinematics.A follow-up work towards this IRDC (Law et al., in prep.)  ) to the magnetic field strength assuming the core is at equilibrium with correction of higher magnetic field support (Tan et al. 2013) based on mass derived from the FIR mass surface density map.(b) Top Right: Similar to (a), but the predicted virial magnetic field is based on the background subtracted mass.(c) Bottom Left: Scatter plot of the observed total magnetic field strength (Btot = √ 3Bx) to the magnetic field strength assuming the core is at equilibrium with correction of higher magnetic field support but based on mass derived from the MIREX map .(d) Bottom Right: Similar to (c), but here the predicted virial magnetic field is based on the background subtracted mass

B. EFFECTS ON NUMBER OF PIXELS IN MEASURING DISPERSION IN POLARIZATION ANGLES
We have presented maps of angle dispersion and the corresponding magnetic field strength maps.We select the smallest window size of 4 pixels (equivalent to two SOFIA-HAWC+ beams).However, the low number statistics in the number of vectors in computing the field strengths via (refined) DCF methods may underestimate the angle dispersion, thus overestimating the field strength.We present the mean polarisation angle and the dispersion in the polarization angle based on a 3x3 window centred at each pixel.This analysis shows similar results as the map using the 2×2 window, but with a slightly higher level of angle dispersion, given it is averaging over a larger region.We notice that the magnetic field strength is thus on average slightly lower.Thus, we conclude that there is no significant impact on our main findings by the use of a 2 × 2 window function compared to a 3 × 3 window function.

Figure 2 .)Figure 3 .
Figure 2. The SOFIA-HAWC+ Stokes I intensity map of G28.37+0.07 at 214 µm.Inferred magnetic field position angles are shown with black line segments, separated from each other by a beam size of 18 ′′ .2,i.e., 0.441 pc.The size of the beam is shown by the black circle in the lower left corner.The contours show constant levels of extinction (same as Fig. 1).

Figure 4 .Figure 5 .
Figure 4. Map of polarization angle uncertainties toward G28.37+0.07(note scale is truncated at 45 • ).Generally, the typical uncertainties in the main spine of the IRDC range from a few to 30 • .

)
Figure 6.(a) Top left: Integrated intensity map of C 18 O(1-0) emission of IRDC G28.37+0.07,evaluated from 70 to 86 km s −1 .(b) Top right: As (a), but now showing the first moment map.(c) Bottom left: As (a), but now showing the second moment map.(d) Bottom right: As (a), but now showing the uncertainty in the second moment map.
3.4.The B − N and B − n relations

Figure 7 .
Figure 7. (a) Top Left: Map of mean magnetic field position angle, ⟨θ⟩.(b) Top Right: Map of dispersion in magnetic field position angle, σ θ .(c) Bottom Left: Map of 1-D non-thermal velocity dispersion, σv.(d) Bottom Right: Map of mass-weighted H nuclei number density, nH.

Figure 8 .
Figure 8. Plane of sky magnetic field strength maps of G28.37+0.07 using the DCF and r-DCF (or ST) methods (Left column: DCF; right column: r-DCF).The top row shows the field strength maps with a logarithmic scale, while the bottom row shows the same maps but in linear scale.

Figure 10 .
Figure 10.(a) Left: Map of fractional uncertainty in plane of sky magnetic field strength as estimated via the DCF method.(b) Right: As (a), but now for the r-DCF method.

Figure 11 .Figure 12 .
Figure 11.(a) Left: B1D versus NH relation for IRDC G28.37+0.07,with magnetic field estimated via the DCF method and NH derived from the Herschel FIR Σ map.Black points show data from all pixels in the maps, while red points show median values of binned data, with error bars showing the dispersion about the average.A power law fit to the data (see text) is shown by the solid line.(b) Right: As (a), but now for the r-DCF method for magnetic field strength.

Figure 13 .
Figure 13.(a): Mass surface density in pixel scale equivalent to SOFIA beam FWHM.(b): Dense cores from Butler et al. (2014) overlaid on the magnetic field strength G28.37+0.07measured via the ST method.Each circle is centered at the coordinate of each dense core and with a fix aperture radius equivalent to 1 SOFIA-HAWC+ beam FWHM.The ellipse shows the IRDC region defined by Simon et al. (2006) (c): Mass to flux ratio (µ) map toward G28.37+0.07 based on magnetic field strength estimated from classical DCF method.(d): As (c), but now showing mass to flux ratio (µ) based on magnetic field strength estimated from refined DCF method (Skalidis & Tassis 2021).

Figure 14 .
Figure 14.(a) Top Left: Test of virial equilibrium of IRDC dense cores.Here masses are estimated from the FIR Σ map and no background subtraction is done when evaluating source mass, velocity dispersion and B-field strength.(b) Top Right: As (a), but now with masses estimated from the MIREX Σ map.(c) Bottom Left: As (a), but now with background subtraction when evaluating source mass.(d) Bottom Right: As (c), but now with masses estimated from the MIREX Σ map.

Figure 15 .
Figure 15.(a) Top Left: As Fig. 14a, but now surface terms are not considered when evaluating EK and EB.(b) Top Right: As (a), but now with masses estimated from the MIREX Σ map.(c) Bottom Left: As (a), but now with background subtraction when evaluating source mass.(d) Bottom Right: As (c), but now with masses estimated from the MIREX Σ map.

Figure 16 .
Figure 16.(a) Top Left: Test of virial equilibrium of IRDC dense cores.Here masses are estimated from the FIR Σ map and no background subtraction is done when evaluating source mass, velocity dispersion and B-field strength.(b) Top Right: As (a), but now with measurements from the MIREX Σ map.(c) Bottom Left: As (a), but now with background subtraction when evaluating source mass.(d) Bottom Right: As (b), but now with measurements from the MIREX Σ map.

Figure 17 .
Figure 17.(a) Top Left: Observed total 1-D velocity dispersion against predicted 1-D velocity dispersion from the TCM prediction based on the mass derived from the FIR mass surface density map.(b) Top Right: Similar to (a), but the predicted velocity dispersions are derived from the background-subtracted mass.(c) Bottom Left:Observed velocity dispersion against predicted 1-D velocity dispersions from the TCM with masses are derived from the MIREX map.(d) Bottom Right: Similar to (c), but the predicted velocity dispersions are derived from the background-subtracted mass.

BFigure 18 .Figure 19 .
Figure18.(a) Top Left: Scatter plot of the observed total magnetic field strength (Btot = √ 3Bx) to the magnetic field strength assuming the core is at equilibrium based on mass derived from the FIR mass surface density map .(b) Top Right: Similar to (a), but here the predicted virial magnetic field is based on the background subtracted mass (c) Bottom Left: Scatter plot of the observed total magnetic field strength (Btot = √ 3Bx) to the magnetic field strength assuming the core is at equilibrium but based on mass derived from the MIREX map .(d) Bottom Right: Similar to (c), but here the predicted virial magnetic field is based on the background subtracted mass

Figure 21 .
Figure 21.Scatter plot of outflow axis position angle (Kong et al. 2019) versus local magnetic field orientation.Protostellar core masses are evaluated from ALMA mm continuum data (see text).

Figure 22 .
Figure 22.Histogram of angle difference between protostellar outflow and local magnetic field orientations.

Figure B1 .
Figure B1.(a) Top Left: Mean polarization angle map constructed based on a 3-by-3 window function.(b) Top Right: The corresponding dispersion map of the magnetic field position angles.(c) Bottom Left: The corresponding magnetic field strength map based on the DCF method.(d) Bottom Left: The corresponding magnetic field strength map based on the r-DCF method.

Table 4 .
Results of PRS and Kuiper tests toward the HROs constructed from the FIR and MIREX Σ maps.