Dark Dragon Breaks Magnetic Chain: Dynamical Substructures of IRDC G28.34 Form in Supported Environments

We have comprehensively studied the multiscale physical properties of the massive infrared dark cloud G28.34 (the Dragon cloud) with dust polarization and molecular line data from Planck, FCRAO-14 m, James Clerk Maxwell Telescope, and Atacama Large Millimeter/submillimeter Array. We find that the averaged magnetic fields of clumps tend to be either parallel with or perpendicular to the cloud-scale magnetic fields, while the cores in clump MM4 tend to have magnetic fields aligned with the clump fields. Implementing the relative orientation analysis (for magnetic fields, column density gradients, and local gravity), velocity gradient technique, and modified Davis–Chandrasekhar–Fermi analysis, we find that G28.34 is located in a trans-to-sub-Alfvénic environment; the magnetic field is effectively resisting gravitational collapse in large-scale diffuse gas, but is distorted by gravity within the cloud and affected by star formation activities in high-density regions, and the normalized mass-to-flux ratio tends to increase with increasing density and decreasing radius. Considering the thermal, magnetic, and turbulent supports, we find that the environmental gas of G28.34 is in a supervirial (supported) state, the infrared dark clumps may be in a near-equilibrium state, and core MM4-core4 is in a subvirial (gravity-dominant) state. In summary, we suggest that magnetic fields dominate gravity and turbulence in the cloud environment at large scales, resulting in relatively slow cloud formation and evolution processes. Within the cloud, gravity could overwhelm both magnetic fields and turbulence, allowing local dynamical star formation to happen.


INTRODUCTION
Magnetic fields and turbulence are two major forces resisting the gravitational collapse of molecular clouds in star formation regions (McKee & Ostriker 2007).Observational studies of the magnetic field are important to understand how it regulates star formation and how it is affected by star formation.Observations of polarized dust emission (produced by dust grain alignment, Lazarian 2007;Lazarian & Hoang 2007;Andersson et al. 2015) have been the most widely used technique to trace the plane-ofsky (POS) magnetic field orientation in star-forming molecular clouds (e.g., Hildebrand et al. 1984).There has been an increasing number of both single-dish and interferometric dust polarization observations in molecular clouds1 (Pattle & Fissel 2019;Hull & Zhang 2019).From the observational magnetic field studies, the recent review papers summarised that magnetically transto-super-critical clumps/cores form in sub-critical and trans-to-sub-Alfvénic clouds (Liu et al. 2022a,b).The substructures of clouds may transit to an averagely trans-to-super-Alfvénic state as density increases (Li 2021;Liu et al. 2022b), but this result is less clear due to the uncertainties of the analysis methods.Despite the progress, most of the previous observational studies of magnetic fields were in more evolved star formation regions where significant star-forming activities have taken place (e.g., Sanhueza et al. 2021;Liu et al. 2023a).The general magnetic field properties of clouds at early star formation stages remain under explored.
The massive star formation process is relatively less understood than low-mass star formation partly due to a lack of observations at early evolutionary stages.Massive infrared dark clouds (IRDCs) are believed to harbor the early phase of massive star formation (Pillai et al. 2006;Rathborne et al. 2006;Sanhueza et al. 2012Sanhueza et al. , 2019)), but their weak polarized dust emission makes it more challenging to detect than more evolved regions.So far, there have been only a handful of single-dish (Pillai et al. 2015;Juvela et al. 2018;Liu et al. 2018;Tang et al. 2019;Soam et al. 2019;Añez-López et al. 2020;Ngoc et al. 2023) and interferometric (Beuther et al. 2018;Cortes et al. 2019;Liu et al. 2020) studies of magnetic fields in IRDCs, and there is a lack of multi-scale studies of magnetic fields in the same IRDC.Since the dynamic role of the magnetic field may vary from large to small scales (Liu et al. 2022a,b), it is essential to comprehensively investigate the multi-scale magnetic field in a single IRDC to advance our understanding of the magnetic field properties in the early stages of massive star formation.
G28.34+0.06(hereafter G28.34, also known as the Dragon cloud) is a well-studied massive filamentary IRDC located at a distance of 4.8 kpc (e.g., Carey et al. 1998;Pillai et al. 2006;Rathborne et al. 2006;Wang et al. 2008;Zhang et al. 2009;Lin et al. 2017;Wang 2018).The majority of G28.34 is 8 µm-dark except for the northern end (Zhang et al. 2009).Rathborne et al. (2006) has identified 18 millimeter dust continuum clumps (MM1-18) within G28.34 and in its vicinity.The clumps in G28.34 are found to further fragment into smaller substructures with higher resolution observations (e.g., Zhang et al. 2009Zhang et al. , 2015;;Kong 2019).All the clumps in G28.34 show signs of star formation activities (e.g., Wang et al. 2006;Kong et al. 2019).The large mass reservoir and infrared dark behavior make G28.34 a perfect place to study the extremely early evolutionary stage of massive star formation.Liu et al. (2020) presented a study of the small-scale magnetic field in three clumps (MM1, MM4, and MM9) in G28.34 with Atacama Large Millimeter/submillimeter Array (ALMA) polarization observations.They found that even considering both the turbulent and magnetic support, core MM4-core4 is still in a non-equilibrium state dominated by gravity.As a follow-up work of Liu et al. (2020), in this paper, we utilize the James Clerk Maxwell Telescope (JCMT) dust polarization observations and the Planck dust polarization data to study the multi-scale magnetic field in G28.34 and to determine the multi-scale energy balance in this IRDC.
The data were calibrated with Common Astronomy Software Applications (CASA, McMullin et al. 2007).We performed two rounds of phase-only self-calibrations for the dust continuum.We imaged the molecular line cubes and Stokes I, Q, and U maps of dust continuum using the CASA task TCLEAN with a Briggs weighting parameter of robust = 0.5.We adopted a pixel size of 0.1 ′′ for the imaging.The synthesized beam of the three configurations-combined images is ∼ 0.7 ′′ × 0.5 ′′ (∼0.016-0.012pc at a distance of 4.8 kpc), which improves the resolution of our previous two configurations-combined images (∼ 0.9 ′′ × 0.7 ′′ , Liu et al. 2020).The maximum recoverable scale2 is ∼ 13 ′′ (∼0.3 pc at 4.8 kpc).Before primary beam correction, the 1σ root-mean-square (RMS) noise is σ I ∼0.03 and 0.03 mJy beam −1 for the Stokes I dust continuum maps and σ QU ∼0.01 and 0.012 mJy beam −1 for the Stokes Q or U dust continuum maps of MM4 and MM9, respectively.The debiased polarized intensity P I and its corresponding uncertainty σ P I are calculated as P I = Q 2 + U 2 − σ 2 QU (Vaillancourt 2006) and σ P I ∼ √ 2σ QU , where σ QU is the 1σ rms noise on the background region of the Q or U maps.The polarization position angle θ p is estimated with θ p = 0.5 arctan(U/Q).The uncertainty on the polarization position angle (Naghizadeh-Khouei & Clarke 1993) is given by δθ = 0.5 σ 2 QU /(Q 2 + U 2 ) ∼ 20 • .26(σP I /P I) ∼ 28 • .65(σQU /P I).We only adopt the N 2 D + (3-2) line data in this study.The RMS noises of the N 2 D + line cubes (before primary beam correction) are ∼1.6 and 1.9 mJy beam −1 at a 0.16 km s −1 channel for MM4 and MM9, respectively.All the ALMA images shown in this paper are before primary beam correction.All the continuum fluxes used for calculations of physical parameters are after primary beam correction.

JCMT dust polarization observations and molecular line data
The 850 µm polarized emission of G28.34 was observed by SCUBA-2/POL-2 (Holland et al. 2013;Friberg et al. 2016) on the JCMT between 2022 February 24 and 2022 June 25 under the project M22AP018 (PI: Junhao Liu).The observations were made with the POL-2 DAISY mode with low noise levels in a central region of a 3 ′ radius and with increased noises toward the edge.The spatial resolution of JCMT is ∼ 14 ′′ (∼0.33 pc) at 850 µm.The center of our DAISY field is in the southern part of G28.34 near MM4.A similar POL-2 observation centered at MM1 has been conducted by the JCMT large program B-Fields in STar-Forming Region Observations (BISTRO).A detailed analysis of the multi-scale magnetic field in the infrared bright clump MM1 will be presented by Hwang, J. et al. in prep.as part of the BISTRO survey.
Additionally, we collect the archival 13 CO (3-2) (program: M10AC06) and HCO + (4-3) (programs: M16BP081 and M17BP087) line data observed with the Heterodyne Array Receiver Program (HARP) and Auto-Correlation Spectrometer and Imaging System (ACSIS, Buckle et al. 2009) on the JCMT.The spatial and spectral resolutions are ∼ 14 ′′ and 0.055 km s −1 , respectively, for both lines.To increase the signal-to-noise ratio (S/N), we smooth the two lines to a spectral resolution of 0.2 km s −1 .We estimate the main beam temperature (T mb ) from the corrected antenna temperature (T * A ) adopting a main beam efficiency of 0.64 (Buckle et al. 2009).The typical observation uncertainties of 13 CO (3-2) and HCO + (4-3) are 0.5-1 K and 0.2-0.4K (in T mb ), respectively, per 0.2 km s −1 channel within our studied region.

Planck 353 GHz dust polarization data
We adopt the Planck High Frequency Instrument (HFI, Lamarre et al. 2010) 353 GHz Stokes I, Q, and U maps of the thermal dust emission (version R3.00, Planck Collaboration et al. 2020) toward G28.34 and its surrounding area constructed with the Generalized Needlet Internal Linear Combination method (GNILC, Remazeilles et al. 2011).We also adopt the earlier released dust optical depth (τ 353 ) and temperature maps (version R1.02, Planck Collaboration et al. 2014).The resolution of the Planck maps is 5 ′ (∼7 pc) at 353 GHz.The pixel size of the Planck maps is 1.71 ′ .Within our considered map area (1 • × 1 • ), the mean values for the uncertainties of Q and U (i.e., δQ and δU ) are ∼3 and 4 µK CMB , respectively.The debiased polarized intensity P I and its corresponding uncertainty δP I are calculated as , respectively.We compute polarization position angles in the equatorial coordinates with θ p = 0.5 arctan(U/Q) − ∆θ g−e p , where ∆θ g−e p = arctan( cos(l − 32.9 • ) cos b cot 62.9 • − sin b sin(l − 32.9 is the angle between the galactic and equatorial reference directions (Corradi et al. 1998).For G28.34 at l = 28.34• and b = 0.06 • , we adopt ∆θ g−e p ≈ 62.82 • .The uncertainty on the Planck polarization position angle is given by δθ ∼ 28 • .65(δPI/P I).The Planck and JCMT observations are at the same frequency, so we could compare their intensities to investigate the consistency of the two datasets.To compare the Planck and JCMT data, we convolve the JCMT I and P I maps to the same resolution as the Planck maps.At 5 ′ resolution, the JCMT peak intensity is 98 and 3.1 mJy beam −1 for I and P I, respectively.At the same position, the Planck intensity is 240 and 3.7 mJy beam −1 for I and P I, respectively.The JCMT peak I and P I values are ∼41% and ∼84% of the Planck values, respectively.This comparison suggests that the JCMT POL-2 data filters out a significant amount of the total intensity3 but recovers the majority of the polarized intensity.

Dust polarization and magnetic fields
Here we briefly overview the multi-scale magnetic field structures in G28.34 revealed by Planck, JCMT, and ALMA.Assuming that the observed linear polarization of the dust continuum is due to dust grain alignment, we rotate the dust polarization position angle by 90 • to infer the magnetic field orientation.
Figure 1(a) shows the large-scale magnetic field orientation surrounding G28.34 revealed by Planck.The well-ordered largescale magnetic field shows a predominant northeast-southwest orientation along the galactic plane.
Figure 1(b) shows the magnetic field orientation in G28.34 and several nearby massive clumps revealed by JCMT.The magnetic field morphology in G28.34 is complex.Along the ridge of the main straight dark filament (containing MM10, MM14, MM4, and MM9), the magnetic field seems to be perpendicular to the filament spine, which is in agreement with previous observational studies of magnetic fields in filamentary IRDCs (Liu et al. 2018;Soam et al. 2019;Tang et al. 2019;Añez-López et al. 2020) and might be a result of gravitational accretion flows (Gómez et al. 2018;Li et al. 2018).In the northwestern clumps (MM1, MM2, MM11, and MM16) and in the diffuse region to the southeast of MM4 and MM9, the magnetic field tends to be parallel to the main straight dark filament, which may suggest that the magnetic field has kept its initial configuration inheriting from the large-scale magnetic field in these regions.
Figure 2 shows the magnetic field orientation in clumps MM4 and MM9 revealed by ALMA.Overall, our three configurationscombined (C-1, C-3, and C-4) ALMA images exhibit magnetic field morphologies similar to the two configurations-combined (C-1 and C-3) ALMA images reported by Liu et al. (2020).

Comparing multi-scale magnetic fields
One question to be addressed is how the small-scale magnetic field is correlated with the large-scale magnetic field (e.g., Li et al. 2009;Zhang et al. 2014;Li et al. 2015a).To study this, we first compare the JCMT magnetic field orientation (θ JCMT ) in every pixel of the JCMT detection area (S/N(P I)>2) with the Planck magnetic field orientation (θ Planck ) at the nearest pixel of the Planck map and calculate their angular difference with approaches similar to Zhang et al. (2014).Figure 3 overlays the magnetic field orientations revealed by Planck and JCMT. Figure 4(a) shows the spatial distribution of the absolute angular difference between the Planck and JCMT magnetic field orientations.Figure 5(a) shows the histogram of the absolute angular difference between the Planck and JCMT magnetic field orientations.Although overall there is no strong relation between the Planck and JCMT magnetic field orientations (Figure 5(a)), the spatial distribution for their angular difference is not random (as mentioned above and see Figure 4   averaged magnetic fields (e.g., Li et al. 2009).To do this, we average the Planck polarization data (i.e., Q and U ) toward the cloud within a circle of 10 pc and the JCMT polarization detection toward each clump within a circle of 1 pc, and recalculate the polarization position angles.We find that half of the clumps (MM1, MM2, MM11, MM16, and MM17) within our studied region have averaged magnetic fields aligned within 30 • of the cloud-scale magnetic field, while the other half of the clumps (MM4, MM6, MM9, MM10, and MM14) have averaged magnetic fields misaligned at 60 • -90 • with respect to the cloud-scale magnetic field (see Figure 5(b)).The bimodal distribution implies that the clump-scale magnetic field is organized with respect to the uniform cloud-scale Planck field.Otherwise, a randomly-orientated small-scale magnetic field cannot produce the observed bimodal distribution (Zhang et al. 2014).Future larger-sample observational studies are required to understand whether the bimodal distribution for angular differences between cloud-and clump-scale averaged magnetic fields is ubiquitous, and future theoretical and numerical studies will be needed to understand the physical mechanism behind this bimodal distribution.
Similarly, we investigate the angular difference between the polarization position angles revealed by JCMT and ALMA.Figures 4(b) and 5(b) show the spatial distribution and histogram, respectively, for the angular difference between the JCMT and ALMA magnetic field orientations toward MM4.The local spatial distribution of the angular difference (Figure 4(b)) is complex and we refrain from describing them in detail.The histogram of angular difference (Figure 5(c)) shows that the magnetic field on scales of cores and condensations in MM4, as revealed by ALMA, is preferentially aligned with the clump-scale magnetic field revealed by JCMT.We further compare the averaged clump-scale and core-scale magnetic fields for core1−5 in MM4.To calculate the averaged magnetic fields of cores, we adopt the core positions in Zhang et al. (2009) and average the polarization detections of each core within a radius of 0.1 pc in a way similar to the derivation of the averaged clump-scale magnetic fields.We find  that most cores (i.e., core1, 2, 4, and 5) have averaged magnetic fields aligned within 20 • of the clump-scale magnetic field of MM4 (see Figure 5(d)).The preservation of magnetic field orientation could suggest that the magnetic field plays a crucial role in the collapse of this clump and in the formation of dense cores within (Li et al. 2015b).The clump-scale magnetic field is perpendicular to the chain of cores in MM4, which could suggest that the fragmentation in this clump is regulated by the magnetic field in the parental clump (Nakamura & Li 2008).We refrain from comparing the averaged condensation-scale magnetic fields with large-scale magnetic fields in MM4 because the fragmentation in this clump at condensation level is complex (Zhang et al. 2015;Kong 2019).On the other hand, the marginal polarization detection in MM9 does not support any statistics.With a rough comparison, we find that the angle between the clump-scale magnetic field of MM9 and the condensation-scale magnetic field in C1-Sa is ∼65 • .The magnetic field orientation revealed by Planck is shown as grey line segments.The average magnetic field orientation within each 1-pc clump is indicated as black line segments.Contour starts at 50 mJy beam −1 and continues at 150 mJy beam −1 .(b).Map of absolute angular differences (color scale) between magnetic field orientations revealed by JCMT and ALMA in MM4.The magnetic field orientation revealed by JCMT is shown as grey line segments.The average magnetic field orientation within each 0.1-pc core is indicated as black line segments.Contour levels are (±3, 6, 10, 20, 30, 40, 50, 70, 90, 110, 150, 180, 210, 250, 290, 340, 390, 450) ×σI .

Density and mass
Here we briefly describe the estimation of the gas density and mass from dust continuum and molecular line data.Then we investigate the mass-radius relation and density-radius relation.

Core
For the dense cores revealed by ALMA observations, we estimate the dust mass with where F ν is the flux density at frequency ν, d is the distance, κ ν = (ν/1THz) β is the dust opacity (Hildebrand 1983) in m 2 kg −1 , and B ν (T ) is the Planck function at temperature T .We adopt a dust emissivity index (β) of ∼1.5 (e.g., Beuther et al.Histogram of angular differences between magnetic field orientations revealed by Planck and JCMT in G28.34. (b).Histogram of angular differences between averaged cloud-and clump-scale magnetic fields in G28.34.(c).Histogram of angular differences between magnetic field orientations revealed by JCMT and ALMA in MM4. (d).Histogram of angular differences between averaged clump-and core-scale magnetic fields in MM4.
2007; Chen et al. 2007) and a gas temperature of 15 K (Wang 2018).Adopting a gas-to-dust ratio of Λ = 100 (Savage & Jenkins 1972), the gas mass can be estimated with M gas = ΛM dust .Then the gas column density is estimated with where µ H2 = 2.8 is the mean molecular weight per hydrogen molecule (Kauffmann et al. 2008), m H is the atomic mass of hydrogen, and A is the area.

Clump and cloud
We adopt the column density map (hereafter the COMB map) at intermediate scales from Lin et al. (2017).The COMB map was made by performing an iterative SED fitting procedure for multi-wavelength (70 µm to 850 µm) continuum data from Herschel, Caltech Submillimeter Observatory (CSO), JCMT, and Planck, using image combination techniques to recover extended emission for ground-based telescopes while preserving the high angular resolution.Details about basic data combinations and SED fitting procedures can be found in Lin et al. (2016).In particular, we used an updated combination procedure similar to Jiao et al. (2022), where the Planck map is first deconvolved based on an extrapolated model image.The updated combined image benefits from a better dynamical range to recover more extended emissions.The resolution of the resulting column density maps is 10 ′′ .The size of the COMB map is ∼600 ′′ (∼7 pc).With the column density map, the gas mass is derived using Equation 3.

Environmental gas
For the environmental gas surrounding IRDC G28.34, we calculate the gas column density from the Planck τ 353 map and the FCRAO-14m 13 CO (1-0) data.
For the Planck data, we convert τ 353 to the column density of hydrogen atoms (N H ) with the relation (Planck Collaboration et al. 2014, 2016) τ Note that N H accounts for the column density of both the atomic gas and molecular gas (Planck Collaboration et al. 2014, 2016) along the line of sight (LOS).The gas mass is given by where we assume the mean atomic weight per hydrogen atom is µ H ∼ 1.4.
For the 13 CO (1-0) data, we estimate the integrated line intensity with W 13CO = Σ N ch i I i ∆v ch , where I i , ∆v ch , and N ch are the line intensity, channel width, and the number of integrated channels, respectively.The 13 CO (1-0) line is integrated between 71 and 86 km s −1 to cover the main velocity component of the cloud (Beuther et al. 2020).Because the 13 CO (1-0) line peak brightness temperature is much smaller than the gas temperature, this line is very likely optically thin (Beuther et al. 2020).In the optically thin case, the upper state column density is given by (Goldsmith & Langer 1999) where k B , h, and A ul are the Boltzmann constant, Planck constant, and Einstein coefficient, respectively.Assuming local thermal equilibrium, the total column density of 13 CO is given by where g u , E u , and Z are the statistical weight of the upper state, the upper energy level, and the partition function, respectively.
Here we adopt T = 15 K as well.The values of A ul , g u , E u , and Z are adopted from the CDMS (Müller et al. 2001) and LAMDA (Schöier et al. 2005) databases.Thus, the column density of H 2 can be estimated using the standard abundance (Wilson et al. 2013) The standard abundance is only valid when N H2 < 5 × 10 21 cm −2 and could present a scatter of factor two to five (Wilson et al. 2013).Because CO isotopes could be depleted in dark clouds (Bergin & Tafalla 2007), N H2 toward IRDC G28.34 estimated with 13 CO may only be a lower limit.With the estimated gas column density, the gas mass is derived with Equation 3. We find that the ratio between 0.5N H estimated with the Planck τ 353 map and N H2 estimated with the 13 CO data is ∼10 toward the center of G28.34.Their ratio gradually increases to ∼10 3 toward the edge of the Planck map.As the Planck observations trace all the atomic and molecular gases along the LOS, it is reasonable that the column density traced by Planck should be higher than the molecular gas integrated within a specific velocity range.On the other hand, the systematic lower gas column densities of 13 CO might be due to the 13 CO depletion, the FCRAO-14m observation filtering out the large-scale emission at the extent of the off-position, and the uncertainties in the adopted 13 CO abundance.

Scaling relations
The mass-radius (M − r) and density-radius (n − r) relations are important properties of star formation regions (Larson 1981).With the estimated column density and gas mass, we investigate the M − r and n − r relations for the cores and clumps in IRDC G28.34 and in its surrounding gas.
With the gas mass derived from the Planck τ 353 map and from the FCRAO-14m 13 CO (1-0) data, we obtain the mass profile for the environmental gas within circles of different radii (from beam size to 20 pc) centered at G28.34.With the COMB mass map, we obtain the mass profile for each infrared dark clump within circles of different radii (from beam size to 2 pc) centered at the local continuum peak.For simplicity, we only consider clumps with both molecular line detection and sufficient dust polarization detection (i.e., MM4, MM6, MM9, MM10, and MM14) in our analysis.MM10 and MM14 are unseparated at the resolution of the COMB map and are thus considered as one clump in our analysis.We also obtain a mass estimation for the whole COMB map area (r ∼7 pc).At the resolution of ALMA, the fragmentation status is complicated in MM4 and MM9.Thus, we only report the total mass and effective radius ( A/π) for the ALMA area with S/N(I)>5.Similarly, we obtain the profile for the average column density within circles of different radii.Assuming the studied structures are spherical, the average number density is estimated from the average column density with n = 0.75N/r.Mass M (M ) M La rs on = 4 6 0 r 1 .9M = 1843.2r 1.59 10 1 10 0 10 1 10 2 Radius r (pc) Figure 6(a) shows the obtained mass-radius relation.An obvious trend is that the mass and radius for the ALMA data and the COMB data are positively correlated.With a simple minimum chi-squared power-law fit, we obtain M ∝ r 1.59±0.14 .Here we assume an uncertainty of a factor of 2 for the mass during the fitting.The ALMA data for MM4 is higher than the fitted power law, which might be due to the deviation from the assumed spherical structure in this clump.The fitted power law is higher than the empirical threshold for massive star formation (M lim = 870M ⊙ (r/pc) 1.33 , Kauffmann et al. 2010), indicating that IRDC G28.34 is potentially forming massive stars.The obtained power-law index of 1.59 is smaller than the value of 1.9 reported by Larson (1981).Note that our M − r relation presents the relation in an individual cloud, while Larson's third law presents the relation for an ensemble of cloud samples, so the two relations may not be comparable.The mass derived from the Planck map is higher than the fitted power law, while the mass derived from the FCRAO-14m 13 CO data is lower than the fitted power law.As discussed in Section 3.2.3,this deviation might be because the molecular mass is overestimated by the Planck map and underestimated by the FCRAO-14m 13 CO data.
Figure 6(b) shows the obtained average number density-radius relation.With a power-law fit for the ALMA data and the COMB data, we obtain n ∝ r −1.41±0.14 .Here we assume an uncertainty of a factor of 2 for the density during the fitting.The obtained power-law index of -1.41 is steeper than the value of -1.1 reported by Larson (1981).Similarly, our n − r relation in an individual cloud may not be comparable to Larson's third law which presents the relation for an ensemble of cloud samples.

Molecular line and velocity fields
Here we briefly overview the multi-scale velocity structures revealed by the FCRAO- 14m 13 CO (1-0), JCMT 13 CO (3-2) and  HCO + (4-3), and ALMA N 2 D + (3-2) line data.The information of each line is summarised in Table 1.The value of line parameters are adopted from the CDMS (Müller et al. 2001) and LAMDA (Schöier et al. 2005) databases.No information on the collision rate coefficient exists for N 2 D + in LAMDA, so we adopt the values for N 2 H + instead.A particular line tracer is only sensitive to gas with densities above its critical density (n c ) but no more than a factor of 2 orders of magnitude (Goodman et al. 1998).We calculate the intensity-weighted velocity centroid V c (x) at position x with where I i (x), v i , ∆v ch , and N ch are the line intensity, line-of-sight velocity, velocity channel width, and number of channels, respectively.The propagated uncertainty of where σ ch is the noise of one spectral channel (see Section 2).For FCRAO-14m 13 CO (1-0), we only consider line emissions from 71 to 86 km s −1 which covers the main velocity component of G28.34 at large scales (Beuther et al. 2020).There are no apparent LOS distant gas structures at this velocity superposed at the same POS position for this IRDC (Simon et al. 2006;Beuther et al. 2020).For JCMT observations toward the cloud, we only consider velocities within ∼5 km s −1 with respect to the local-standard-of-rest (LSR) velocity (V lsr ∼ 79 km s −1 ) of G28.34 to avoid potential contamination from high-velocity outflows.For ALMA observations toward MM4 and MM9, we only consider velocities within ∼3 km s −1 with respect to the V lsr of each clump (∼79 and 80 km s −1 for MM4 and MM9, respectively, Zhang et al. 2015) because there is nearly no line emission at >3 km s −1 .All lines are likely to be optically thin as the main-beam temperature of the line peaks is smaller than the gas/dust temperature of ∼15 K (see Appendix B).

Velocity centroid maps
Figures 7 and 8 show the multi-scale velocity structure surrounding and within G28.34.The FCRAO-14m 13 CO (1-0) data (see Figure 7(a)) reveals a velocity gradient across Galactic latitudes near G28.34 (Beuther et al. 2020).The origin of this global gradient is still unclear.At higher resolution, the JCMT 13 CO (3-2) data reveals a velocity gradient perpendicular to the main dark filament (see Figure 7(b)) that is consistent with the global velocity gradient.Previous studies have interpreted this gradient as gas flows converging at the position of the filament (Beuther et al. 2020).The origin of the converging flow may be related to gravitationally driven collapse or external compression (Beuther et al. 2020).The HCO + (4-3) line with a higher critical density than the 13 CO (3-2) transition reveals the velocity structure of most clumps in G28.34 except for MM17 (see Figure 7(c)).Within the dark filament, the higher-density gas traced by the HCO + (4-3) line is more redshifted than the lower-density gas traced by the 13 CO (3-2) transition.With the NH 3 observations, Zhang et al. (2015) have identified a longitudinal velocity gradient of 0.6 km s −1 along the main dark filament from MM10 to MM14 to MM4, which reverses from MM4 to MM9.This velocity gradient is also seen in the JCMT 13 CO (3-2) and HCO + (4-3) data.Toward MM1, there seems to be a north-south velocity gradient deviating from the trend of the large-scale gradient, which may be an indicator of local gravitational infall in this clump.At smaller scales, the variation of V c is overall small in dense cores in MM4 and MM9 as seen by the ALMA N 2 D + (3-2) observations.There is a redshifted velocity component toward the streamer-like continuum structure in the north of MM4-core5, which may indicate that the core is accreting gas from the clump gas through the streamer.There is a possible sign of a redshifted velocity component in the north of C1-N.Other than these redshifted velocity components, there are no apparent signs of ordered velocity gradients in the major part of cores in MM4 and MM9.

Multi-scale statistics
Here we use the velocity dispersion-radius (σ v − r) relation and velocity centroid dispersion-radius (σ c − r) relation to explore the multi-scale velocity statistics in G28.34.
In a similar way to the derivation of the M − r relation (Section 3.2.4),we derive the σ v − r relation for the average spectra (see Appendix B) of each line within circles of different radii, except for the ALMA data where the spectra are averaged within the whole emission area.Similar to the approaches that most previous studies have adopted, we fit a single-component Gaussian profile for each line to obtain the Full Width at Half Maximum (FWHM.i.e., linewidth).The relation between the FWHM and velocity dispersion is F W HM = √ 8 ln 2σ v .For the JCMT data, we focus on the infrared dark clumps MM4, MM6, MM9,  , 6, 10, 20, 30, 40, 50, 70, 90, 110, 150, 180, 210, 250, 290, 340, 390, 450) ×σI .The synthesized beam and a scale bar are indicated in the lower left and lower right corner of each panel, respectively.
and MM10 (+MM14).The thermal velocity linewidth (<0.06 km s −1 for our studied molecules) is much smaller than the nonthermal linewidth in our studied regions due to the low temperature, so we regard the total linewidth as the non-thermal linewidth and do not subtract the thermal contribution.It should be noted that the radius of each considered region does not necessarily reflect the corresponding spatial scale of the measured linewidth.This is because the linewidth is also dependent on the LOS depth, which cannot be easily determined from observations.As we only fit a single-component Gaussian profile for all the lines, the linewidth in some regions could be broadened by the superposition of substructures with different LOS local-standard-of-rest (LSR) velocities at small scales (i.e., indistinguishable multiple velocity components in the line profiles.See Appendix B).Due to these effects, the measured linewidth from the observed spectra likely presents the upper limit of the true non-thermal linewidth at the corresponding scale of the radius.As seen in Figure 9(a), the σ v − r relation for each region seems to be relatively flat with small slopes or even negative slopes for the single-dish (FCRAO-14m and JCMT) observations.The flat slopes suggest that our measured velocity dispersion is likely dominated by the LOS integration up to the largest depth of the gas traced by specific lines, which is similar to what was found in the low-mass cloud Polaris Flare (Ossenkopf & Mac Low 2002).The JCMT velocity dispersions are lower than the FCRAO-14m 13 CO (1-0) velocity dispersions.This might be due to the higher n c of the 13 CO (3-2) and HCO + (4-3) lines, which trace fewer substructures along the LOS.On the other hand, the single-dish observations filter out the diffuse emission at the extent of the off position.The particular JCMT observations may have set a nearer off position than the FCRAO-14m observations, which could filter out more diffuse emissions and is also plausible to account for the lower velocity dispersions revealed by JCMT.For the JCMT observations, the velocity dispersions of HCO + (4-3) are lower than that traced by 13 CO (3-2), which might be due to the higher n c and fewer LOS substructures traced by HCO + (4-3).In addition, the ALMA velocity dispersions are much lower than those traced by JCMT and FCRAO-14m, which may mainly be because the emissions from lower-density substructures are filtered out by the interferometric observations.
We also study the scaling relation of velocity fields with the dispersion (i.e., standard deviation) of velocity centroids (σ c ).The velocity centroid is insensitive to thermal broadening and thus primarily reveals non-thermal motions (Lazarian et al. 2022).Moreover, velocity centroids present averaged statistics along the LOS and are not affected by the broadening effect of LOS depth as for the linewidth.Thus the radius in the σ c − r relation correctly reflects the true scale in the scaling relation.Due to the LOS averaging, for a single-component turbulent structure, the turbulent velocity centroid dispersion is smaller than the turbulent velocity dispersion by a factor of (Dickman & Kleiner 1985), where N turb is the number of independent turbulent cells along the LOS, L is the effective depth, and λ c is the turbulent correlation length.Previous observational V c studies usually implicitly assume that N turb is constant or does not vary too much within an individual cloud (e.g., Ossenkopf & Mac Low 2002;Liu et al. 2023b).Considering that each line tracer mainly traces a layer of gas with densities between n c and 100n c , the validity of this assumption could be unclear.
Figure 9(b) shows the velocity centroid dispersion-radius relation.The σ c − r relation in each region is complex and we refrain from interpreting them in detail in this work.The σ c − r relations from different instruments and lines are not continuous as well.The σ c in most clumps traced by JCMT is comparable to σ c in the environmental gas traced by FCRAO-14m, but lower than the σ c of the dense cores traced by ALMA.If N turb = L/λ c is similar for different line data or N turb is larger for higher-resolution data, the observed σ c − r relation could suggest that the early-stage massive star formation activities have already increased the non-thermal motions in small-scale and high-density regions in IRDC G28.34.However, we cannot rule out the possibility that the velocity centroid dispersion at a smaller r is less averaged due to a smaller N turb .
In summary, we find the velocity statistics are not continuous at different scales revealed by different instruments and lines, nor universal in different regions.However, due to the unsolved issues on the LOS length of linewidth and the averaging of velocity centroids, we cannot rule out the possibility that the actual multi-scale velocity statistics in G28.34 could still follow a continuous and universal power-law relation.Solving those issues requires the application of more advanced analysis methods (e.g., Lazarian & Pogosyan 2006).

Magnetic field strength
The Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951;Chandrasekhar & Fermi 1953) and its modified forms have been widely used to estimate the plane-of-sky (POS) magnetic field strength in molecular clouds.The basic assumptions of the DCF method are: there is an equipartition between the transverse (i.e., perpendicular to the ordered field) turbulent magnetic and kinetic energies; the turbulence is isotropic; and the turbulent-to-ordered or turbulent-to-total magnetic field ratio can be traced with the statistics of magnetic field orientations.Under these assumptions, the POS ordered and total magnetic field strengths are given by respectively, where µ 0 is the permeability of vacuum, ρ = µ H2 m H n H2 is the gas density, n H2 is the volume density, and σ v is the line-of-sight turbulent velocity dispersion.In small angle approximation (i.e., the ordered magnetic field is prominent), the POS turbulent-to-ordered field ratio (B t /B 0 ) or turbulent-to-total field ratio (B t /B) are usually estimated with the angular dispersion (σ θ ) of POS magnetic field orientations: A recent review of the DCF method can be found in Liu et al. (2022a,b).

Environmental gas
We estimate the magnetic field strength in the environmental gas of IRDC G28.34 with the Planck polarization data and the FCRAO-14m 13 CO (1-0) data.As demonstrated later (see Section 3.5.2), the Planck polarization detection toward G28.34 mainly traces the emission from the surrounding gas of this IRDC.Thus it is appropriate to use the Planck polarization map to obtain the magnetic field information of the cloud environment.The radius of the polarization detection region must be ≳2 of the beam size to obtain meaningful statistics of the turbulent magnetic field (Liu et al. 2021).Thus we consider a circular region with a radius of r =15 pc centered at G28.34.The angular dispersion of the Planck magnetic field is σ θ = 3.9 • within the considered region.The almost straight POS magnetic field lines at 15-pc scale and the small angular dispersion imply B 0 ∼ B. As mentioned in Section 3.2.3, the density estimated from the 13 CO (1-0) data tends to be underestimated while the density estimated from the Planck map tends to be overestimated.Thus we adopt the density extrapolated from the power-law relation (n = 6398 cm −3 (r/pc) −1.41 ) fitted from the COMB data and ALMA data.At r = 15 pc, we obtain n H2 ∼ 144 cm −3 .Similar to the turbulent kinetic field, the turbulent magnetic field as well as its ratio with respect to the ordered or total field could be underestimated by a factor of √ N B if there are N B independent magnetic turbulent cells along the LOS (Zweibel 1990).The 3D unaveraged turbulent magnetic and kinetic fields within the considered region are not directly measurable from observations.As the DCF method assumes equipartition between turbulent magnetic and kinetic energies, it is reasonable to further assume that the turbulent magnetic and kinetic fields are averaged by a similar extent (i.e., N B ∼ N turb ) due to the LOS signal integration.Although N turb and N B may not be accurately measured from our observational data, the corrections of the LOS signal integration effect for the turbulent kinetic field (as traced by σ c ) and turbulent magnetic field (as traced by σ θ ) could cancel out with ( (Cho & Yoo 2016).Within our studied region, the velocity centroid dispersion is σ c = 0.54 km s −1 .Incorporating the correction for the LOS signal integration, the magnetic field strength of the environmental gas within the 15-pc region is given by B ∼ B 0 ∼ √ µ 0 ρσ c /σ θ ∼ 0.074 mG.Alternatively, using a model with n H2 ∼ 100 cm −3 and l = 8 pc, Ostriker et al. (2001) numerically investigated the uncertainties in the DCF estimation of clouds and derived a correction factor of 0.5 (hereafter the Ost01 correction factor).As the density of their simulation is comparable to that of the region of our interest, we apply the Ost01 correction factor for an alternative magnetic field strength estimation.
Adopting σ v = 3.93 km s −1 at r = 15 pc, we obtain B ∼ B 0 ∼ 0.5 √ µ 0 ρσ v /σ θ ∼ 0.27 mG.Note that the critical density of 13 CO (1-0) is an order of magnitude higher than the density of our considered region, so the turbulent velocity dispersion here may be underestimated, which might lead to underestimation of the magnetic field strength as well.
At the considered scale, additional uncertainty on the estimated field strength could come from the anisotropy of MHD turbulence when gravity is negligible (Lazarian et al. 2022).The correction for the turbulence anisotropy requires knowledge of the mean-field inclination and the fraction of turbulence modes.However, both parameters cannot be easily measured with our existing data.Thus we are unable to correct this effect in this work.Some recent theoretical and numerical studies suggest that in a non-gravitational and very sub-Alfvénic environment, the turbulent kinetic energy is in equipartition with the coupling-term magnetic field energy fluctuation rather than with the turbulent magnetic energy (Federrath 2016;Skalidis & Tassis 2021;Beattie et al. 2022), which could bring another uncertainty to the DCF method.However, it has been debated whether their proposed coupling-term energy equipartition is valid for molecular clouds (Lazarian et al. 2022;Li et al. 2022;Liu et al. 2022a,b).Firstly, the physical interpretation for including the coupling-term field in the energy equipartition is unclear.Even if the reference velocity is set as the cloud average velocity, the substructures moving at a different velocity from the cloud average velocity or the ordered non-turbulent motions (e.g., infall or rotation) due to star formation activities could still generate a significant amount of coupling-term velocity field fluctuation within molecular clouds.There is no reason to only consider the coupling-term magnetic field but not consider the coupling-term velocity field.Secondly, all the numerical works supporting the coupling-term energy equipartition have adopted the whole simulation-averaged mean magnetic field, which represents the interstellar medium (ISM) magnetic field at approximately the turbulence injection scale (∼100 pc), in the calculation of the coupling-term magnetic field.While it is probably fine to adopt the global mean field as the local mean field in an extremely sub-Alfvénic case (e.g., M A = 0.01, Beattie et al. 2022), the star-forming 10-pc clouds are only trans-to-sub-Alfvénic (e.g., M A ≳ 0.6, Hu et al. 2019).In such cases, the 10-pc local mean field could have significantly deviated from the 100-pc global mean field and thus it would be inappropriate to adopt the global mean field of the ISM in the calculation of the cloud energetics.In summary, more work needs to be done to understand the validity of the coupling-term energy equipartition assumption in molecular clouds.

Clump
For the dense clumps revealed by JCMT observations, the magnetic field shows non-linear ordered field structures due to gravitational effects and star formation activities, which overestimates the angular dispersion that should be only attributed to turbulence.In G28.34, all clumps have σ θ > 25 • within r = 1 pc, but using σ θ to estimate B is only valid when σ θ < 25 • in both low and high density regions (Ostriker et al. 2001;Liu et al. 2021).To account for the ordered field contribution, we use the angular dispersion function (ADF) method (Hildebrand et al. 2009;Houde et al. 2009Houde et al. , 2016)), a modified DCF method, to estimate the magnetic field strength of these clumps.
The ADF accounting for the ordered field contribution, beam-smoothing effect, and turbulent correlation effect is given by where ∆Φ(l) is the angular difference of two position angles separated by l, l W = l beam / √ 8 ln 2 is the standard deviation of the beam size, a ′ 2 l 2 is the second-order term of the Taylor expansion for the ordered field, and l δ is the turbulent correlation length.Here the POS turbulent-to-total magnetic energy ratio ⟨B 2 t ⟩/⟨B 2 ⟩ does not consider the effect of LOS signal integration.Note that the second-order term a ′ 2 l 2 is only valid to represent the ordered field at small l when the higher-order terms are negligible.Also note that polarization observations can only trace the magnetic field orientation with a 180 • ambiguity, so ∆Φ(l) is constrained to be within [−90, 90] degrees (Hildebrand et al. 2009;Houde et al. 2009).Because the actual magnetic field direction angle in the range of -180 • to 180 • is approximated by the position angle in the range of -90 • to 90 • , the ADF method implicitly assumes that the turbulent field is smaller than the ordered field (i.e., sub-Alfvénic).
Figure 10 shows the ADFs for clumps MM4, MM6, MM9, and MM10.We adopt a binning interval of l bin = l beam /2 for the distance lag.We fit each ADF via reduced χ 2 minimization.We fit the ADFs over different maximum l ranges and obtain the best fit with the smallest reduced χ 2 .The best fitted results are ⟨B 2 t ⟩/⟨B 2 ⟩ 0.5 = 0.25, 0.26, 0.31, 0.22 for MM4, MM6, MM9, and MM10, respectively.
Due to the effect of LOS signal integration, the statistics of POS polarization position angles could underestimate the turbulent magnetic field.The ADF method proposes that the turbulent magnetic energy is averaged by a factor of for single-dish data due to the LOS integration (Houde et al. 2009), where l ∆ is the LOS cloud effective depth.Houde et al. (2009) suggested that l ∆ could be estimated as the width at half of the maximum of the normalized autocorrelation function for the integrated normalized polarized flux.However, the normalized autocorrelation function of the integrated normalized polarized flux is greater than half of the maximum for all the clumps, so l ∆ in our case cannot be derived in this way.Moreover, the numerical study by Liu et al. (2021) found that the ADF method may not work well for the effect of LOS signal integration.So we refrain from adopting the correction for this effect suggested by Houde et al. (2009).Alternatively, we adopt the analytical corrections suggested by Cho & Yoo (2016) (CY16) or numerical corrections by Liu et al. (2021) (Liu21).Cho & Yoo (2016) suggested using the turbulent velocity centroid dispersion instead of the turbulent velocity dispersion in the DCF equation to account for the LOS averaging effect.Here we derive σ c from the 13 CO (3-2) line data because its critical density is closer to the clump densities.The velocity centroid dispersions are σ c =0.34, 0.30, 0.22, and 0.23 km s −1 for MM4, MM6, MM9, and MM10, respectively.Adopting n = 6.3 × 10 3 , 5.8 × 10 3 , 5.7 × 10 3 , and 6.0 × 10 3 cm −3 , we obtain B ∼ √ µ 0 ρσ c (⟨B 2 t ⟩/⟨B 2 ⟩) −0.5 ∼0.081, 0.065, 0.039, and 0.060 mG for MM4, MM6, MM9, and MM10, respectively.Here we assume that the velocity centroid dispersion is dominated by turbulent motions rather than non-turbulent motions.The magnetic field strength could be overestimated if σ c is mostly non-turbulent.On the other hand, Liu et al. (2021) has numerically derived the average correction factors for the clump-to-core scale magnetic field strength estimated with the ADF method using the velocity dispersion.The velocity dispersions are σ v =2.32, 2.54, 1.91, and 2.70 km s −1 for MM4, MM6, MM9, and MM10, respectively.Adopting a numerical correction factor of 0.21 (with a 45% uncertainty, Liu et al. 2021), we obtain B ∼ 0.21 √ µ 0 ρσ v (⟨B 2 t ⟩/⟨B 2 ⟩) −0.5 ∼0.120, 0.120, 0.075, and 0.154 mG for MM4, MM6, MM9, and MM10, respectively.The magnetic field strengths of the clumps estimated with either the CY16 or Liu21 corrections are comparable to the magnetic field strength of the environmental gas, which means the magnetic field does not significantly scale with density in the low-density gas.This behavior is consistent with previous magnetic field strength estimations in other regions (Pattle & Fissel 2019;Liu et al. 2022a,b).Using the DCF methods to derive the ordered field strength could have large uncertainties in self-gravitating regions (Liu et al. 2021), so we only derive the total field strength for the clumps.
Within self-gravitating clumps and cores, the uncertainty from anisotropic turbulence on the DCF method is negligible (Liu et al. 2022b).Thus we do not consider the correction for this effect.Because there is no evidence supporting coupling-term energy equipartition in self-gravitating regions (Liu et al. 2022b), we keep our assumption of turbulent energy equipartition.Note that if the clumps are super-Alfvénic, the turbulent magnetic energy could be smaller than the turbulent kinetic energy, and our derived magnetic field strengths for the clumps could be overestimated.

Core
From our ALMA observations, only one core (MM4-core4) has sufficient polarization detections to derive the magnetic field strength.Liu et al. (2020) has estimated the ordered magnetic field strength of MM4-core4 using the ADF method.The kinetic information adopted in Liu et al. (2020) was from the EVLA NH 3 line data (Wang et al. 2012), which had a different resolu- tion and filtering scale larger than the ALMA data.Here we recalculate the magnetic field strength with the CY16 and Liu21 corrections and with the velocity information from the ALMA N 2 D + line data.In MM4-core4, the velocity dispersion and velocity centroid dispersion of N 2 D + are 0.61 and 0.37 km s −1 , respectively.We adopt the radius (r = 0.053 pc), density (n = 1.1 × 10 6 cm −3 ), and turbulent-to-total magnetic field strength ratio without correction for LOS integration ((⟨B 2 t ⟩/⟨B 2 ⟩) 0.5 ∼ 1. i.e., the total field is dominated by the turbulent field) estimated in Liu et al. (2020).Adopting the CY16 and Liu21 corrections, we obtain B ∼ √ µ 0 ρσ c (⟨B 2 t ⟩/⟨B 2 ⟩) −0.5 ∼0.29 mG and B ∼ 0.21 √ µ 0 ρσ v (⟨B 2 t ⟩/⟨B 2 ⟩) −0.5 ∼0.10 mG, respectively.It should be noted that both the CY16 and Liu21 corrections may not well account for the LOS integration effect at <0.1 pc and n > 10 6 cm −3 (Liu et al. 2021).Also, note that the turbulent velocity motions may be overestimated due to the existence of non-turbulent non-thermal motions.And if the core is super-Alfvénic, the energy equipartition assumption could break down.Those above-mentioned factors could all overestimate the B strength for MM4-core4.On the other hand, the ALMA observations filter out the large-scale emission, but Liu et al. (2021) did not perform the filtering when deriving the velocity dispersion from the simulations.This inconsistency could underestimate the B strength estimated from the Liu21 correction when using the line velocity dispersion from interferometric observations, which could explain the smaller B value from the Liu21 correction than that from the Cy16 correction.Overall, the B value in core MM4-core4 is comparable to or only slightly larger than that in clump MM4, which may suggest that gravity does not significantly compress and amplify the magnetic field strength in the clumps in the early massive star formation stage.

Alfvénic Mach number
Using the angular dispersion of polarization position angles to derive the Alfvénic Mach number M A requires corrections for the LOS signal integration and other effects.It is inappropriate to just equal the angular dispersion and M A as some previous studies have done.Based on the property of MHD turbulence (Goldreich & Sridhar 1995), Lazarian et al. (2018) proposed an alternative method (Velocity Gradient Technique, VGT) to estimate M A in non-self-gravitating regions with the statistics of velocity gradient orientation (θ VG ) from molecular line observations.The advantage of VGT is that the velocity gradient orientation histogram is independent of the LOS integration (Lazarian et al. 2018).It was numerically shown that the top-tobottom ratio (N top /N bot ) of the velocity gradient orientation histogram has the relation with M A in a sub-Alfvénic region (Hu et al. 2021) The Top (N top ) and Bottom (N bot ) values of the histogram can be obtained by fitting a Gaussian profile (i.e., (N top − N bot ) exp(−α(θ VG − θ 0 ) 2 ) + N bot , where α and θ 0 are coefficients) for the V c histogram.
Figure 11 shows the histogram of the centered velocity gradient orientation (θ VG − θ 0 ) for the FCRAO-14m 13 CO (1-0) line data within the r = 15 pc region centered at G28.34.We fit the histogram with a Gaussian profile and obtain M A = 0.74, which is slightly sub-Alfvénic.
Within molecular clouds, the velocity gradient could be significantly affected by gravitational effects and star formation activities, which makes the property of velocity gradient statistics overall deviate from what is expected for pure MHD turbulence.Thus we refrain from applying the VGT to smaller scales.

Magnetic field orientation
The Planck polarization map probes all the LOS structures, which makes it hard to separate the emission from the cloud environment and the foreground/background.Here we independently derive the magnetic field orientation from the FCRAO-14m 13 CO (1-0) data with the VGT to confirm whether the Planck polarization toward G28.34 mainly traces emission from the cloud environment.The basic principle behind VGT is that the magnetized turbulence eddy is anisotropic and elongated along the magnetic field, so the velocity gradient of MHD turbulence should be perpendicular to the magnetic field (Goldreich & Sridhar 1995) in the absence of gravity.By selecting different velocity ranges, this technique can separate LOS velocity components corresponding to different distances, thus avoiding contamination from the foreground/background.
We implement the velocity channel gradients (VChG, Lazarian & Yuen 2018) on the FCRAO-14m data with approaches similar to Hu et al. (2022) and Hu & Lazarian (2023): (1) We convolve the line intensity map at each velocity channel with a Sobel kernel to obtain the intensity gradient map of each channel (i.e., channel gradient).Pixels with S/N<5 are not considered in this step.
(2) In the POS, we bin the Planck map over 3×3 pixels and perform sub-block averaging (Yuen & Lazarian 2017a) for the line channel gradients within the area of each binned Planck pixel.The ratio between the block size and the line spatial resolution is ∼6.7, which is similar to previous VGT studies (e.g., Hu et al. 2022) and is sufficient for the statistics required by sub-block averaging.With a circular Gaussian fit for the histogram of channel gradients, we find the averaged gradient orientation Ψ gs for each block.Note that due to sub-block averaging, the resolution of the spectral data is degraded to the block size (∼5 ′ ).( 3) Integrating channels along the LOS, we calculate pseudo-Stokes parameters Q g and U g with where I i is the block-averaged intensity.Here we integrate the velocity channels within three different velocity ranges: the cloud velocity range (from 71 to 86 km s −1 ), ±5 km s −1 (from 74 to 84 km s −1 ), and ±3 km s −1 (from 76 to 82 km s −1 ).
(4) The orientation of VChG is given by θ VChG = 0.5 arctan(U g /Q g ).The pseudo-magnetic field orientation is then given by θ B,VGT = θ VChG + 90  Figure 12 compares the magnetic field orientations inferred from the VChG with those traced by the Planck dust polarization.The two orientations are well aligned with each other in bright regions near G28.34.As we only select the cloud velocity range in the VGT analysis, the field orientation inferred from VChG should not contain contributions from the foreground/background.Moreover, we convolve the JCMT maps to the same resolution as Planck, then estimate the average polarization position angle within 3 ′ of the peak of the convolved total intensity map.The magnetic field orientation measured from the convolved JCMT map is only 9 • different from the Planck measurement and 7 • different from the VGT measurement at the nearest pixel (see Figure 12).The close alignment between measurements from VChG, JCMT, and Planck tends to suggest that the Planck dust polarization toward G28.34 mainly traces the emission from G28.34 and its close environment, but not the foreground/background.Yuen & Lazarian (2017b) suggested that the magnetic field orientation and the block-averaged velocity gradient could transit to a perpendicular alignment (i.e., re-rotation) when gravity becomes important in high-density regions.Although there are some slight angular differences between the field orientations revealed by VChG and Planck, overall there is a lack of evidence for perpendicular alignments toward the brightest pixels.This suggests that gravitational motions do not dominate in the large-scale diffuse gas.On the other hand, the field orientations traced by VChG and Planck show larger differences in weak-emission regions, which may suggest that the Planck dust polarization mainly traces the materials outside of the cloud velocity range in those regions or those weak-emission areas are significantly influenced by noise.
In Section 3.4.1,we used the angular dispersion measured from the Planck polarization data to estimate the magnetic field strength of the cloud environment.Alternatively, if we use the angular dispersion from the VGT measurements in the DCF calculations, the estimated magnetic field strength will decrease by a factor of ∼3.However, this alternative field strength estimation should be regarded as a lower limit.This is because the angular dispersion from VGT is expected to be affected by ordered velocity fields.i.e., although not dominant, the gravitational motions may have already contributed to the orientations of VChG, which increases the angular dispersion of VGT measurements.

Relative orientation analysis
Here we investigate the multi-scale physical properties of G28.34 with a synergistic local relative orientation analysis (Liu et al. 2023a) combining the approaches of the polarization-intensity gradient (Koch-Tang-Ho or KTH) method (Koch et al. 2012a) and the Histogram of Relative Orientation (HRO) analysis (Soler et al. 2013).Basically, we characterize the relative angle between magnetic fields and other orientations (column density gradient and local gravity) with the alignment measure (AM) parameter (Lazarian & Yuen 2018).The AM is given by AM = ⟨cos(2ϕ o2 o1 )⟩, where ϕ o2 o1 = |θ o1 − θ o2 | is the angle between two orientations and is in the range of 0 to 90 • .The uncertainty of ϕ o2 o1 is given by δϕ o2 o1 = δθ 2 o1 + δθ 2 o2 .Data points with δϕ > 10 • are excluded from our analysis.The uncertainty of AM is given by (Liu et al. 2023a) where n ′ is the number of data points considered.For the Planck data, we consider area within r = 15 pc centered at G28.34.
For the JCMT data, we consider area with S/N(I)>5 and S/N(P I)>2.For the ALMA data, we consider area with S/N(I)>3 and S/N(P I)>2.We calculate AM in 6, 10, and 15 different column density bins for the Planck, JCMT, and ALMA data, respectively.The typical numbers of pixels are ∼20, ∼20-80, and ∼30-180 per bin for the Planck, JCMT, and ALMA data, respectively.

Magnetic field versus column density gradient
Figure 13 shows the AM-N relation for the angle between the magnetic field (θ B ) and column density gradient (θ NG ).The column density gradient is obtained by applying a 3×3 Sobel kernel to the column density map.The uncertainty of θ NG is derived following Liu et al. (2023a).As θ NG is perpendicular to the column density contour, studying the AM-N relation for ϕ N G B is equivalent to the HRO analysis (Soler et al. 2013).As seen in Figure 13, the relative orientation changes from a statistically slightly more perpendicular alignment (AM N G B ≲ 0) in the environment at low column densities as revealed by Planck to a slightly more parallel alignment (AM N G B ≳ 0) in G28.34 at intermediate column densities as revealed by JCMT.This transition can only be reproduced in trans-to-sub-Alfvénic simulations in numerical HRO studies (see a review of the HRO studies in Liu et al. 2022b), which suggests that G28.34 is in a trans-to-sub-Alfvénic environment.This is in agreement with the result of our VGT analysis in Section 3.5.The reasons for the different alignment at different column densities and the transition of alignment are still under debate (Liu et al. 2022b).As discussed below (see Section 3.6.2),the transition to AM N G B ≳ 0 may be related to the distortion of gravity.At high column densities revealed by ALMA, the two angles tend to transit back to AM N G B ∼ 0 as N increases, which may be due to the influence of early massive star formation activities (e.g., infall, rotation, outflows, accretion, and et al.).Overall, the alignment between θ B and θ NG in IRDC G28.34 at different column densities is very similar to what we have found in the evolved massive star formation region NGC 6334 (Liu et al. 2023a).18) between magnetic field (θB) and column density gradient (θNG) as a function of column density for Planck (left), JCMT (middle), and ALMA (right) observations.Due to the filtering of large-scale emissions for JCMT and ALMA, the absolute column densities from different instruments are not comparable.AM > 0 and AM < 0 indicate a preferentially parallel and perpendicular alignment, respectively.

Magnetic field versus local gravity
Figure 14 shows the AM-N relation for the angle between the magnetic field (θ B ) and local gravity (θ LG ).Considering gas mass at S/N(I)>3, the 2D local gravity direction (θ LG ) at position r j is calculated with the standard formula of gravitation: where e ij is the direction between position r i and r j , G is the gravitational constant, m j is the mass at position r j , and θ LG is the direction of g j (r).Overall, the AM-N trend for ϕ LG B is similar to the trend for ϕ N G B , where AM LG B transits from a statistically more perpendicular alignment to a statistically more parallel alignment, then transit back to AM LG B ∼ 0 as N increases.As ϕ LG B directly traces the correlation between magnetic fields and gravity, the observed AM-N trend for ϕ LG B suggests that the magnetic field is resisting gravitational distortion at lower density, but is dragged by gravity as density increases within the cloud, and might be affected by early massive star formation activities near the central young stellar objects at even higher densities (Liu et al. 2023a).LIU ET AL.

Normalised mass-to-flux ratio
Based on ideal MHD equations, Koch et al. (2012b) suggested that the normalized mass-to-flux ratio can be estimated with where Σ B is the local ratio between the magnetic field force (F B ) and the gravitational force (F G ). Σ B can be estimated with LG sin(90 if the hydrostatic gas pressure is negligible. Figure 15 shows the λ KTH derived from the KTH method.For the Planck data, there is λ KTH ∼ 1, which suggests a magnetically trans-critical state in the environment.For the JCMT data, the λ KTH increases with N and transits from λ KTH < 1 to λ KTH > 1.The discrepancy between the Planck and JCMT data might be because the JCMT data filters out the large-scale emission and underestimates the gravitational force at lower column densities.For the ALMA data, overall there is λ KTH ≳ 1.However, the magnetic field may be affected by star formation feedback from central young stellar objects, which could make the KTH method not applicable to the ALMA data.Note that the systematic uncertainty of the λ KTH value estimated with the KTH method is unclear due to the lack of numerical tests.Different star formation theories (e.g., Bonnell et al. 1997;McKee & Tan 2002) have distinct predictions over the energy balance between gravity, magnetic fields, and turbulence in different scales and evolutionary stages of massive star formation.Observationally determining the energy budget of massive star formation regions is key to distinguishing between those theories.Here we use the virial theorem to characterize the multi-scale equilibrium state of G28.34.Neglecting the surface term, tension term, and ordered velocity motion, the virial theorem is written as where I is the moment of inertia, E G is the gravitational energy, E th is the thermal energy, E turb is the turbulent kinetic energy, and E B is the magnetic energy.
where k i = (5−2i)/(3−i), R is the radius, and G is the gravitational constant.The thermal energy is given by E th = 1.5nkB T V , where V = 4πR 3 /3 is the Volume.The turbulent kinetic energy is given by E turb = 1.5M σ 2 turb , where σ turb is the 1D turbulent velocity dispersion.The magnetic energy is given by E B = B 2 V /(2µ 0 ).

Gravity vs Thermal force
The relative importance between gravity and the thermal term in the virial theorem can be characterized with the ratio 2E th /|E G |. Using the fitted mass profile M ∝ r 1.59 and density profile n ∝ r −1.41 (see Section 3.2) and assuming a constant temperature (T = 15 K), we have derived this thermal-to-gravitational ratio as a function of radius.Figure 16(a) shows the derived thermal-to-gravitational ratio as a function of radius r.The thermal force is negligible compared to the gravitational force throughout the scales of our interest which is critical for star formation.This is reasonable because the thermal energy is usually less dominant compared to other forces in massive star formation (Tan et al. 2014).There is a trend that 2E th /|E G | increases as r decreases.At a very small scale near the central Young Stellar Objects, the thermal force might surpass the gravitational force and establish an equilibrium state solely supported by thermal pressure.Because the power-law indexes for the mass profile and density profile may change at small scales, higher-resolution observations are required to study whether a small-scale equilibrium between thermal force and gravity could be achieved, and if so, at what scale this equilibrium might occur.

Gravity vs Turbulence
The relative importance between gravity and turbulence is usually characterized by the turbulent virial parameter (Bertoldi & McKee 1992).For a spherical structure with density profile n ∝ r −i , the turbulent virial mass M turb is given by α turb < 1 suggests the turbulence cannot solely resist gravitational collapse (i.e., sub-virial), and vice versa.
We have calculated the turbulent virial parameter α turb for the multi-scale structures of G28.34.For the density profile, we adopt i = 1.41 (see Section 3.2).For the clumps and cores, we adopt the mass estimated in Section 3.2.For the environmental gas, we adopt the mass extrapolated from the fitted mass-radius relation M ∝ r 1.59 (see Section 3.2).As discussed in Section 3.3, the LOS velocity dispersion σ v tends to overestimate the actual 1D velocity dispersion at scale r because the LOS substructures are superposed and the LOS depth is larger than r, while the velocity centroid dispersion σ c tends to underestimate the actual 1D velocity dispersion at scale r due to the LOS averaging.As it is hard to derive the pure turbulent velocity dispersion corresponding to a specific scale (r) in molecular clouds, we adopt the LOS velocity dispersion σ v as the upper limit of the 1D turbulent velocity dispersion and the velocity centroid dispersion σ c as the lower limit, with the assumption that the non-turbulent part of σ c is much weaker than the turbulent part.For the JCMT line data, we only adopt σ v and σ c estimated from 13 CO (3-2) data as its critical density is closer to the clump densities.
Figure 16(b) shows the estimated turbulent virial parameter α turb as a function of radius r.It is clear that the cores in MM4 and MM9 are sub-virial.For the clumps, the cloud, and the environmental gas, the α turb estimated with σ v is super-virial, while the α turb estimated with σ c is sub-virial.Due to the uncertainties in the turbulent velocity dispersion, we are unable to determine whether these large-scale structures are turbulence-supported or not.

Gravity vs Magnetic fields
The relative significance of gravity and magnetic fields can be assessed using their energy ratio E B /|E G |.Alternatively, this comparison is more usually expressed by the normalized mass-to-flux ratio λ (Crutcher et al. 2004).For a spherical structure with density profile n ∝ r −i , λ is given by (Liu et al. 2022a) λ > 1 indicates magnetic fields cannot solely resist gravitational collapse (i.e., magnetically super-critical), and vice versa.For a spherical structure, the relation between λ and the magnetic-to-gravitational energy ratio is With the magnetic field strengths estimated in Section 3.4, we have calculated the energy ratio E B /|E G | as well as the normalized mass-to-flux ratio λ for the environmental gas within 15-pc, the clumps, and MM4-core4.In the calculation, the POS total magnetic field B is converted to the 3D total magnetic field B 3d with the relation B 3d ∼ B × 1.25 (Liu et al. 2022b).Similarly, we adopt i = 1.41 for the density profile and adopt the extrapolated mass and density for the environmental gas.Table 2 summarizes the calculated λ values.
Figure 16(c) shows the energy ratio E B /|E G | as a function of radius.As most previous studies have used λ to compare magnetic fields with gravity, we also show the estimated λ values as functions of r, n, and N in Figures 17(a)-(c).We see that λ increases with increasing density and decreases with increasing radius.This trend is consistent with that derived from the KTH method (see Figure 15).The environmental gas tends to be magnetically sub-critical, which means the magnetic field is strong enough to resist gravitational collapse in the large-scale diffuse gas.In the clumps/cores (as revealed by JCMT and ALMA), the state transits to super-critical, allowing gravitational collapse to happen.The transition of the magnetic critical state occurs at the cloud-to-clump scale.At the core scale, the magnetic support is notably weaker, even weaker than the thermal support (see Figures 16(a) and (c)).Note that although the magnetic field strengths estimated with the CY16 and Liu21/Ost01 corrections are different, the critical states derived with the two different corrections are consistent.Although it is hard to assess the uncertainty of the λ values estimated from the modified DCF methods due to the many assumptions adopted, it is reasonable to think that the systematic uncertainties from those assumptions may shift λ toward larger or smaller values, but does not significantly change the general trend of λ.The trend of increasing λ with density in G28.34 is in agreement with the λ − N trend from the previous DCF estimations in the literature (reviewed in Liu et al. 2022a,b), which is shown in Figure 17(d) for comparison.Note that the sample of the DCF compilation statistics in Liu et al. (2022a,b) are mostly from different regions, while our analysis presents the multi-scale magnetic critical state in an IRDC for the first time.The similarity between the multi-scale magnetic critical state in IRDC G28.34 and in other star formation regions may suggest that early massive star formation regions follow a similar evolution route as in other star formation regions of different evolutionary stages and masses.The general trend of λ appears to be consistent with magnetic field-controlled star formation theories (Mouschovias et al. 2006), where magnetically sub-critical clouds gradually form super-critical substructures that collapse.The increase of λ at higher densities may be due to the dissipation of magnetic flux (e.g., ambipolar diffusion or magnetic reconnection, Mouschovias & Ciolek 1999;Lazarian & Vishniac 1999), or mass accumulation along magnetic field lines.

Gravity vs Sum of competing forces
The total equilibrium state of star formation regions can be determined by comparing |E G | with the competing terms (2E th + 2E turb + E B ) in the virial equation.Similar to the turbulent virial parameter, we could define a total virial parameter as the ratio of energies: Alternatively, the total virial parameter can be defined as the ratio between the total virial mass and mass: Note that the values of α total,E and α total,M are not equal.The total virial mass is given by (Liu et al. 2020) where the thermal virial mass is estimated with and the magnetic virial mass4 is estimated with The environmental gas should have α total > 1 since the magnetic field itself is stronger than gravity.i.e., gravity is unimportant in the large-scale diffuse gas.At the intermediate scale, although it is clear that magnetic fields are dominated by gravity in the clumps (i.e., λ > 1), it is hard to assess their total equilibrium states (considering both the turbulent and magnetic supports).This is because there are some uncertainties in the derivation of the pure turbulent velocity dispersion at specific scales for the singledish data (see Section 3.3).Approximately adopting an average value between the LOS velocity dispersion and the velocity centroid dispersion as the turbulent velocity dispersion, we obtain that some clumps have α total < 1 while some clumps have α total > 1, with the average state being near quasi-equilibrium (i.e., α total not far from 1).Alternatively, assuming that the magnetic support is comparable to the turbulent support at the clump scale, we obtain similar results that the average state of clumps is near quasi-equilibrium.At the core scale, even when considering the upper limit of supporting forces, the total virial parameter for core MM4-core4 is only α total,E ∼ 0.5.Thus, it is safe to say that gravity dominates over the combination of competing forces for MM4-core4.
In summary, we suggest that the G28.34 cloud is located in a globally-supported environment and its clumps are likely in an approximate quasi-equilibrium state, but the cores therein are undergoing dynamic collapse.

Implications on massive star formation
It has been long debated whether molecular clouds and their substructures are in equilibrium or not.The Planck dust polarization survey (Planck Collaboration et al. 2016) has found that the Gould Belt Clouds, including one massive star formation region (Orion), are in magnetically sub-critical (i.e., λ < 1) states with DCF estimations.Later, Liu et al. (2023a) analyzed the Planck dust polarization data in another massive star formation region NGC 6334 and found that this region is also in a magnetically sub-critical state at large scale.Combined with our finding of magnetically sub-critical state in the environment of IRDC G28.34, we suggest that massive star-forming clouds may be globally supported by magnetic fields both in the early and late evolutionary stages.The quasi-equilibrium state of clouds is essential to explain the low star formation rate observed and to allow ambipolar diffusion to happen (McKee & Ostriker 2007), but contradicts the idea that clouds are short-lived and are undergoing global dynamical collapse (i.e., the global hierarchical collapse model, Vázquez-Semadeni et al. 2019).On the other hand, the sub-critical state does not necessarily mean the region will expand or disperse.Sub-critical clouds could still create overdense regions via other mechanisms such as large-scale turbulent inertial flows (i.e., the inertial-flow model, Padoan et al. 2020) or local infall through magnetic channels (Koch et al. 2018(Koch et al. , 2022)), instead of via symmetric gravitational collapse.
Within the G28.34 cloud, while the clumps may not be far from equilibrium, the cores are dominated by gravity.This is consistent with the findings of previous studies that gravity is more important in higher-density regions, and that cores in both early and evolved massive star formation regions tend to be averagely dominated by gravity even considering both the magnetic and turbulent supports (see reviews in Liu et al. 2022a,b).The dynamic state of cores is inconsistent with the turbulent core accretion model (McKee & Tan 2002), which predicts α total ∼ 1 across different scales.The near-equilibrium state of clumps tends to be inconsistent with the competitive accretion model (Bonnell et al. 1997), which requires α total < 1 for cloud substructures (Krumholz et al. 2005).Thus, both the two major massive star formation models may need some modifications to be consistent with the observational results.

Magnetic fields or turbulence?
Different star formation theories have distinct explanations for the controlling factor of the formation and evolution of molecular clouds and cloud substructures: some emphasize the role of magnetic fields (Mouschovias et al. 2006), while others highlight the role of turbulence (Mac Low & Klessen 2004).Qualitative and quantitative comparisons between magnetic fields and turbulence are needed to determine their relative importance in star formation.
There is no way to observationally compare the turbulent kinetic energy and the turbulent magnetic energy yet (Liu et al. 2022b).People usually assume an equipartition between turbulent kinetic and magnetic energies, which implicitly assumes that the total magnetic energy is larger than the turbulent kinetic energy.On the other hand, the relative importance between the ordered magnetic field and the turbulence can be directly derived from the dust polarization maps, without the need for molecular line observations.The relation between the 3D Alfvénic Mach number and the POS angular dispersion is M A ∼ (f t /f u /Q c /f o )σ θ (Liu et al. 2022b), where f u and f t are factors for the 3D-to-POS conversion of B 0 and B t , f o is a correction factor for the ordered field contribution, and Q c is a correction factor to account for other effects (e.g., the LOS signal integration, the difference between the orientation and direction, et al.).However, the actual value of f t , f u , and f o is usually unclear in individual regions.Statistical values exist for the correction factors (e.g., Crutcher et al. 2004;Liu et al. 2021), but those statistical values (especially for f u and f o ) may only be appropriate for statistical studies (e.g., Liu et al. 2022a;Pattle et al. 2023).Thus, we refrain from deriving M A from polarization angular dispersions in G28.34.
The relative orientation analysis (Section 3.6) offers an alternative way to qualitatively compare magnetic fields and turbulence.With an approach equivalent to the HRO analysis, we find that the magnetic field and column density gradient transits from statistically more perpendicular to more parallel as column density increases for the Planck and JCMT observations (Section 3.6.1),which is a sign of trans-to-sub-Alfvénic turbulence at large scales.This transition of alignment is consistent with previous studies in low-mass star formation regions and in evolved massive star formation regions (Planck Collaboration et al. 2016;Liu et al. 2023a).
The statistical tool VGT offers another way to compare magnetic fields and turbulence at large scales.Implementing the VGT, we find M A = 0.74 within r = 15 pc, which implies that G28.34 is located in a slightly sub-Alfvénic environment.The sub-Alfvénic state is consistent with the results from previous VGT studies in nearby low-mass star formation regions (Hu et al. 2019).
In summary, we conclude that both low-mass and high-mass star formation happens in a trans-to-sub-Alfvénic environment, which means magnetic fields play a more important role than turbulence in controlling star formation at large scales.Since gravity is not important in large-scale diffuse gas, cloud formation and evolution at the cloud scale should be mainly controlled by the property of trans-to-sub-Alfvénic MHD turbulence.Within the cloud, the situation is more complicated at small scales.
Figures 20 and 21 show the integrated intensity maps of the FCRAO-14m 13 CO (1-0), JCMT 13 CO (3-2) and HCO + (4-3), and ALMA N 2 D + (3-2) data.The integrated velocity ranges are indicated in each panel.Figures 22,23,and 24 show examples of the average spectra for different line data.Many lines show signs of multiple velocity components, which provides evidence for the superposition of substructures along the LOS.We do not try to identify or separate the multiple velocity components because the combination of those components could appear as a single Gaussian line profile at a coarser spatial resolution, and there are always multiple velocity components seen with higher resolution observations inside a "single" velocity component.We perform a single-peak Gaussian fit for each line and obtain the best-fitting result.The fitted velocity dispersion should be the upper limit of the pure turbulent velocity dispersion., 6, 10, 20, 30, 40, 50, 70, 90, 110, 150, 180, 210, 250, 290, 340, 390, 450) (a)).The difference between global and local statistics signifies the importance of investigating the local distributions of the angular difference.Thus, we further investigate the relation between the cloud-and clump-scale 18

Figure 1 .
Figure 1.(a).Magnetic field orientations revealed by Planck (black line segments) 0.85 mm dust polarization data overlaid on the Planck dust optical depth map (color scales) toward the IRDC G28.34 and its surrounding materials.Line segments are of arbitrary length.The 5 ′ (∼7 pc) beam (white circle) and a scale bar of 10 pc are indicated in the lower left and right corners, respectively.A white dashed line indicates the galactic plane (b = 0 • ).The dashed rectangle indicates the JCMT map area in (b).(b).Magnetic field orientations revealed by JCMT POL-2 0.85 mm dust polarization observations overlaid on the JCMT 0.85 mm total intensity map (color scales.Only area with S/N(I)>25 is shown) toward the IRDC G28.34.Black and grey line segments indicate S/N(P I)>3 and 2<S/N(P I)<3, respectively.Line segments are of arbitrary length.The 14 ′′ (∼0.33 pc) beam (black circle) and a scale bar of 1 pc are indicated in the lower left and right corner, respectively.The infraredbright and infrared-dark molecular clumps identified by Rathborne et al. (2006) are labeled with green and brown numbers, respectively.Red dashed contours indicate the ALMA fields of MM4 and MM9 corresponding to primary-beam responses of 0.5.Orange dashed circles mark the regions for MM4, MM6, MM9, and MM10 within which we use the polarization position angles to calculate the magnetic field strength in Section 3.4.2.

Figure 2 .
Figure 2. Magnetic field orientations revealed by ALMA 1.3 mm dust polarization observations overlaid on the ALMA 1.3 mm total intensity map (color scales) toward MM4 and MM9.Black and grey line segments indicate S/N(P I)>3 and 2<S/N(P I)<3, respectively.Line segments are of arbitrary length.The ∼ 0.7 ′′ × 0.5 ′′ (∼0.016-0.012pc) synthesized beam (black ellipse) and a scale bar of 0.1 pc are indicated in the lower left and right corners, respectively.

Figure 4 .
Figure 4. (a).Map of absolute angular differences (color scale) between magnetic field orientations revealed by Planck and JCMT in G28.34.The magnetic field orientation revealed by Planck is shown as grey line segments.The average magnetic field orientation within each 1-pc clump is indicated as black line segments.Contour starts at 50 mJy beam −1 and continues at 150 mJy beam −1 .(b).Map of absolute angular differences (color scale) between magnetic field orientations revealed by JCMT and ALMA in MM4.The magnetic field orientation revealed by JCMT is shown as grey line segments.The average magnetic field orientation within each 0.1-pc core is indicated as black line segments.Contour levels are(±3, 6, 10, 20, 30, 40, 50, 70, 90, 110, 150, 180, 210, 250, 290, 340, 390, 450)  ×σI .

Figure 5 .
Figure 5. (a).Histogram of angular differences between magnetic field orientations revealed byPlanck and JCMT in G28.34.(b).Histogram of angular differences between averaged cloud-and clump-scale magnetic fields in G28.34.(c).Histogram of angular differences between magnetic field orientations revealed byJCMT and ALMA in MM4.(d).Histogram of angular differences between averaged clump-and core-scale magnetic fields in MM4.

Figure 6 .
Figure 6.(a).Mass-radius relation.The red solid line indicates the empirical threshold for massive star formation (M lim = 870M⊙(r/pc) 1.33 , Kauffmann et al. 2010).The blue solid line indicates the Larson's third law (MLarson = 460M⊙(r/pc) 1.9 , Larson 1981).(b).Average number density-radius relation.The blue solid line indicates the Larson's third law (nLarson = 1586 cm −3 (r/pc) −1.1 , Larson 1981).Different colors represent different regions.Different symbols indicate different datasets.The black dashed lines present the results of a simple power-law fit for the ALMA data and the COMB data.

Figure 7 .
Figure 7. (a) Velocity centroid map of the FCRAO-14m 13 CO (1-0) data in the surrounding of G28.34.Black contours correspond to the Planck τ353 map, starting from 0.0005 and continuing with an interval of 0.0005.A dashed line indicates the galactic plane (b = 0 • ).The dashed rectangle indicates the JCMT map area in (b) and (c).(b)-(c) Velocity centroid maps of the 13 CO (3-2) and HCO + (4-3) emission toward the IRDC G28.34 obtained with JCMT.Black contours correspond to the JCMT 0.85 mm dust continuum map.Contours start at 50 mJy beam −1 and continue at a step of 150 mJy beam −1 .The beam and a scale bar are indicated in the lower left and lower right corner of each panel, respectively.
Figure 9(a) shows the velocity dispersion-radius relation.At first glimpse, it is obvious that the σ v − r relations from different instruments and lines are not continuous.The discontinuity of σ v − r relation is similar to what Peretto et al. (2023) have recently found in a large sample of 27 infrared dark clouds.While Peretto et al. (2023) interpreted the discontinuous σ v − r relation as clump dynamics decoupling from those of the parental clouds, we suggest that the discontinuity is more likely due to the fact that the physical conditions traced by different line data are different and that the σ v − r relation does not correctly present the true velocity dispersion-size relation of the gas in 3D.

Figure 9 .
Figure 9. (a).Velocity dispersion-radius relation.The blue solid line indicates Larson's first law (σ1D,Larson = 0.83 km s −1 (r/pc) 0.38 , Larson 1981).The original 3D velocity dispersion in Larson's first law is converted to 1D velocity dispersion with the relation σ3D,Larson = √ 3σ1D,Larson.(b).Velocity centroid dispersion-radius relation.Only data points with l greater than the beam resolution and smaller than the effective radius of the considered region are shown.Different colors represent different regions.Different symbols indicate different instruments and lines.

Figure 10 .
Figure 10.Angular dispersion functions for MM4, MM6, MM9, and MM10.The diamond symbols represent the observed data points.The error bars indicate the statistical uncertainties propagated from the observational uncertainty (Houde et al. 2009).The blue dashed line indicates the best-fitted results.The cyan dashed line shows the large-scale component (⟨B 2 t ⟩/⟨B 2 ⟩ + a ′ 2 l 2 ) of the best fit.The horizontal dashed line indicates the ADF value corresponding to a random field (0.36, Liu et al. 2021).

Figure 11 .
Figure 11.Histogram of the centered velocity gradient orientation for the FCRAO-14m 13 CO (1-0) line data within the r = 15 pc region.The blue line indicates the best-fitted result.

Figure 12 .
Figure 12.Pseudo magnetic field orientations inferred from VChG (red lines) within different velocity ranges (indicated in each panel) overlaid on the column density map.The Magnetic field orientations revealed by the Planck dust polarization are shown in black lines.The blue line indicates the average magnetic field orientation of the convolved JCMT polarization map.Line lengths are arbitrary.The white dashed circle marks the r = 15 pc region surrounding G28.34 within which we perform the DCF analysis.Within r = 15 pc, the pseudo magnetic field orientations derived from VGT within different velocity ranges are very similar.

Figure 13 .
Figure13.Relative orientations (characterized by AM .See Equation18) between magnetic field (θB) and column density gradient (θNG) as a function of column density for Planck (left), JCMT (middle), and ALMA (right) observations.Due to the filtering of large-scale emissions for JCMT and ALMA, the absolute column densities from different instruments are not comparable.AM > 0 and AM < 0 indicate a preferentially parallel and perpendicular alignment, respectively.

Figure 14 .
Figure 14.Relative orientations (characterized by AM ) between magnetic field (θB) and local gravity (θLG) as a function of column density for Planck (left), JCMT (middle), and ALMA (right) observations.

Figure 15 .
Figure 15.Normalised mass-to-flux ratio derived from the KTH method as a function of column density for Planck (left), JCMT (middle), and ALMA (right) observations.

Figure 16 .
Figure 16.(a).Thermal-to-gravitational ratio 2E th /|EG| as a function of radius.The horizontal dashed line indicates a balance between the thermal term 2E th and gravitational term EG in the virial equation.(b).Turbulent-to-gravitational ratio 2E turb /|EG| = α turb (i.e., turbulent virial parameter) as a function of radius.Different colors represent different regions.Different symbols indicate different instruments.Close symbols indicate α turb estimated with LOS velocity dispersion σv, while open symbols indicate α turb estimated with velocity centroid dispersion σc.The horizontal dashed line indicates a balance between the turbulent term 2E turb and gravitational term EG in the virial equation.(c).Magnetic-to-gravitational ratio EB/|EG| = 1/λ 2 as a function of radius.Different colors represent different regions.Different symbols indicate different instruments and modified DCF methods.The horizontal dashed line indicates a balance between magnetic energy and gravitational energy.

Figure 17 .
Figure 17.(a)-(c).Normalised mass-to-flux ratio derived from the DCF method as functions of radius, number density, and column density for Planck, JCMT, and ALMA observations toward G28.34.Different colors represent different regions.Different symbols indicate different instruments and modified methods.The horizontal dashed line indicates an energy balance between magnetic energy and gravitational energy (i.e., λ = 1).(d).Normalised mass-to-flux ratio derived from the compilation of DCF estimations of different star formation regions in the literature before June 2021 (reproduced with permission.See reviews in Liu et al. 2022a,b, and see references therein).Different colors represent different DCF variants and different symbols indicate different observational instruments (see Figure 3 of Liu et al. 2022b).

Figure 24 .
Figure 24.Average ALMA N2D + (3-2) spectra of MM4 and MM9 within the whole emission area in the ALMA field.The dashed Gaussian profiles indicate the best-fitting results.The vertical dashed lines indicate the velocity ranges used for the analyses.Our computing ability for ALMA data reduction is limited, so we only imaged the line within velocity ranges from 75 to 83 km s −1 .

Table 1 .
Summary of molecular line data Critical density at 15 K.For N2D + , we adopt the collision rate coefficient of N2H + .
d Maximum recoverable scale for ALMA.
For a spherical structure with density profile n ∝ r −i , the gravitational energy is given by E

Table 2 .
Summary of physical parameters