ABSTRACT

The relative importance of magnetic fields, turbulence, and gravity in the early phases of star formation is still not well understood. We report the first high-resolution dust polarization observations at 850 |$\mu$|m around the most massive clump, located at the hub of the Giant Molecular Cloud G148.24+00.41, using SCUBA-2/POL-2 at the James Clerk Maxwell Telescope. We find that the degree of polarization decreases steadily towards the denser portion of the cloud. Comparing the intensity gradients and local gravity with the magnetic field orientations, we find that local gravity plays a dominant role in driving the gas collapse as the magnetic field orientations and gravity vectors seem to point towards the dense clumps. We also find evidence of U-shaped magnetic field morphology towards a small-scale elongated structure associated with the central clump, hinting at converging accretion flows towards the clump. Our observation has resolved the massive clump into multiple substructures. We study the magnetic field properties of two regions, central clump (CC) and northeastern elongated structure (NES). Using the modified Davis–Chandrasekhar–Fermi method, we determine that the magnetic field strengths of CC and NES are ∼24.0 ± 6.0 |$\mu$|G and 20.0 ± 5.0 |$\mu$|G, respectively. The mass-to-flux ratios are found to be magnetically transcritical/supercritical, while the Alfv|$\acute{\text{e}}$|n Mach number indicates a trans-Alfv|$\acute{\text{e}}$|nic state in both regions. These results, along with Virial analysis, suggest that at the hub of G148.24+00.41, gravitational energy has an edge over magnetic and kinetic energies.

1 INTRODUCTION

Molecular clouds, the enigmatic birthplaces of stars, are now well established to consist of filamentary structures, revealed by Herchel observations (Molinari et al. 2010; André et al. 2013; Zavagno et al. 2023). Within these immense structures of cold gas and dust, the gravity, turbulence, and magnetic fields dictate the process of star formation at different scales, from large scale (cloud and filaments) to small scale (clumps and dense cores) (Klessen, Heitsch & Low 2000; Ballesteros-Paredes et al. 2007; Federrath 2015; Tang et al. 2019; Wang et al. 2020b; Pattle et al. 2022). However, their relative role during the different stages of cloud evolution is still unclear and a topic of debate (Li et al. 2014).

The magnetic field exists throughout star-forming molecular clouds across various scales (see the review article by Pattle et al. 2022) and plays a crucial role in the formation of molecular clouds and filamentary structures (Soler et al. 2013; Hennebelle & Inutsuka 2019). Numerical simulations show that the strong magnetic fields play an important role in the magnetically channelled gravitational collapse of clouds (Nakamura & Li 2008; Gómez, Vázquez-Semadeni & Zamora-Avilés 2018), and can channel turbulent flows along the filaments (Li & Houde 2008; Soler et al. 2013; Zamora-Avilés, Ballesteros-Paredes & Hartmann 2017), guide the accreting matter (Seifried & Walch 2015; Shimajiri et al. 2019), and dynamically influence the formation of cores along the dense ridges of the filamentary clouds (e.g. Koch et al. 2014; Zhang et al. 2014; Cox et al. 2016; Pattle et al. 2017; Soam et al. 2018; Liu et al. 2019; Eswaraiah et al. 2021). At parsec (or few parsec) scale, the magnetic field typically shows an ordered structure, mostly aligned with the long axis of the low-density elongated gas structures such as striations, while in high-density filaments, the magnetic field lines are preferentially perpendicular to the long axes of filaments (e.g. Cox et al. 2016; Planck Collaboration XXXV 2016; Soler et al. 2017; Ward-Thompson et al. 2017; Soam et al. 2019; Tang et al. 2019; Doi et al. 2020). At sub-parsec scales, the magnetic field can be very complex depending upon the turbulent nature of the magnetic field and stellar feedback (e.g. Hull et al. 2017; Eswaraiah et al. 2020, 2021). Dust polarization studies at small scales have revealed a variety of magnetic field morphologies in dense clumps and cores, and the results suggest that the magnetic field is scale-dependent and varies with the environment (e.g. Girart et al. 2013; Hull et al. 2017; Ward-Thompson et al. 2017; Pattle et al. 2018; Eswaraiah et al. 2020, 2021).

Moreover, it is not clear whether it is turbulence or magnetic field along with gravity that dominates the star formation process. The role of gravity has long been recognized as the primary factor driving the collapse of dense regions and initiating the birth of protostellar cores. On the other hand, turbulence, the chaotic, and ubiquitous motion of gas within these massive clouds, influences the fragmentation of the collapsing gas. However, the role of magnetic field, in comparison to turbulence and gravity, is relatively less understood at various stages of star formation.

The ‘strong magnetic field’ theory of star formation stresses the importance of the magnetic field in the formation and evolution of clouds and subsequent structures (Mouschovias, Tassis & Kunz 2006; Tan et al. 2013; Hennebelle 2018). Conversely, the ‘weak magnetic field’ theory suggests that turbulent flows control the formation and evolution of clouds and cores, and create the compressed regions where stars form (Padoan & Nordlund 2002; Mac Low & Klessen 2004; Federrath & Klessen 2012). Observationally also, in some high-mass star-forming regions, it has been found that the turbulent energy is more dominant or comparable to magnetic energy (e.g. Beuther et al. 2010, 2020; Girart et al. 2013; Wang et al. 2020b). In contrast, other studies found magnetic energy to be more dominant than turbulent energy (e.g. Girart et al. 2009; Beuther et al. 2018; Eswaraiah et al. 2020; Chung et al. 2023). Therefore, more observational evidence is required at the early stages of star formation to better constrain the theoretical models and relative roles of gravity, magnetic field, and turbulence in the star-forming regions.

Dust polarization observation is a key tool for tracing the plane-of-sky (POS) magnetic field geometry in star-forming regions. The polarization is caused by asymmetric dust grains that preferentially align their shorter axis along the magnetic field direction (Lazarian 2007; Hoang & Lazarian 2008). There are different mechanisms for dust grain alignment (see the review article by Andersson, Lazarian & Vaillancourt 2015), but the most widely accepted one being is the radiative alignment torque (RAT) mechanism (Lazarian 2007; Hoang & Lazarian 2014; Andersson, Lazarian & Vaillancourt 2015).

Located at a distance of ∼3.4 kpc, G148.24+00.41 is a massive cloud having mass ∼105 Mand dust temperature ∼14.5 K, and it resembles a hub-filament morphology (Rawat et al. 2023). The cloud is still at the early stages of its evolution such that the stellar feedback mechanisms like stellar winds and expanding H ii regions are not yet significant (Rawat et al. 2023), which makes G148.24+00.41 a potential candidate to study the early stages of the star and cluster formation process. Fig. 1(a) shows the 13CO intensity map of G148.24+00.41, integrated in the velocity range of −37.0 to −30.0 km s−1, where one can see a bright spot in the centre of the map. This location corresponds to the most massive clump of the cloud, which is nomenclatured as C1 in Rawat et al. (2024). Rawat et al. (2024), using CO molecular data (spatial resolution ∼52 arcsec or 0.9 pc at distance ∼3.4 kpc, velocity resolution ∼0.17 km s−1), also found filamentary features and seven massive clumps in G148.24+00.41. The C1 clump with a mass of ∼2100 Mand an effective size of ∼1.8 pc based on C18O data, is found to be at the hub/central location of the cloud. The filamentary features attached to the C1 clump can also be seen in Herschel 250 |$\mu$|m image, shown in Fig. 1(b) (adapted from Rawat et al. 2023). The figure clearly reveals the hub-filament system morphology of the cloud. Rawat et al. (2023) found that the clump hosts an infrared cluster, seen in Spitzer bands, and hypothesized that the cluster might grow to a richer cluster by accumulating gas from the extended reservoir. Such hub filamentary systems with the clump being located at the nexus or junction of these filaments are of particular interest because these are the sites where cluster formation would take place, as advocated in simulations and observations (e.g. Naranjo-Romero et al. 2012; Gómez & Vázquez-Semadeni 2014; Gómez, Vázquez-Semadeni & Zamora-Avilés 2018; Vázquez-Semadeni et al. 2019; Kumar et al. 2020).

(a) 13CO (J = 1–0) intensity map of G148.24+00.41, integrated in the velocity range −37.0 to −30.0 km s−1. (b) The central region (encompassed by the solid green box in panel-a) of G148.24+00.41as seen in Herschel 250 $\mu$m band, showing the hub-filamentary morphology of the cloud. The figure is adapted from Rawat et al. (2023), in which the blue circle shows the JCMT scanned region of diameter ∼12 arcmin (∼12 pc). The green dashed box marks the central area of the hub, where an infrared cluster is seen (Rawat et al. 2023), and the cross sign indicates the position of a massive young stellar object. (c) The 850 $\mu$m Stokes I intensity map of the central region of G148.24+00.41 mapped by JCMT SCUBA-2/POL-2, along with the contours of 13CO integrated intensity emission, drawn from 1.5 to 15 K km s−1with a step size of ∼0.96 K km s−1. The rms noise of the 4 arcsec pixel-size Stokes I map is around ∼5 mJy beam−1. In panel-c, the yellow ellipse shows the position of the C1 clump, identified by Rawat et al. (2024) using C18O data (spatial resolution ∼52 arcsec). The beam sizes of the 13CO integrated intensity map, Herschel 250 $\mu$m map, and JCMT 850 $\mu$m map are ∼52, 18, and 14 arcsec, respectively, shown as a framed-blue dot at the bottom left of each panel.
Figure 1.

(a) 13CO (J = 1–0) intensity map of G148.24+00.41, integrated in the velocity range −37.0 to −30.0 km s−1. (b) The central region (encompassed by the solid green box in panel-a) of G148.24+00.41as seen in Herschel 250 |$\mu$|m band, showing the hub-filamentary morphology of the cloud. The figure is adapted from Rawat et al. (2023), in which the blue circle shows the JCMT scanned region of diameter ∼12 arcmin (∼12 pc). The green dashed box marks the central area of the hub, where an infrared cluster is seen (Rawat et al. 2023), and the cross sign indicates the position of a massive young stellar object. (c) The 850 |$\mu$|m Stokes I intensity map of the central region of G148.24+00.41 mapped by JCMT SCUBA-2/POL-2, along with the contours of 13CO integrated intensity emission, drawn from 1.5 to 15 K km s−1with a step size of ∼0.96 K km s−1. The rms noise of the 4 arcsec pixel-size Stokes I map is around ∼5 mJy beam−1. In panel-c, the yellow ellipse shows the position of the C1 clump, identified by Rawat et al. (2024) using C18O data (spatial resolution ∼52 arcsec). The beam sizes of the 13CO integrated intensity map, Herschel 250 |$\mu$|m map, and JCMT 850 |$\mu$|m map are ∼52, 18, and 14 arcsec, respectively, shown as a framed-blue dot at the bottom left of each panel.

The magnetohydrodynamic simulations suggest that the filamentary converging flows would impact the magnetic field morphology of the star-forming regions (Gómez, Vázquez-Semadeni & Zamora-Avilés 2018). In G148.24+00.41, Rawat et al. (2024) found converging gas flows along the filaments towards the hub. So, it would be interesting to investigate the influence of gas flows on the B-field morphology in the hub of the cloud. In addition, it is equally important to examine the role of magnetic fields in the stability of such hub systems, as it is expected that gravity would play a dominant role in the onset of star formation within such systems. In this work, we investigate the morphology and strength of the magnetic field of the C1 clump of G148.24+00.41 for the first time, in order to understand the relative importance of the magnetic field in comparison to gravity and turbulence in the overall star formation process of the clump. This paper is organized as follows. In Section 2, we describe the observations, data reduction, and other data sets used in this work. In Section 3, we present the analysis and results related to the magnetic field morphology, its relative orientation in comparison to intensity gradients and local gravity, and its strength. We also discuss the dust and gas properties of the studied region. In Section 4, we discuss the results and compare the strength of gravity, magnetic field, and turbulence. In Section 5, we summarize the findings of this work.

2 OBSERVATIONS AND DATA

2.1 Dust continuum polarization observations using JCMT SCUBA-2/POL-2

We observed the C1 clump/hub region of G148.24+00.41with SCUBA-2/POL-2 mounted on the James Clerk Maxwell Telescope (JCMT), a single-dish sub-millimeter telescope in Mauna Kea, Hawaii, USA. The POL-2 instrument is a linear polarimetry module (Friberg et al. 2016) for the Submillimetre Common User Bolometer Array-2 (SCUBA-2), a 10 000 bolometer camera on the JCMT (Holland et al. 2013). The data was acquired between 2022 November 25 and 2023 January 03 (project code: M22BP055; PI: Vineet Rawat) in the band 2 weather conditions under an atmospheric optical depth at 225 GHz (τ225) of 0.04 to 0.06. The observations were taken in 10 sets with an integration time of 30 minutes each, resulting in a total integration time of around 5.5 h. The POL-2 DAISY scan mode (Holland et al. 2013; Friberg et al. 2016) was adopted, which generates a map of high signal-to-noise ratio (SNR) within a central region spanning a diameter of 3 arcmin, and the noise level gradually increases towards the edges of the map. The region is observed in both the 450 and 850 |$\mu$|m continuum polarizations simultaneously, with a resolution of 9|${_{.}^{\prime\prime}}$|6 and 14|${_{.}^{\prime\prime}}$|1, respectively. Due to the low sensitivity of the 450 |$\mu$|m data, this paper presents the analyses and results based on only 850 |$\mu$|m dust polarization data.

The data reduction was carried out using the pol2map1script in the smurf package (Chapin et al. 2013) of Starlink (Currie et al. 2014). The POL-2 is characterized by linear polarization that produces Stokes I, Q, and U vector maps. The Skyloop mode was utilized to minimize the uncertainty associated with map creation, while the MAPVARS mode was enabled to asses the total uncertainty from the standard deviation among individual observations. The details of the data reduction process of POL-2 can be found in Pattle et al. (2017) and Wang et al. (2019). Finally, the I, Q, and U maps, along with their variance maps, are used to create a debiased polarization vector catalogue. The catalogue consists of total intensity (I), stokes vectors (Q and U), polarization intensity (PI), polarization fraction (P), polarization angle (θP), and their associated uncertainties (δI, δQ, δU, δPI, δP, and δθP, respectively).

The I, Q, and U maps are produced with 4 arcsec pixel size, while the polarization catalogue is binned to 12 arcsec, for better sensitivity. A flux calibration factor of 668.25 Jy beam−1 pW−1 is used for 850 |$\mu$|m Stokes I, Q, and U map to convert them from pW to mJy beam−1 and to account for the flux-loss due to POL-2 insertion into the telescope. This calibration factor comprises 495 Jy/beam/pW for reductions using 4 arcsec pixels of SCUBA-2, multiplied by the standard 1.35 factor for POL-2 losses (Mairs et al. 2021).

The polarised intensity is defined to be positive, so the uncertainties of the Q and U Stokes vector would bias the polarized intensities towards larger values (Vaillancourt 2006; Kwon et al. 2018). The debiased polarization intensity and its uncertainty are calculated as

(1)

and

respectively. The debiased polarization fraction and its uncertainty are then calculated as

(2)

and

respectively. The polarization angle and its uncertainty are calculated as

(3)

and

respectively. The polarization angle increases from the north toward the east, following the IAU convention. The mean rms noises in the Stokes I, Q, U, and PI measurements with 12″ bin size are 1.4, 1.1, 1.1, and 1.1 mJy beam−1, respectively. Following the standard convention, for magnetic field, hereafter, B-field orientations, the polarization angles are rotated by 90 deg.

2.2 Molecular line data from PMO

We used the 13CO (J = 1–0) and C18O (J = 1–0) molecular line data, which were taken as a part of the Milky Way Imaging Scroll Painting survey (MWISP; Su et al. 2019) using the 13.7-m radio telescope at Purple Mount Observatory (PMO). The observation of CO isotopologues was taken simultaneously using a 3 × 3 beam sideband-separating Superconducting Spectroscopic Array Receiver (SSAR) system (Shan et al. 2012) and using the position-switch on-the-fly mode. The spatial resolution (Half Power Beam Width; HPBW) of 13CO and C18O is around ∼52 arcsec (∼0.9 pc at 3.4 kpc). And the spectral resolution of 13CO and C18O is ∼0.17 km s−1with a sensitivity of ∼0.3 K (for details, see Su et al. 2019).

3 ANALYSES AND RESULTS

Fig. 1(c) shows the Stokes I map of the region, where the location of the C1 clump is also shown. From the figure, it can be seen that the 850 |$\mu$|m JCMT data (beam size ∼14 arcsec) has resolved multiple sub-structures in the central region of the cloud. We found sub-structures like a central clump, a clump located on the western side, and a prominent elongated structure on the northeastern side of the central clump. In addition to these prominent sub-structures, a few compact structures are also visible in the image.

3.1 B-field morphology

In order to select the significant polarization detections, we set the following criteria for selecting data: I/δI > 10, P/δP > 2, and P < 30  per cent. By doing this, we got 69 polarization measurements in our target region. The P values range from ∼2 per cent to ∼29 per cent with a mean and standard deviation around ∼11 ± 8 per cent. The B-field orientations are widely distributed, ranging from ∼6to 180with a mean and standard deviation around ∼91± 48, suggesting a complex B-field morphology in the region. The mean uncertainties in polarization fraction and polarization angle are ∼3.5 per cent and ∼9, respectively. Fig. 2(a) and (b) show the distribution of polarization vectors and B-field orientations, respectively, over the 850 |$\mu$|m Stokes I dust continuum emission map of the region. The contour levels in the map are shown above 3σ from the background, where σ is the mean rms noise (5 mJy beam−1) of the Stokes I map.

(a) Polarization vector map with lengths proportional to polarization fraction, and (b) magnetic field orientation map with fixed lengths. The background is the Stokes I image at 850 $\mu$m, and the contour levels are drawn at 3σ above the rms noise level of 5 mJy beam−1, starting from 15 to 300 mJy beam−1. The segments shown are binned to a 12 arcsec pixel grid and correspond to polarization data with I/δI > 10 and P/δP > 2. The light cyan and green vectors in panel-b show the measurements with 2 < P/δP < 3 and P/δP > 3, respectively. The regions used for B-field calculation are also shown in panel-b by a white circle and yellow ellipse for the CC and NES, respectively.
Figure 2.

(a) Polarization vector map with lengths proportional to polarization fraction, and (b) magnetic field orientation map with fixed lengths. The background is the Stokes I image at 850 |$\mu$|m, and the contour levels are drawn at 3σ above the rms noise level of 5 mJy beam−1, starting from 15 to 300 mJy beam−1. The segments shown are binned to a 12 arcsec pixel grid and correspond to polarization data with I/δI > 10 and P/δP > 2. The light cyan and green vectors in panel-b show the measurements with 2 < P/δP < 3 and P/δP > 3, respectively. The regions used for B-field calculation are also shown in panel-b by a white circle and yellow ellipse for the CC and NES, respectively.

In this work, based on 850 |$\mu$|m Stokes I intensity and magnetic field orientations, we defined the central clump, clump located on the western side, and northeastern elongated structure as CC, WC, and NES, respectively, as marked in Fig. 2(b). The approximate extents of these regions are defined by considering the outermost closed contours of the 850 |$\mu$|m Stokes I map. From Fig. 2(b), it can be seen that the B-field orientations in the CC are mostly oriented along the east-west direction (PA ∼90), while some of them are at smaller position angles. There exist mixed B-field orientations in the NES region, some in the low-density area are nearly perpendicular to the major axis of the NES, while some closer to the CC are parallel to it. In the WC, most of the B-fields are converging towards the centre, aligned along the south-east direction, which may be influenced by gravity (see Section 3.3.2 and 4.1). Overall, the B-field morphology around the central region of G148.24+00.41is complex, which is probably due to hierarchical fragmentation and a network of filamentary flows towards the hub, as found by Rawat et al. (2023, 2024).

Fig. 3(a) shows the histogram of the B-field orientations in the central region of G148.24+00.41, which is broadly distributed. The CC is showing mixed morphology, having two peaks, one at ∼38(i.e. with position angles close to north-east), and the second is at ∼80(i.e. with position angles parallel to east). The NES shows a flat distribution over a broad range, but a slightly higher distribution at a position angle around ∼180. The WC, though, has a small number of segments, shows a peak around ∼125, i.e. mostly in the southeast direction. All these orientations are also clearly evident in Fig. 3(b), which shows the distribution of B-field position angles. From Fig. 3(b), it can be seen that the B-field angles change roughly from ∼180to ∼70while going from the elongated structure towards the central clump.

(a) Histogram of B-field position angles for the whole region, CC, NES, and WC. (b) Distribution of B-field position angles over the contours of 850 $\mu$m Stokes I map. The contour levels are the same as in Fig. 2.
Figure 3.

(a) Histogram of B-field position angles for the whole region, CC, NES, and WC. (b) Distribution of B-field position angles over the contours of 850 |$\mu$|m Stokes I map. The contour levels are the same as in Fig. 2.

3.2 Variation of polarization fraction: depolarization effect

Fig. 4 shows the distribution of polarization fraction over the contours of Stokes I emission. From the figure, it can be seen that the polarization fraction is lower in the high-intensity regions compared to the low-intensity regions, which shows the decreasing trend of polarization fraction with the total intensity, known as depolarization. The depolarization effect has been reported in several studies (Girart, Rao & Marrone 2006; Tang et al. 2013; Sadavoy et al. 2018; Soam et al. 2018; Liu et al. 2019, 2020), and is mainly explained by the inefficient radiative alignment of dust grains in high-density regions or integration effect across a complex magnetic fields. In high-density regions, the radiative alignment torques decrease due to the attenuation of interstellar radiation that results in poor grain alignment, and hence the decrease in polarization fraction. However, the grain characteristics like size, shape, composition, and grain growth can also affect the dust grain alignment. The turbulent nature of the B-field and unresolved complex and tangled B-fields within the JCMT beam, being averaged across the beam, can also give low dust polarization (Planck Collaboration XXXV 2016, XII 2020). We want to point out that in the hub/C1 clump of G148.24+00.41, supersonic non-thermal motions have been found (sonic Mach number = 3–4.4, see Rawat et al. 2024). Therefore, the turbulent nature of the B-field can also be the cause of depolarization in our target.

Distribution of dust polarization fraction (P in %) over the contours of 850 $\mu$m Stokes I map. The contour levels are the same as in Fig. 2.
Figure 4.

Distribution of dust polarization fraction (P in %) over the contours of 850 |$\mu$|m Stokes I map. The contour levels are the same as in Fig. 2.

The relation between polarization fraction and intensity is expected to follow a power law, PI−α (Whittet et al. 2008). A range of α values has been found in molecular clouds from ∼0.5 to 1 (Chung et al. 2023, and references therein). The α value is often used as an indicator of the dust grain alignment efficiency. When α = 0, it implies a constant grain alignment efficiency, α = 0.5 implies that the alignment decreases linearly with the increasing optical depth, while α = 1 implies an alignment limited to the outer regions of the cloud, and at higher density, there is no preferred alignment of grains relative to the magnetic field (Whittet et al. 2008). We fit the P–I relation with a single power law (weighted-fit) and found an index, α = 0.95 ± 0.04, which shows that the dust grain alignment efficiency is decreasing in the central dense region of G148.24+00.41. However, Pattle et al. (2019) shows that the conventional approach of fitting a single power law over the polarization measurements debiased with Gaussian noise is only applicable above a high SNR cut. But in low polarized intensity regions, a high SNR would discard more data, and therefore the α index will be overestimated (Pattle et al. 2019; Chung et al. 2023; Lin et al. 2023).Wang et al. (2019) also found that the value of α depends upon the cut of SNR and tends to −1 if we put a constraint on P/δP. Hence, to obtain the true value of the α index, it is recommended to use the non-debiased polarization measurements, including both the low- and high-SNR data, and should not put constraints on P/δP (Pattle et al. 2019; Wang et al. 2019).

We followed the Bayesian method of Wang et al. (2019) to determine the true value of α by using the non-debiased polarization data, which follows well the Rice distribution (see Pattle et al. 2019; Wang et al. 2019, and references therein)

(4)

where P and P0 are the observed and true polarization fraction, respectively, σP is the uncertainty in the polarization fraction, and I0 is the zeroth-order modified Bessel function. We used non-debiased data with an SNR of 2 (i.e. I/δI > 2) to include most of the data points and used the power-law model, P0 = βI−α, with uncertainty σP = σQU/I, where σQU represents the rms noise in Q and U measurements, I is the total observed intensity, and α, β, and σQU are the free model parameters. We employed the Markov Chain Monte Carlo method and used a python package PyMC3 (Salvatier, Wiecki & Fonnesbeck 2016) to fit the Rician model to the data. We set the uniform priors on all three model parameters: 0 < α < 2, 0 < β <100, and 0 < σQU < 5, and otherwise a value of 0 for all the parameters. The details of the methodology are given in Wang et al. (2019). Fig. 5 shows the derived posterior of each model parameter, along with their 95 per cent highest density interval (HDI), depicting the uncertainty in each parameter. The mean values of α, β, and σQU are ∼0.6, 37, and 1.7, respectively. The α value derived from the non-debiased polarization data is smaller than the α value derived from the conventional approach (i.e. ∼0.95). Fig. 6 shows the non-debiased polarization fraction versus total intensity plot with 50, 68, and 95 per cent confidence interval. The derived α value suggests that the grain alignment is still persisting in the hub of G148.24+00.41, but with decreasing efficiency in the dense regions.

The probability distribution function of the fitted model parameters derived using the Bayesian method over the non-debiased polarization data. The mean values of the parameters are shown along with the 95 per cent HDI intervals to represent the uncertainties. The 95 per cent confidence intervals are marked as horizontal bars.
Figure 5.

The probability distribution function of the fitted model parameters derived using the Bayesian method over the non-debiased polarization data. The mean values of the parameters are shown along with the 95 per cent HDI intervals to represent the uncertainties. The 95 per cent confidence intervals are marked as horizontal bars.

Non-debiased polarization fraction versus total intensity. The blue line shows the mean, and the coloured regions show the 95, 68, and 50 per cent confidence limits, as predicted by the posteriors of α = 0.6, β = 37, and σQU = 1.7.
Figure 6.

Non-debiased polarization fraction versus total intensity. The blue line shows the mean, and the coloured regions show the 95, 68, and 50 per cent confidence limits, as predicted by the posteriors of α = 0.6, β = 37, and σQU = 1.7.

3.3 Relative orientations of magnetic fields, intensity gradients, and local gravity

In star-forming regions, various forces interact, like gravity, magnetic field, and turbulence, which shape the geometry of these regions and drive the star-formation process (Ballesteros-Paredes et al. 2007; Koch, Tang & Ho 2012a; Pattle et al. 2022). Along with the overall strength of these individual factors for an entire region (discussed in Section 4.2), it is also important to investigate their localized relative orientations in the map, as it would give insight into the localized effect of these factors (Koch, Tang & Ho 2012a, b, 2013; Tang et al. 2019; Liu et al. 2020; Wang et al. 2020b). Koch, Tang & Ho (2012a) developed a technique, ‘the polarization-intensity gradient-local gravity,’ using Magnetohydrodynamics (MHD) force equations to measure the local magnetic field strengths. Following the approach of Koch, Tang & Ho (2012a, b), we find out the angular difference between magnetic field, intensity gradient, and local gravity, and discuss their relative importance at different positions.

3.3.1 Intensity gradient versus magnetic field

We used the 850 |$\mu$|m dust continuum intensities of all pixels in the map to determine the directions of intensity gradients. All the pixels that have values above a certain threshold (i.e. 3σ above the mean rms noise in the Stokes I map) are considered for computing the direction of gradients, except those that exist at the edges. For a pixel at position (αi, δj), the position angle (θ′IG) of the intensity gradient is calculated as

(5)

where |$\Delta I_{\delta _j} = I_{\delta _{j+1}}-I_{\delta _{j-1}}$| and |$\Delta I_{\alpha _i} = I_{\alpha _{i+1}}-I_{\alpha _{i-1}}$|⁠.

The θ′IG values are then converted to gradient directions (θIG) by doing the quadrant corrections, i.e. arranging the angles between 0 and 360 (for details, see Eswaraiah et al. 2020). In order to plot the gradient orientations instead of directions, we folded the θIG between 0 and 180. For comparison of the gradient orientations (θIG) with the B-field orientations (θB), we took the average of all the θIG values within a diameter of ∼14 arcsec (corresponds to the beam size of JCMT at 850 |$\mu$|m) around each B-field position. We calculated the circular mean to get the average of intensity gradients. In this approach, the angles are treated as unit vectors, which is adequate for broad distributions and ambiguity in angles (Tang et al. 2019). Fig. 7(a) shows the orientations of intensity gradients relative to B-field orientations over the 850 |$\mu$|m Stokes I map.

(a) The orientations of the B-fields (green segments) and intensity gradients (red segments) are overlaid on the 850 $\mu$m Stokes I map. (b) The distribution of the offset between the position angles of the B-fields and intensity gradients, i.e. ΔθB, IG = |(θB − θIG)| over the contours of 850 $\mu$m Stokes I map. The contour levels are same as in Fig. 2.
Figure 7.

(a) The orientations of the B-fields (green segments) and intensity gradients (red segments) are overlaid on the 850 |$\mu$|m Stokes I map. (b) The distribution of the offset between the position angles of the B-fields and intensity gradients, i.e. ΔθB, IG = |(θB − θIG)| over the contours of 850 |$\mu$|m Stokes I map. The contour levels are same as in Fig. 2.

We find that the local differences between these orientations are overall widely distributed. However, it can be seen that the intensity gradients are mostly aligned with the B-fields in the CC and WC regions, while the differences in orientations are relatively higher in the NES region. In the observed central region of G148.24+00.41, we found a moderate correlation between θB and θIG, in the CC and WC (see Figs 7a, b). Fig. 7(b) shows the distribution of ΔθB, IG = |(θB − θIG)| over the contours of 850 |$\mu$|m Stokes I map. The ΔθB, IG values lie between 0 and 90 after considering them to be the acute angle. A stronger correlation between θB and θIG tells that the material is following the B-field lines (Koch, Tang & Ho 2013; Tang et al. 2019, more discussion in Section 4.1).

3.3.2 Local gravitational field versus magnetic field

In order to investigate the localized effect of gravity on the B-field morphology of the structures, we used the 850 |$\mu$|m dust continuum map, to compute the projected gravitational field vectors. The gravitational force at any pixel (FG, i) is the vector sum of the forces from all the surrounding pixels and is expressed as (Wang et al. 2020b).

(6)

where Ii and Ij are the intensity of the pixel at position i and j, respectively, and k is the term that takes care of the conversion of emission to total column density and also includes the gravitational constant. N is the total number of pixels within the selected area, rij is the projected distance between the pixels i and j, and |$\hat{r}$| is the unit vector. Considering only the directions of the local gravitational forces, we take k to be 1 in the above equation by assuming that the spatial distribution of dust will be analogous to the spatial distribution of mass. Similar to the intensity gradient map, we selected those pixels which have intensity values above the threshold, and obtained the local gravity vectors (θLG) at each B-field position, by taking an average of all vectors within the 14 arcsec beam size. Fig. 8(a) shows the orientations of local gravity vectors relative to B-field orientations over the 850 |$\mu$|m Stokes I map. From the figure, it can be seen that similar to intensity gradients, the local gravity vectors are also mostly aligned with the B-fields in the CC and WC region, whereas they deviate from the B-fields in the NES region. Fig. 8(b) shows the distribution of ΔθB, LG = |(θB − θLG)| values over the contours of 850 |$\mu$|m Stokes I emission, which are treated to be acute angles.

(a) The orientations of the B-fields (green segments) and local gravity (magenta vectors) are overlaid on the 850 $\mu$m Stokes I map. (b) The distribution of the offset between the position angles of the B-fields and local gravity, i.e. ΔθB, LG = |(θB − θLG)| over the 850 $\mu$m Stokes I map. The contour levels are same as in Fig. 2.
Figure 8.

(a) The orientations of the B-fields (green segments) and local gravity (magenta vectors) are overlaid on the 850 |$\mu$|m Stokes I map. (b) The distribution of the offset between the position angles of the B-fields and local gravity, i.e. ΔθB, LG = |(θB − θLG)| over the 850 |$\mu$|m Stokes I map. The contour levels are same as in Fig. 2.

The relative orientations of intensity gradient and local gravity have a similar distribution like that of ΔθB, IG and ΔθB, LG, as shown in Fig. A1 in appendix A. Figs 9(a)–(c) show the histogram distribution of the relative position angle differences of B-field, intensity gradients, and local gravity, i.e. ΔθB, IG = |(θB − θIG)|, ΔθB, LG = |(θB − θLG)|, and ΔθIG, LG = |(θIG − θLG)|. From the figure, it can be seen that the differences in the offset angles are mostly distributed towards the smaller angles. The median of ΔθB, IG, ΔθB, LG, and ΔθIG, LG is 34, 32, and 30 with median absolute deviation of 22, 18, and 15, respectively. The higher angular deviations in the histograms are primarily due to position angles in the elongated structures on the northeastern side, as well as from some structures located south of the CC.

Distribution of difference in position angles of (a) magnetic field (θB) and intensity gradient (θIG), (b) magnetic field (θB), and local gravity (θLG), and (c) intensity gradient (θIG) and local gravity (θLG). The red, blue, and green histograms show the difference in position angles within the CC, NES, and WC regions, respectively.
Figure 9.

Distribution of difference in position angles of (a) magnetic field (θB) and intensity gradient (θIG), (b) magnetic field (θB), and local gravity (θLG), and (c) intensity gradient (θIG) and local gravity (θLG). The red, blue, and green histograms show the difference in position angles within the CC, NES, and WC regions, respectively.

As discussed previously, Koch, Tang & Ho (2012a, b) developed a ‘polarization-intensity gradient method’ that can estimate the local field-to-gravity force ratio ΣB. This method is based on the assumption that the emission intensity gradients reflect the direction of matter flow due to the combined influences of magnetic pressure force and gravitational force. Using the MHD force equations and geometrically solving them by incorporating the angle between the magnetic field and intensity gradient (ΔB, IG), and between intensity gradient and local gravity (ΔIG, LG), the magnetic field (FB)-to-gravity force (FG) ratio can be obtained as

(7)

In the above equation, the hydrostatic gas pressure is assumed to be negligible. Fig. 10 shows the ΣB distribution plot over the Stokes 850 |$\mu$|m intensity map. From the figure, it can be seen that ΣB is mostly ≤1, with a median around ∼0.6, which shows that the magnetic field is not solely enough to balance the gravitational force (Koch, Tang & Ho 2012a). This implies that gravity dominates over the magnetic field to govern the gas motion towards the centre. However, we note that the POL-2 images generally filter out large-scale structures, and so here, the intensity gradient only traces the local structure on a 4 arcsec pixel-scale. Therefore, to check the effect of large-scale structures on the intensity gradient and local gravity, we generated similar maps from Herschel 250 |$\mu$|m image of G148.24+00.41and found that the maps are comparable with the JCMT maps. The use of 250 |$\mu$|m map is an optimal choice because compared to Herschel′s longer wavelength (i.e. 350 and 500 |$\mu$|m) bands, its resolution (∼18 arcsec) is comparable to the resolution of the JCMT 850 |$\mu$|m (14 arcsec) map and also it is a better tracer of cold dust compared to Herschel′s shorter wavelength (i.e. 70 and 160 |$\mu$|m) bands.

ΣB distribution over the contours of 850 $\mu$m Stokes I map. The contour levels are same as in Fig. 2.
Figure 10.

ΣB distribution over the contours of 850 |$\mu$|m Stokes I map. The contour levels are same as in Fig. 2.

Due to the relatively low number of B-field segments in the WC region, we focused our further analysis, like the study of structure, dust properties, and B-field strength calculation, towards the CC and NES regions.

3.4 Column and number densities

We calculate the molecular hydrogen column density using the 850 |$\mu$|m dust continuum emission. Assuming the dust emission to be optically thin, the column density can be calculated using the relation (Kauffmann et al. 2008),

(8)

where Sν is the flux density in Jy at frequency ν, λ is the wavelength (0.85 mm), Td is the mean dust temperature in K, dust opacity κν = 0.1(ν/1 THz)β = 0.0125 cm|$\rm {^2 g^{-1}}$| for ν = 0.353 THz and dust opacity index, β = 2 (Battersby et al. 2011; Deharveng et al. 2012), and θHPBW is the beam size (14 arcsec at 850 |$\mu$|m).

The mean Td is taken from the Herschel Hi-GAL survey (Molinari et al. 2010) dust temperature map (Schisano et al. 2020) of the G148.24+00.41 region. Within the boundaries of CC and NES, the mean Td is around ∼16.8 and 13.0 K, respectively. The total column density, ∑N(H2), for CC and NES, are found to be (3.4 ± 1.3) × 1024 cm−2 and (1.9 ± 0.7) × 1024 cm−2, respectively. Then, assuming the spherical geometry for CC, its number density (⁠|$n_{H_2}$|⁠) can be estimated using the relation

(9)

where M and V are the mass and volume of the region, respectively, μ is the mean molecular weight assumed to be 2.8 (Kauffmann et al. 2008), mH is the mass of hydrogen, Apixel is the pixel area, and Reff is the effective radius. The Reff for CC is calculated as (Area/π)0.5, and is around ∼0.9 ± 0.1 pc. Though NES is assumed as an elliptical structure with a semi-minor axis, r1 = 0.35 ± 0.03 pc and a semi-major axis, r2 = 1.10 ± 0.09 pc (see Fig. 2a), in 3-dimension, this elongated structure could be better described by a cylindrical geometry. Therefore, to calculate the number density of NES, we adopted its radius (r) and length (L) to be r1 and 2r2, respectively. Under this approximation, we estimated the |$n_{H_2}$| for NES by using V = πr2L in equation (9). The gas mass within the regions is estimated from the total integrated molecular hydrogen column density, using equation (1) given in Rawat et al. (2023). The total column density, |$n_{H_2}$|⁠, and the mass of the regions are given in Table 1. The uncertainties in the estimated cloud parameters are mainly due to uncertainty in the gas-to-dust ratio (23 per cent), the dust opacity index (30 per cent), and the distance of the cloud (9 per cent) (Rawat et al. 2023, and references therein).

Table 1.

Parameters estimated for CC and NES.

NoParameterUnitCCNES
1Effective radius (Reff)pc0.9 ± 0.1semiminor axis (r1) = 0.35 ± 0.03, semimajor axis (r2) = 1.10 ± 0.09
2Mean column density (N(H2))cm−2(5.8 ± 2.2)× 1021(7.1 ± 2.7) × 1021
3Number density (⁠|$n_{H_2}$|⁠)cm−31560 ± 7803290 ± 1645
4Mass (M)M330 ± 148188 ± 85
5Mean dust temperature (Td)K16.813.0
6Observed velocity dispersion (σobs)km s−10.69 ± 0.050.41 ± 0.03
7Thermal velocity dispersion (σth)km s−10.060.05
8Non-thermal velocity dispersion (σnt)km s−10.69 ± 0.050.41 ± 0.03
Structure function analysis
1Turbulent-to-ordered magnetic field energy ratio |$\left(\frac{\mathinner {\langle {\delta B^2}\rangle }^{1/2}}{B_0}\right)$|0.43 ± 0.010.46 ± 0.04
2Angular dispersion (σθ)deg22.7 ± 0.623.9 ± 1.5
3Plane-of-sky magnetic field strength (Bpos)μG24.0 ± 6.020.0 ± 5.0
4Mass-to-flux ratio (λB)1.8 ± 0.82.7 ± 1.2
5Alfv|$\acute{\text{e}}$|n velocity (VA)km s−11.0 ± 0.40.6 ± 0.2
6Alfv|$\acute{\text{e}}$|n mach number (⁠|$\mathcal {M}_A$|⁠)1.2 ± 0.41.2 ± 0.4
7Magnetic pressure (⁠|$P\rm {_B}$|⁠)|$\rm {dyne\, cm^{-2}}$|(3.7 ± 1.9) × 10−11(2.6 ± 1.3) × 10−11
8Turbulent pressure (⁠|$P\rm {_{turb}}$|⁠)|$\rm {dyne\, cm^{-2}}$|(5.2 ± 2.7) × 10−11(2.5 ± 1.3) × 10−11
Virial balance
1Kinetic energy (EK)J(4.7 ± 2.2) × 1038(6.3 ± 3.0) × 1037
2Magnetic energy (EB)J(3.3 ± 3.0) × 1038(7.0 ± 5.0) × 1037
3Gravitational energy (EG)J(10.4 ± 9.0) × 1038(14 ± 12) × 1037
4Kinetic virial parameter (αvir, k)0.9 ± 0.40.9 ± 0.4
5Total virial parameter (αvir, tot)1.2 ± 0.61.4 ± 0.7
NoParameterUnitCCNES
1Effective radius (Reff)pc0.9 ± 0.1semiminor axis (r1) = 0.35 ± 0.03, semimajor axis (r2) = 1.10 ± 0.09
2Mean column density (N(H2))cm−2(5.8 ± 2.2)× 1021(7.1 ± 2.7) × 1021
3Number density (⁠|$n_{H_2}$|⁠)cm−31560 ± 7803290 ± 1645
4Mass (M)M330 ± 148188 ± 85
5Mean dust temperature (Td)K16.813.0
6Observed velocity dispersion (σobs)km s−10.69 ± 0.050.41 ± 0.03
7Thermal velocity dispersion (σth)km s−10.060.05
8Non-thermal velocity dispersion (σnt)km s−10.69 ± 0.050.41 ± 0.03
Structure function analysis
1Turbulent-to-ordered magnetic field energy ratio |$\left(\frac{\mathinner {\langle {\delta B^2}\rangle }^{1/2}}{B_0}\right)$|0.43 ± 0.010.46 ± 0.04
2Angular dispersion (σθ)deg22.7 ± 0.623.9 ± 1.5
3Plane-of-sky magnetic field strength (Bpos)μG24.0 ± 6.020.0 ± 5.0
4Mass-to-flux ratio (λB)1.8 ± 0.82.7 ± 1.2
5Alfv|$\acute{\text{e}}$|n velocity (VA)km s−11.0 ± 0.40.6 ± 0.2
6Alfv|$\acute{\text{e}}$|n mach number (⁠|$\mathcal {M}_A$|⁠)1.2 ± 0.41.2 ± 0.4
7Magnetic pressure (⁠|$P\rm {_B}$|⁠)|$\rm {dyne\, cm^{-2}}$|(3.7 ± 1.9) × 10−11(2.6 ± 1.3) × 10−11
8Turbulent pressure (⁠|$P\rm {_{turb}}$|⁠)|$\rm {dyne\, cm^{-2}}$|(5.2 ± 2.7) × 10−11(2.5 ± 1.3) × 10−11
Virial balance
1Kinetic energy (EK)J(4.7 ± 2.2) × 1038(6.3 ± 3.0) × 1037
2Magnetic energy (EB)J(3.3 ± 3.0) × 1038(7.0 ± 5.0) × 1037
3Gravitational energy (EG)J(10.4 ± 9.0) × 1038(14 ± 12) × 1037
4Kinetic virial parameter (αvir, k)0.9 ± 0.40.9 ± 0.4
5Total virial parameter (αvir, tot)1.2 ± 0.61.4 ± 0.7
Table 1.

Parameters estimated for CC and NES.

NoParameterUnitCCNES
1Effective radius (Reff)pc0.9 ± 0.1semiminor axis (r1) = 0.35 ± 0.03, semimajor axis (r2) = 1.10 ± 0.09
2Mean column density (N(H2))cm−2(5.8 ± 2.2)× 1021(7.1 ± 2.7) × 1021
3Number density (⁠|$n_{H_2}$|⁠)cm−31560 ± 7803290 ± 1645
4Mass (M)M330 ± 148188 ± 85
5Mean dust temperature (Td)K16.813.0
6Observed velocity dispersion (σobs)km s−10.69 ± 0.050.41 ± 0.03
7Thermal velocity dispersion (σth)km s−10.060.05
8Non-thermal velocity dispersion (σnt)km s−10.69 ± 0.050.41 ± 0.03
Structure function analysis
1Turbulent-to-ordered magnetic field energy ratio |$\left(\frac{\mathinner {\langle {\delta B^2}\rangle }^{1/2}}{B_0}\right)$|0.43 ± 0.010.46 ± 0.04
2Angular dispersion (σθ)deg22.7 ± 0.623.9 ± 1.5
3Plane-of-sky magnetic field strength (Bpos)μG24.0 ± 6.020.0 ± 5.0
4Mass-to-flux ratio (λB)1.8 ± 0.82.7 ± 1.2
5Alfv|$\acute{\text{e}}$|n velocity (VA)km s−11.0 ± 0.40.6 ± 0.2
6Alfv|$\acute{\text{e}}$|n mach number (⁠|$\mathcal {M}_A$|⁠)1.2 ± 0.41.2 ± 0.4
7Magnetic pressure (⁠|$P\rm {_B}$|⁠)|$\rm {dyne\, cm^{-2}}$|(3.7 ± 1.9) × 10−11(2.6 ± 1.3) × 10−11
8Turbulent pressure (⁠|$P\rm {_{turb}}$|⁠)|$\rm {dyne\, cm^{-2}}$|(5.2 ± 2.7) × 10−11(2.5 ± 1.3) × 10−11
Virial balance
1Kinetic energy (EK)J(4.7 ± 2.2) × 1038(6.3 ± 3.0) × 1037
2Magnetic energy (EB)J(3.3 ± 3.0) × 1038(7.0 ± 5.0) × 1037
3Gravitational energy (EG)J(10.4 ± 9.0) × 1038(14 ± 12) × 1037
4Kinetic virial parameter (αvir, k)0.9 ± 0.40.9 ± 0.4
5Total virial parameter (αvir, tot)1.2 ± 0.61.4 ± 0.7
NoParameterUnitCCNES
1Effective radius (Reff)pc0.9 ± 0.1semiminor axis (r1) = 0.35 ± 0.03, semimajor axis (r2) = 1.10 ± 0.09
2Mean column density (N(H2))cm−2(5.8 ± 2.2)× 1021(7.1 ± 2.7) × 1021
3Number density (⁠|$n_{H_2}$|⁠)cm−31560 ± 7803290 ± 1645
4Mass (M)M330 ± 148188 ± 85
5Mean dust temperature (Td)K16.813.0
6Observed velocity dispersion (σobs)km s−10.69 ± 0.050.41 ± 0.03
7Thermal velocity dispersion (σth)km s−10.060.05
8Non-thermal velocity dispersion (σnt)km s−10.69 ± 0.050.41 ± 0.03
Structure function analysis
1Turbulent-to-ordered magnetic field energy ratio |$\left(\frac{\mathinner {\langle {\delta B^2}\rangle }^{1/2}}{B_0}\right)$|0.43 ± 0.010.46 ± 0.04
2Angular dispersion (σθ)deg22.7 ± 0.623.9 ± 1.5
3Plane-of-sky magnetic field strength (Bpos)μG24.0 ± 6.020.0 ± 5.0
4Mass-to-flux ratio (λB)1.8 ± 0.82.7 ± 1.2
5Alfv|$\acute{\text{e}}$|n velocity (VA)km s−11.0 ± 0.40.6 ± 0.2
6Alfv|$\acute{\text{e}}$|n mach number (⁠|$\mathcal {M}_A$|⁠)1.2 ± 0.41.2 ± 0.4
7Magnetic pressure (⁠|$P\rm {_B}$|⁠)|$\rm {dyne\, cm^{-2}}$|(3.7 ± 1.9) × 10−11(2.6 ± 1.3) × 10−11
8Turbulent pressure (⁠|$P\rm {_{turb}}$|⁠)|$\rm {dyne\, cm^{-2}}$|(5.2 ± 2.7) × 10−11(2.5 ± 1.3) × 10−11
Virial balance
1Kinetic energy (EK)J(4.7 ± 2.2) × 1038(6.3 ± 3.0) × 1037
2Magnetic energy (EB)J(3.3 ± 3.0) × 1038(7.0 ± 5.0) × 1037
3Gravitational energy (EG)J(10.4 ± 9.0) × 1038(14 ± 12) × 1037
4Kinetic virial parameter (αvir, k)0.9 ± 0.40.9 ± 0.4
5Total virial parameter (αvir, tot)1.2 ± 0.61.4 ± 0.7

3.5 Velocity dispersion

We used C18O molecular line data to determine the non-thermal velocity dispersion along the line-of-sight. Rawat et al. (2024) found that in comparison to 12CO and 13CO, C18O emission is optically thin in the G148.24+00.41 cloud. Fig. 11 shows the C18O spectra averaged within the boundary of the two regions, CC and NES. The Gaussian fitting over the spectra gives the mean velocity as −34.15 ± 0.04 and −33.80 ± 0.04 km s−1, and velocity dispersion (σobs) as 0.69 ± 0.05 and 0.41 ± 0.03 km s−1, for CC and NES, respectively.

C18O average spectral profile for (a) central clump and (b) northeastern elongated structure.
Figure 11.

C18O average spectral profile for (a) central clump and (b) northeastern elongated structure.

The non-thermal velocity dispersion (σnt) can be calculated using the relation, |$\sigma _{\mathrm{ nt}} = \sqrt{\sigma _{\mathrm{ obs}}^2-\sigma _{\mathrm{ th}}^2}$|⁠, where σth is the thermal velocity dispersion, defined as |$\sqrt{\frac{K_\mathrm{ B} T_\mathrm{ k}}{M_{\hbox{C$^{18}$O} }}}$|⁠. Here, KB is the Boltzmann constant, Tk is the gas kinetic temperature, and |$M_{\hbox{C$^{18}$O} }$| is the mass of C18O molecule (30 a.m.u). We have used the excitation temperature map from Rawat et al. (2024) to get the approximate values of Tk, as 11.4 and 10.3 K for CC and NES, respectively. The estimated σth for the regions is ∼ 0.06 and 0.05 km s−1, respectively, and hence negligible, leading σnt ∼ σobs. This is an indication of the presence of turbulence in CC and NES, which has also been reported in Rawat et al. (2024) for the C1 clump.

3.6 Magnetic field strength

Davis (1951) and Chandrasekhar & Fermi (1953) proposed a method to estimate the plane-of-sky component of the magnetic field (Bpos), known as the Davis–Chandrasekhar–Fermi (DCF) method, which is based on the assumption that the turbulence-induced Alfv|$\acute{\text{e}}$|n waves perturb the ordered B-field structure. Therefore, there will be a distorted component of the B-field that would appear as an irregular scatter in polarization angles in comparison to those that are produced by large-scale ordered B-field. Thus, the DCF method implies that the ratio of turbulent (δB) to ordered B-field (B0) is proportional to the ratio of non-thermal velocity dispersion (σnt) to Alfv|$\acute{\text{e}}$|n velocity (⁠|$V_A = B_0/\sqrt{4 \pi \rho }$|⁠, ρ is the gas mass density), i.e. |$\frac{\delta B}{B_0} = \frac{\sigma _{nt}}{V_A}$|⁠. Also, the dispersion in the B-field position angles (σθ) about the large-scale ordered B-field is assumed as |$\sigma _\theta = \frac{\delta B}{B_0}$|⁠. Using these relations, the plane-of-sky magnetic field component, Bpos can be estimated as

(10)

where Q is the correction factor for the line-of-sight and beam-integration effects (Ostriker, Stone & Gammie 2001). The studies show that the beam-integration effect can lead to an underestimation of angular dispersion in polarization angles, resulting in an overestimation of the magnetic field strength (Ostriker, Stone & Gammie 2001; Padoan et al. 2001; Houde et al. 2009). To determine the angular dispersion, there are different statistical methods (see Hildebrand et al. 2009; Houde et al. 2009; Pattle et al. 2017; Liu, Zhang & Qiu 2022), which are used in the literature. In this work, we use the structure-function (SF; Hildebrand et al. 2009) method, which gives the |$\frac{\delta B}{B_0}$| ratio by accounting for the spatial variation of position angles.

3.6.1 Structure function analysis

In the SF method, the magnetic field is assumed to be composed of a large-scale structured field and a turbulent field that are statistically independent. The distinctive behaviour of the two components enables them to distinguish and extract the turbulent component, facilitating the computation of σθ. The SF method computes the difference in position angles, |$\Delta \phi (l) \equiv \phi (\bf x)-\phi (x + {\bf {\it l}})$|⁠, between the N(l) pairs of pixels separated by l = |l|, using the following function:

(11)

This function is referred to as the ‘angular dispersion function’. We want to point out that the polarization position angles are used here for the dispersion function. The angular dispersions, Δϕ, are kept ≤90, to avoid the effect of the ±180 ambiguity of the magnetic field lines. Under the limit, δ < l < < d, the square of the angular dispersion function, known as the ‘structure function’ is characterized by (Hildebrand et al. 2009)

(12)

where δ is the correlation length of the turbulent component, and d is the typical length for variation in large-scale B-field. The quadratically added terms in the dispersion function, m2l2 and b2, are the contribution from the B0 and δB, respectively. The B0 is expected to increase almost linearly with slope m for l < < d, and b is a constant turbulent contribution for l > δ (for details, see Hildebrand et al. 2009). The |$\sigma ^2_M(l)$| is the contribution from the measured uncertainty in the position angles.

The turbulent to large-scale magnetic field strength is given by (Hildebrand et al. 2009)

(13)

and B0 can be estimated by using the modified DCF relation:

(14)

Using the Q correction factor, we can determine the plane-of-sky magnetic field strength:

(15)

where Q is taken to be 0.5 (Heitsch et al. 2001; Ostriker, Stone & Gammie 2001).

We calculated the dispersion function corrected by measurement uncertainty, i.e. |$\mathinner {\langle {\Delta \phi ^2(l)}\rangle }_{tot}-\sigma ^2_M(l)$| with a bin size of 12 arcsec. We used various bin sizes and found that the fit is converged, and fitting errors are the least for 12 arcsec. Fig. 12 shows the dispersion in position angles as a function of the length scale for CC and NES. We fitted the dispersion function with the model defined in equation (12) using least square fit, over the first few data points to ensure the limit, l < < d. The best-fits turbulent component, b, for CC and NES are 32.1 ± 0.9 and 33.8 ± 2.1, respectively. The dispersion in position angles can be obtained as |$\sigma _\theta = b/\sqrt{2}$|⁠, which is around ∼22.7 ± 0.6 and 23.9 ± 1.5 for CC and NES, respectively. These values are close to the maximum value at which the DCF methods would give reliable results (σθ ≤ 25, Ostriker, Stone & Gammie 2001) if a correction factor of 0.5 is applied. Using equation (13), (14), and (15), we calculate the |$\frac{\delta B}{B_0}$| ratio to be around ∼0.43 ± 0.01 and 0.46 ± 0.04, and Bpos to be around ∼24.0 ± 6.0 and 20.0 ± 5.0 |$\mu$|G, for CC and NES, respectively. All the estimated parameters are given in Table 1.

The angular dispersion function for (a) CC and (b) NES. The solid curve represents the best-fitting model to the data, and the points used for fitting are shown in black encircled circles. The intercept of the best-fitting model (l = 0) gives the turbulent contribution to the total angular dispersion. The error bars denote the statistical uncertainties after binning and propagating the individual measurement uncertainties.
Figure 12.

The angular dispersion function for (a) CC and (b) NES. The solid curve represents the best-fitting model to the data, and the points used for fitting are shown in black encircled circles. The intercept of the best-fitting model (l = 0) gives the turbulent contribution to the total angular dispersion. The error bars denote the statistical uncertainties after binning and propagating the individual measurement uncertainties.

4 DISCUSSION

The plane-of-sky magnetic field strength for the CC and NES regions is around ∼24 and 20 |$\mu$|G, respectively. Given the uncertainty of at least a factor of 2 associated with the B-field estimation by the DCF method (Crutcher 2012), the estimated B-field strengths are within the range of ∼10–100 |$\mu$|G observed in star-forming regions (Chapman et al. 2011; Crutcher 2012; Pattle et al. 2022). The values are also consistent (i.e. within a factor of 1.5) with the upper limits of the B-field values from Crutcher et al. (2010) relation for the respective density of the regions.

4.1 Correlation between magnetic fields, intensity gradients, and local gravity

Due to the low statistics in the polarization data, it is difficult to conclusively comment on the overall morphology of the B-field, intensity gradients, and local gravity. However, through this comparison, some inferences can be drawn over their relative spatial variance. Figs 7 and 8 shows a correlation between the B-field, intensity gradients, and local gravity in the CC and WC regions, whereas the differences in their position angles are relatively higher in the NES region. This correlation can be due to the collapse of the clumps where gravity has pulled in and aligned B-field lines with the intensity gradients, in the CC and WC regions (Koch, Tang & Ho 2013; Wang et al. 2019). However, in the outer diffused regions, like the elongated structures, the field lines are not yet that much affected by the local gravity.

Some simulations of the global collapse of magnetized clouds have found that the gravitational flows from large scale to small scale, i.e. from filaments to clumps/cores, can drag the magnetic field lines along the flow, causing a ‘U’ shaped geometry of the field lines across the filament spine (Gómez, Vázquez-Semadeni & Zamora-Avilés 2018; Vázquez-Semadeni et al. 2019). We also found evidence of ‘U’ shaped geometry of B-field lines, at the bottom of the elongated part of the central clump. Fig. 13 shows the zoomed-in view of the central clump region, in which the ‘U’ shaped geometry is sketched from the observed B-field orientations and is shown by magenta curves. High-resolution and sensitivity observations would be required to ascertain the observed morphology. A similar effect of gravity over the B-field morphology has also been found in other observational works (Tang et al. 2019; Beuther et al. 2020; Busquet 2020; Pillai et al. 2020; Wang et al. 2020a, b).

Stokes I map of the CC region of G148.24+00.41, over which the observed B-field segments are shown (green segments). The B-field segments show a ‘U’ shaped geometry, sketched with observed segments and shown by magenta dashed curves, which may be caused by the drag of gravitational converging flows towards the centre of the cloud.
Figure 13.

Stokes I map of the CC region of G148.24+00.41, over which the observed B-field segments are shown (green segments). The B-field segments show a ‘U’ shaped geometry, sketched with observed segments and shown by magenta dashed curves, which may be caused by the drag of gravitational converging flows towards the centre of the cloud.

In G148.24+00.41, Rawat et al. (2024), using MWISP CO data, identified filamentary gas flows exhibiting noticeable velocity gradients as they move towards the hub. They found that the velocity gradient increases towards the hub, with a measured value of ∼0.2 km s−1 pc−1 in the proximity of the hub. The aforementioned findings give evidence of accreting gas flows along the filaments towards the cloud’s centre, which can affect the B-field morphology by dragging them along the flow. Future high-resolution dust continuum and polarimetric observations might be able to reveal a significant number of polarization vectors at the clump/core scale to better resolve the substructures and the B-field morphology around the hub of G148.24+00.41.

4.2 Gravitational stability

4.2.1 Mass-to-flux ratio

In order to investigate whether the magnetic field can provide stability to the regions against gravitational collapse, we determine the mass-to-flux ratio (Crutcher et al. 2004). It is generally calculated as a dimensionless critical stability parameter, λB, which is basically the comparison of the mass to magnetic flux ratio with the critical ratio (Crutcher et al. 2004):

(16)

where N(H2) is the mean column density of the region. A clump or dense core is magnetically supercritical (λB > 1) if the magnetic field is not strong enough to support the system against gravitational collapse, whereas a strong magnetic field would make the system magnetically subcritical (λB < 1), i.e. stable against collapse. Using the mean N(H2) ∼(5.8 ± 2.2) × 1021 and (7.1 ± 2.7) × 1021 cm−2, we calculated the critical parameter to be around ∼1.8 ± 0.8 and 2.7 ± 1.2 for CC and NES, respectively. The critical parameter shows that both regions are magnetically supercritical. In general, a statistical correction factor is applied to λB to account for bias due to geometric effects (Crutcher et al. 2004). Different correction factors are suggested for different geometry of clumps with respect to magnetic field, e.g. π/4 for spherical clump and 1/3 for oblate spheroid with major-axis perpendicular to the mean B-field (Crutcher et al. 2004), and 3/4 for prolate spheroid with major-axis parallel to the mean B-field (Planck Collaboration XXXV 2016). Using these correction factors, the mass-to-flux ratios become subcritical to transcritical/supercritical in the range of 0.6 to 1.4 for CC with mean ∼ 1.1, and 0.9 to 2.1 for NES with mean ∼ 1.7. However, if we consider the overestimation of Bpos in the DCF method itself, our estimated mass-to-flux ratios could also be a lower limit. We want to acknowledge here that the estimated B-field strengths and mass-to-flux ratios are only averaged values over the selected regions. The regions can still be supercritical inside but sub-critical in the outer part, as also shown in a recent simulation by Gómez, Vázquez-Semadeni & Palau (2021).

4.2.2 Turbulence versus magnetic field

Simulations suggest that turbulence plays a dual role in the clouds and their substructures by providing turbulent support against the gravitational collapse at a large scale while producing compressions and shocks at small scales, that create density enhancements and trigger the star formation process (Mac Low & Klessen 2004; Ballesteros-Paredes et al. 2007; Hennebelle & Falgarone 2012; Klessen & Glover 2016). Also, whether it is turbulence (weak B-field models) or magnetic field (strong B-field models) or both, is a subject of investigation regarding their respective roles in the formation of clumps, cores, and subsequent star formation within these structures. In order to investigate their impact, we need to find out the relative strength of turbulence in comparison to the magnetic field.

The Alfv|$\acute{\text{e}}$|n Mach number (⁠|$\mathcal {M}_\mathrm{ A}$|⁠) infers the relative importance of turbulence and magnetic field in molecular clouds and is defined as |$\mathcal {M}_\mathrm{ A}$|= |$\sqrt{3} \sigma _{\mathrm{ nt}}/V_\mathrm{ A}$|⁠. The Alfv|$\acute{\text{e}}$|n velocity is calculated as, |$V_\mathrm{ A} = B_{\mathrm{ tot}}/\sqrt{4\pi \rho }$|⁠. The total B-field strength (Btot) of the regions can be determined by using a statistical relation, Btot = (4/π)Bpos (Crutcher et al. 2004). For CC and NES, the Btot is found to be around 31 ± 8 and 26 ± 6 |$\mu$|G, respectively. Using the mass density estimated within the dimensions of two regions (given in Table 1) and corresponding Btot values, the VA for CC and NES are found to be ∼1.0 ± 0.4 and 0.6 ± 0.2 km s−1, respectively. Then, using the corresponding σnt values, the |$\hbox{$\mathcal {M}_\mathrm{ A}$}$| is calculated to be around ∼1.2 ± 0.4 for both the two regions. A star-forming region is super-Alfv|$\acute{\text{e}}$|nic if |$\mathcal {M}_\mathrm{ A}$| > 1, which means that the turbulent pressure is higher than the magnetic pressure. Conversely, it will be sub-Alfv|$\acute{\text{e}}$|nic if |$\mathcal {M}_\mathrm{ A}$| < 1, which means that the magnetic pressure is higher than the turbulent pressure. In the present case, both the regions are trans Alfv|$\acute{\text{e}}$|nic.

We also calculated the magnetic and turbulent pressures using the relations:

(17)

here, a factor of 3/2 is included to estimate the total turbulent pressure by assuming the non-thermal velocity dispersion to be isotropic in spherical geometry. For the CC region, the magnetic and turbulent pressure is found to be ∼(3.7 ± 1.9) × 10−11 and ∼(5.2 ± 2.7) × 10−11 dyne |$\rm {cm^{-2}}$|⁠, respectively. For the NES region, PB and Pt are ∼(2.6 ± 1.3) × 10−11 and (2.5 ± 1.3) × 10−11 dyne |$\rm {cm^{-2}}$|⁠, respectively. In the CC region, the turbulent pressure is higher than the magnetic pressure by a factor of ∼1.4, while in NES, the pressure values are similar. Also, the thermal pressure (⁠|$\sim \rho \sigma _{\mathrm{ th}}^2$|⁠) in both regions, is much smaller than the turbulent and magnetic pressure, which shows that the thermal energy plays a negligible role in energy balance. So, from the Alfv|$\acute{\text{e}}$|n Mach number and pressure estimation, it seems that the turbulence is slightly more dominant in comparison to the magnetic field in CC, i.e. in the centremost part of G148.24+00.41, while it is similar in the NES.

4.2.3 Virial analysis

The virial theorem is a principle that relates the average kinetic energy (EK) and magnetic field energy (EB) of a system to its average gravitational potential energy (EG), which provides insights into the stability and energy distribution of the system. The Virial theorem is written as:

(18)

where |$\mathcal {I}$| is the moment of inertia. The surface energy terms here are neglected. The kinetic energy term is given by (Fiege & Pudritz 2000)

(19)

where M is the mass and |$\sigma _{obs}^2 = \sigma _{th}^2 + \sigma _{nt}^2$|⁠. The magnetic field energy is given by

(20)

and the gravitational potential energy is given by

(21)

where G is the gravitational constant, R, is the effective radius of the sphere, and L is the length of the cylinder. Here a is the density profile index of the sphere (ρ ∝ ra). The energy values for the two regions are listed in Table 1. We found that for CC, |EG| > EK > EB (i.e. 10.4:4.7:3.3), while for NES, |EG| > EKEB (i.e. 14:6.3:7.0), which restates the results of Alfv|$\acute{\text{e}}$|n Mach number and mass-to-flux ratio calculation, i.e. turbulence is slightly more dominant than the magnetic field in CC, while it is similar in NES, but overall the gravity is the dominant factor in both the regions. Nevertheless, the calculation of individual energy terms is an important exercise, enabling the direct comparison of various forces that govern the evolution of clumps/structures, expressed in the same units.

For a non-magnetized system (EB = 0), the stability criteria is given by 2EK + EG < 0, and based on this condition, the kinetic virial parameter is defined as (Bertoldi & McKee 1992)

(22)

Using the M, R, L, and σobs values of the regions (given in Table 1) in equation (22), and adopting a = 2 for spherical case, we calculate the αvir, k to be around 0.9 ± 0.4 and 0.9 ± 0.4, respectively for CC and NES. The derived αvir, k < 2 means that in the case of a non-magnetized sphere (EB = 0), the thermal plus turbulent contribution is not enough to provide stability to the regions against the gravitational collapse (Kauffmann, Pillai & Goldsmith 2013; Mao, Ostriker & Kim 2020). Including the magnetic energy term in the stability criteria, i.e. 2EK + EB + EG < 0, the total virial parameter is calculated using the modified relation (Bertoldi & McKee 1992; Pillai et al. 2011; Sanhueza et al. 2017).

(23)

With the magnetic support, the αvir, tot value is estimated to be around ∼1.2 ± 0.6 and 1.4 ± 0.7 for CC and NES, respectively. The αvir, tot values for the regions are < 2, which shows that the two regions are bound by gravity and thus can collapse to form stars. The virial analysis shows that the total kinetic energy (EK, i.e. thermal plus turbulent) in both regions is not sufficient to support them against the gravitational collapse. While magnetic energy, combined with kinetic energy, is found to be comparable to gravitational potential energy.

In the present work, we have estimated magnetic field strengths and derived various parameters for the CC and NES regions. However, we want to stress that these results must be taken with caution as the measurements are uncertain within a factor two due to the inherent large uncertainty in the mass and density of the studied regions. Moreover, the modified DCF methods can be uncertain up to a factor of 2 or more (Crutcher 2012), and also, the B-field strength in the studied regions can be biased due to the limited number of B-field segments traced by our observations. Due to all these uncertainties, the mass-to-flux ratio and virial status of the regions should be considered as qualitative indicators of the stability of the region. In the present case, we have used the generally accepted Q value of 0.5 (Crutcher 2012) in our estimations, however, if we use Q = 0.4, suggested for parsec scale clumps (Padoan et al. 2001), similar to our studied regions, the B-field strengths will reduce by a factor of ∼1.2. As a consequence, the CC and NES regions would become more magnetically supercritical. Future high-resolution and more sensitive observations would better constrain the magnetic field and turbulence properties of the hub. However, taking the measured mass-to-flux ratio and virial parameters at face value, it can be argued that, at present, gravity has overall an upper hand over magnetic and kinetic energies in CC and NES, which is consistent with the formation of a young cluster noticed by Rawat et al. (2023) in the hub.

4.3 The criticality of magnetic fields in hub-filamentary clumps

Like the CC of G148.24+00.41, the clumps located at the junction of the hub-filamentary systems are known to be potential sites of cluster formation (e.g. Kumar et al. 2020), because such clumps are attached to converging filamentary structures that fuel them with cold gaseous matter (e.g. Myers 2009; Li et al. 2014; Treviño-Morales et al. 2019; Kumar et al. 2022). Although physical processes that govern the star formation are scale-dependent, it is worthwhile to compare the global properties of the parsec or sub-parsec scale hub filamentary clumps studied with similar resolution. In the literature, a few such clumps have been studied with JCMT/POL2, these are: IC 5146 E-hub (Wang et al. 2019), G33.92 + 0.11 (Wang et al. 2020b), Mon R2 (Hwang et al. 2022), IC 5146 W-hub (Chung et al. 2022), and SDC13 (Wang et al. 2022). In the majority of these hub filamentary clumps, except Mon R2, the mass-to-flux ratio and/or virial analysis suggest the dominance or edge of gravitational energy over the magnetic and kinetic energies, similar to the CC region of G148.24+00.41. We found that all the aforementioned clumps are associated with either protostars or an embedded cluster (e.g. Harvey et al. 2008; Gutermuth et al. 2009; Peretto et al. 2014; Liu et al. 2015). Thus, it seems that, at least for those parsec scale hubs that are in the early stages of cluster formation, like the CC of G148.24+00.41, gravity has an upper hand on the energy budget of the system. A large sample of hub-filamentary clumps of various evolutionary stages (e.g. from pre-stellar clumps to clumps hosting emerging clusters of different ages) would be valuable to study the time evolution of various physical processes that govern star formation and its evolution.

5 SUMMARY

We have performed the dust polarization observation of the central part of the G148.24+00.41 cloud to investigate the B-field morphology and its strength relative to gravity and turbulence, using JCMT SCUBA-2/POL-2 at 850 |$\mu$|m. We specifically focused on the C1 clump, located at the hub of G148.24+00.41. The main results are summarised below:

  • The 850 |$\mu$|m Stokes I intensity map reveals the presence of a central clump (CC), northeastern elongated structure (NES), and western clump (WC) around the hub of G148.24+00.41. The B-field segments of CC and NES regions show mixed morphology, while the WC region shows converging B-field segments, mostly aligned along the south-east direction.

  • We found evidence of the depolarization effect, and from the Bayesian analysis over the non-debiased polarization data, found a power-law index, α = 0.6. Although this shows a decreasing level of dust grain alignment, but they can still be aligned with the magnetic field in the central high-density region of the cloud.

  • We compared the relative orientations of B-fields (θB), intensity gradients (θIG), and local gravity (θLG) over the full map. In the CC and WC regions, the three factors are mostly correlated, while the difference in orientations is higher in the NES region. This suggests that gravity is dragging the intensity gradients and aligning them with the B-fields in the CC and WC regions, while the effect of gravity in NES is comparatively less significant.

  • We constructed the ΣB map to see the localized B-field strength in comparison to local gravity, and found that for most of the parts, ΣB < 1, i.e. gravitational force is dominant over the magnetic field force.

  • We used the structure function analysis to determine the B-field strength, and as a result, found Bpos for CC and NES to be around ∼24.0 ± 6.0 and 20.0 ± 5.0 |$\mu$|G, respectively. In both regions, the turbulent component relative to the organized magnetic field structure was determined to be approximately 0.4–0.5.

  • The mass-to-flux ratio and Alfv|$\acute{\text{e}}$|n Mach number calculation for CC and NES shows that both regions are magnetically transcritical/supercritical and trans-Alfv|$\acute{\text{e}}$|nic. The turbulent pressure was found to be higher than the magnetic pressure in CC, while they are similar in NES.

  • The virial analysis shows that for CC, the |EG| > EK > EB (i.e. 10.4:4.7:3.3), while for NES, |EG| > EKEB (i.e. 14:6.3:7.0). The magnetic field and turbulence individually are not strong enough to provide stability to the regions against gravity. Both regions were found to be bound by gravity.

    Overall, we find that currently, gravitational energy has an edge over the other energy terms of the hub region of G148.24+00.41, thereby continue to facilitate the growth of the young cluster in the hub, although we acknowledge that given the large uncertainties associated with our estimates, a conclusive answer would require further precise measurements of magnetic field and cloud properties.

ACKNOWLEDGEMENTS

We thank the anonymous referee for the comments and suggestions that helped to improve the paper. The research work at the Physical Research Laboratory is funded by the Department of Space, Government of India. CE acknowledges the financial support from grant RJF/2020/000071 as a part of the Ramanujan Fellowship awarded by the Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the National Astronomical Research Institute of Thailand; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations in the United Kingdom, Canada, and Ireland. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. The data taken in this paper were observed under the project code M22BP055. DKO acknowledges the support of the Department of Atomic Energy, Government of India, under project identification no. RTI 4002. We thank Junhao Liu for the discussion on the DCF methods. This research made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a northern galactic plane CO survey with the PMO-13.7m telescope. We are grateful to all the members of the MWISP working group, particularly the staff members at PMO-13.7m telescope, for their long-term support. MWISP was sponsored by the National Key R&D Program of China with grant 2017YFA0402701 and CAS Key Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047.

Facilities: JCMT, PMO

DATA AVAILABILITY

The JCMT data used in this work can be shared on reasonable request. We also used the CO molecular data from PMO, which can be shared by the PMO database on reasonable request to the project team.

Footnotes

References

Andersson
B. G.
,
Lazarian
A.
,
Vaillancourt
J. E.
,
2015
,
ARA&A
,
53
,
501

André
P.
,
Könyves
V.
,
Arzoumanian
D.
,
Palmeirim
P.
,
Peretto
N.
,
2013
, in
Kawabe
R.
,
Kuno
N.
,
Yamamoto
S.
eds,
ASP Conf. Ser., Vol. 476, New Trends in Radio Astronomy in the ALMA Era: The 30th Anniversary of Nobeyama Radio Observatory
.
Astron. Soc. Pac
,
San Francisco
, p.
95

Ballesteros-Paredes
J.
,
Klessen
R. S.
,
Mac Low
M. M.
,
Vazquez-Semadeni
E.
,
2007
, in
Reipurth
B.
,
Jewitt
D.
,
Keil
K.
, eds,
Protostars and Planets V
. p.
63
,
preprint
()

Battersby
C.
et al. ,
2011
,
A&A
,
535
,
A128

Bertoldi
F.
,
McKee
C. F.
,
1992
,
ApJ
,
395
,
140

Beuther
H.
et al. ,
2018
,
A&A
,
614
,
A64

Beuther
H.
et al. ,
2020
,
ApJ
,
904
,
168

Beuther
H.
,
Vlemmings
W. H. T.
,
Rao
R.
,
van der Tak
F. F. S.
,
2010
,
ApJ
,
724
,
L113

Busquet
G.
,
2020
,
Nat. Astron.
,
4
,
1126

Chandrasekhar
S.
,
Fermi
E.
,
1953
,
ApJ
,
118
,
113

Chapin
E. L.
,
Berry
D. S.
,
Gibb
A. G.
,
Jenness
T.
,
Scott
D.
,
Tilanus
R. P. J.
,
Economou
F.
,
Holland
W. S.
,
2013
,
MNRAS
,
430
,
2545

Chapman
N. L.
,
Goldsmith
P. F.
,
Pineda
J. L.
,
Clemens
D. P.
,
Li
D.
,
Krčo
M.
,
2011
,
ApJ
,
741
,
21

Chung
E. J.
,
Lee
C. W.
,
Kwon
W.
,
Tafalla
M.
,
Kim
S.
,
Soam
A.
,
Cho
J.
,
2023
,
ApJ
,
951
,
68

Chung
E. J.
,
Lee
C. W.
,
Kwon
W.
,
Yoo
H.
,
Soam
A.
,
Cho
J.
,
2022
,
AJ
,
164
,
175

Cox
N. L. J.
et al. ,
2016
,
A&A
,
590
,
A110

Crutcher
R. M.
,
2012
,
ARA&A
,
50
,
29

Crutcher
R. M.
,
Nutter
D. J.
,
Ward-Thompson
D.
,
Kirk
J. M.
,
2004
,
ApJ
,
600
,
279

Crutcher
R. M.
,
Wandelt
B.
,
Heiles
C.
,
Falgarone
E.
,
Troland
T. H.
,
2010
,
ApJ
,
725
,
466

Currie
M. J.
,
Berry
D. S.
,
Jenness
T.
,
Gibb
A. G.
,
Bell
G. S.
,
Draper
P. W.
,
2014
, in
Manset
N.
,
Forshay
P.
, eds,
ASP Conf. Ser. Vol. 485, Astronomical Data Analysis Software and Systems XXIII
.
Astron. Soc. Pac
,
San Francisco
, p.
391

Davis
L.
,
1951
,
Phys. Rev.
,
81
,
890

Deharveng
L.
et al. ,
2012
,
A&A
,
546
,
A74

Doi
Y.
et al. ,
2020
,
ApJ
,
899
,
28

Eswaraiah
C.
et al. ,
2020
,
ApJ
,
897
,
90

Eswaraiah
C.
et al. ,
2021
,
ApJL
,
912
,
L27

Federrath
C.
,
2015
,
MNRAS
,
450
,
4035

Federrath
C.
,
Klessen
R. S.
,
2012
,
ApJ
,
761
,
156

Fiege
J. D.
,
Pudritz
R. E.
,
2000
,
MNRAS
,
311
,
85

Friberg
P.
,
Bastien
P.
,
Berry
D.
,
Savini
G.
,
Graves
S. F.
,
Pattle
K.
,
2016
, in
Holland
W. S.
,
Zmuidzinas
J.
, eds,
Proc. SPIE Conf. Ser. Vol. 9914, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VIII
.
SPIE
,
Bellingham
, p.
991403

Girart
J. M.
,
Beltrán
M. T.
,
Zhang
Q.
,
Rao
R.
,
Estalella
R.
,
2009
,
Science
,
324
,
1408

Girart
J. M.
,
Frau
P.
,
Zhang
Q.
,
Koch
P. M.
,
Qiu
K.
,
Tang
Y.-W.
,
Lai
S.-P.
,
Ho
P. T. P.
,
2013
,
ApJ
,
772
,
69

Girart
J. M.
,
Rao
R.
,
Marrone
D. P.
,
2006
,
Science
,
313
,
812

Gómez
G. C.
,
Vázquez-Semadeni
E.
,
2014
,
ApJ
,
791
,
124

Gómez
G. C.
,
Vázquez-Semadeni
E.
,
Palau
A.
,
2021
,
MNRAS
,
502
,
4963

Gómez
G. C.
,
Vázquez-Semadeni
E.
,
Zamora-Avilés
M.
,
2018
,
MNRAS
,
480
,
2939

Gutermuth
R. A.
,
Megeath
S. T.
,
Myers
P. C.
,
Allen
L. E.
,
Pipher
J. L.
,
Fazio
G. G.
,
2009
,
ApJS
,
184
,
18

Harvey
P. M.
et al. ,
2008
,
ApJ
,
680
,
495

Heitsch
F.
,
Zweibel
E. G.
,
Mac Low
M.-M.
,
Li
P.
,
Norman
M. L.
,
2001
,
ApJ
,
561
,
800

Hennebelle
P.
,
2018
,
A&A
,
611
,
A24

Hennebelle
P.
,
Falgarone
E.
,
2012
,
A&A Rev.
,
20
,
55

Hennebelle
P.
,
Inutsuka
S.-i.
,
2019
,
Frontiers Astron. Space Sci.
,
6
,
5

Hildebrand
R. H.
,
Kirby
L.
,
Dotson
J. L.
,
Houde
M.
,
Vaillancourt
J. E.
,
2009
,
ApJ
,
696
,
567

Hoang
T.
,
Lazarian
A.
,
2008
,
MNRAS
,
388
,
117

Hoang
T.
,
Lazarian
A.
,
2014
,
MNRAS
,
438
,
680

Holland
W. S.
et al. ,
2013
,
MNRAS
,
430
,
2513

Houde
M.
,
Vaillancourt
J. E.
,
Hildebrand
R. H.
,
Chitsazzadeh
S.
,
Kirby
L.
,
2009
,
ApJ
,
706
,
1504

Hull
C. L. H.
et al. ,
2017
,
ApJ
,
847
,
92

Hwang
J.
et al. ,
2022
,
ApJ
,
941
,
51

Kauffmann
J.
,
Bertoldi
F.
,
Bourke
T. L.
,
Evans
N. J. I.
,
Lee
C. W.
,
2008
,
A&A
,
487
,
993

Kauffmann
J.
,
Pillai
T.
,
Goldsmith
P. F.
,
2013
,
ApJ
,
779
,
185

Klessen
R. S.
,
Glover
S. C. O.
,
2016
, in
Revaz
Y.
,
Jablonka
P.
,
Teyssier
R.
,
Mayer
L.
, eds,
Saas-Fee Advanced Course Vol. 43, Saas-Fee Advanced Course
.
Springer-Verlag
,
Heidelberg
, p.
85

Klessen
R. S.
,
Heitsch
F.
,
Low
M.-M. M.
,
2000
,
ApJ
,
535
,
887

Koch
P. M.
et al. ,
2014
,
ApJ
,
797
,
99

Koch
P. M.
,
Tang
Y.-W.
,
Ho
P. T. P.
,
2012a
,
ApJ
,
747
,
79

Koch
P. M.
,
Tang
Y.-W.
,
Ho
P. T. P.
,
2012b
,
ApJ
,
747
,
80

Koch
P. M.
,
Tang
Y.-W.
,
Ho
P. T. P.
,
2013
,
ApJ
,
775
,
77

Kumar
M. S. N.
,
Arzoumanian
D.
,
Men’shchikov
A.
,
Palmeirim
P.
,
Matsumura
M.
,
Inutsuka
S.
,
2022
,
A&A
,
658
,
A114

Kumar
M. S. N.
,
Palmeirim
P.
,
Arzoumanian
D.
,
Inutsuka
S. I.
,
2020
,
A&A
,
642
,
A87

Kwon
J.
et al. ,
2018
,
ApJ
,
859
,
4

Lazarian
A.
,
2007
,
J. Quant. Spec. Radiat. Transf.
,
106
,
225

Li
H. B.
,
Goodman
A.
,
Sridharan
T. K.
,
Houde
M.
,
Li
Z. Y.
,
Novak
G.
,
Tang
K. S.
,
2014
, in
Beuther
H.
,
Klessen
R. S.
,
Dullemond
C. P.
,
Henning
T.
, eds,
Protostars and Planets VI
.
University of Arizona
,
Tucson
, p.
101
,

Li
H.-b.
,
Houde
M.
,
2008
,
ApJ
,
677
,
1151

Lin
S.-J.
et al. ,
2023
,
preprint
()

Liu
H. B.
,
Galván-Madrid
R.
,
Jiménez-Serra
I.
,
Román-Zúñiga
C.
,
Zhang
Q.
,
Li
Z.
,
Chen
H.-R.
,
2015
,
ApJ
,
804
,
37

Liu
J.
et al. ,
2019
,
ApJ
,
877
,
43

Liu
J.
,
Zhang
Q.
,
Qiu
K.
,
2022
,
Frontiers Astron. Space Sci.
,
9
,
943556

Liu
J.
,
Zhang
Q.
,
Qiu
K.
,
Liu
H. B.
,
Pillai
T.
,
Girart
J. M.
,
Li
Z.-Y.
,
Wang
K.
,
2020
,
ApJ
,
895
,
142

Mac Low
M.-M.
,
Klessen
R. S.
,
2004
,
Rev. Mod. Phys.
,
76
,
125

Mairs
S.
et al. ,
2021
,
AJ
,
162
,
191

Mao
S. A.
,
Ostriker
E. C.
,
Kim
C.-G.
,
2020
,
ApJ
,
898
,
52

Molinari
S.
et al. ,
2010
,
A&A
,
518
,
L100

Mouschovias
T. C.
,
Tassis
K.
,
Kunz
M. W.
,
2006
,
ApJ
,
646
,
1043

Myers
P. C.
,
2009
,
ApJ
,
700
,
1609

Nakamura
F.
,
Li
Z.-Y.
,
2008
,
ApJ
,
687
,
354

Naranjo-Romero
R.
,
Zapata
L. A.
,
Vázquez-Semadeni
E.
,
Takahashi
S.
,
Palau
A.
,
Schilke
P.
,
2012
,
ApJ
,
757
,
58

Ostriker
E. C.
,
Stone
J. M.
,
Gammie
C. F.
,
2001
,
ApJ
,
546
,
980

Padoan
P.
,
Goodman
A.
,
Draine
B. T.
,
Juvela
M.
,
Nordlund
Å.
,
Rögnvaldsson
Ö. E.
,
2001
,
ApJ
,
559
,
1005

Padoan
P.
,
Nordlund
Å.
,
2002
,
ApJ
,
576
,
870

Pattle
K.
et al. ,
2017
,
ApJ
,
846
,
122

Pattle
K.
et al. ,
2018
,
ApJL
,
860
,
L6

Pattle
K.
et al. ,
2019
,
ApJ
,
880
,
27

Pattle
K.
,
Fissel
L.
,
Tahani
M.
,
Liu
T.
,
Ntormousi
E.
,
2022
,
Protostars and Planets VII
,
preprint
()

Peretto
N.
et al. ,
2014
,
A&A
,
561
,
A83

Pillai
T. G. S.
et al. ,
2020
,
Nat. Astron.
,
4
,
1195

Pillai
T.
,
Kauffmann
J.
,
Wyrowski
F.
,
Hatchell
J.
,
Gibb
A. G.
,
Thompson
M. A.
,
2011
,
A&A
,
530
,
A118

Planck Collaboration XII
,
2020
,
A&A
,
641
,
A12

Planck Collaboration XXXV
,
2016
,
A&A
,
586
,
A138

Rawat
V.
et al. ,
2023
,
MNRAS
,
521
,
2786

Rawat
V.
,
Samal
M. R.
,
Walker
D. L.
et al. ,
2024
,
MNRAS
,
0035

Sadavoy
S. I.
et al. ,
2018
,
ApJ
,
869
,
115

Salvatier
J.
,
Wiecki
T. V.
,
Fonnesbeck
C.
,
2016
,
PeerJ Comput. Sci.
,
2
,
e55

Sanhueza
P.
,
Jackson
J. M.
,
Zhang
Q.
,
Guzmán
A. E.
,
Lu
X.
,
Stephens
I. W.
,
Wang
K.
,
Tatematsu
K.
,
2017
,
ApJ
,
841
,
97

Schisano
E.
et al. ,
2020
,
MNRAS
,
492
,
5420

Seifried
D.
,
Walch
S.
,
2015
,
MNRAS
,
452
,
2410

Shan
W.
et al. ,
2012
,
IEEE Trans. Terahertz Sci. Technol.
,
2
,
593

Shimajiri
Y.
,
André
P.
,
Palmeirim
P.
,
Arzoumanian
D.
,
Bracco
A.
,
Könyves
V.
,
Ntormousi
E.
,
Ladjelate
B.
,
2019
,
A&A
,
623
,
A16

Soam
A.
et al. ,
2018
,
ApJ
,
861
,
65

Soam
A.
et al. ,
2019
,
ApJ
,
883
,
95

Soler
J. D.
et al. ,
2017
,
A&A
,
603
,
A64

Soler
J. D.
,
Hennebelle
P.
,
Martin
P. G.
,
Miville-Deschênes
M. A.
,
Netterfield
C. B.
,
Fissel
L. M.
,
2013
,
ApJ
,
774
,
128

Su
Y.
et al. ,
2019
,
ApJS
,
240
,
9

Tan
J. C.
,
Kong
S.
,
Butler
M. J.
,
Caselli
P.
,
Fontani
F.
,
2013
,
ApJ
,
779
,
96

Tang
Y.-W.
,
Ho
P. T. P.
,
Koch
P. M.
,
Guilloteau
S.
,
Dutrey
A.
,
2013
,
ApJ
,
763
,
135

Tang
Y.-W.
,
Koch
P. M.
,
Peretto
N.
,
Novak
G.
,
Duarte-Cabral
A.
,
Chapman
N. L.
,
Hsieh
P.-Y.
,
Yen
H.-W.
,
2019
,
ApJ
,
878
,
10

Treviño-Morales
S. P.
et al. ,
2019
,
A&A
,
629
,
A81

Vaillancourt
J. E.
,
2006
,
PASP
,
118
,
1340

Vázquez-Semadeni
E.
,
Palau
A.
,
Ballesteros-Paredes
J.
,
Gómez
G. C.
,
Zamora-Avilés
M.
,
2019
,
MNRAS
,
490
,
3061

Wang
J.-W.
et al. ,
2019
,
ApJ
,
876
,
42

Wang
J.-W.
et al. ,
2022
,
ApJ
,
931
,
115

Wang
J.-W.
,
Koch
P. M.
,
Galván-Madrid
R.
,
Lai
S.-P.
,
Liu
H. B.
,
Lin
S.-J.
,
Pattle
K.
,
2020b
,
ApJ
,
905
,
158

Wang
J.-W.
,
Lai
S.-P.
,
Clemens
D. P.
,
Koch
P. M.
,
Eswaraiah
C.
,
Chen
W.-P.
,
Pandey
A. K.
,
2020a
,
ApJ
,
888
,
13

Ward-Thompson
D.
et al. ,
2017
,
ApJ
,
842
,
66

Whittet
D. C. B.
,
Hough
J. H.
,
Lazarian
A.
,
Hoang
T.
,
2008
,
ApJ
,
674
,
304

Zamora-Avilés
M.
,
Ballesteros-Paredes
J.
,
Hartmann
L. W.
,
2017
,
MNRAS
,
472
,
647

Zavagno
A.
et al. ,
2023
,
A&A
,
669
,
A120

Zhang
Q.
et al. ,
2014
,
ApJ
,
792
,
116

APPENDIX A: INTENSITY GRADIENTS VERSUS LOCAL GRAVITATIONAL FIELD

Fig. A1(a) shows the relative orientations of intensity gradient and local gravity, and Fig. A1 b shows the distribution of the angular difference between their orientations, i.e. ΔθIG, LG. Similar to ΔθB, IG and ΔθB, LG, the correlation between the angles (θIG and θLG) is better in the CC and WC regions in comparison to NES region.

(a) The orientations of the intensity gradients (red segments) and local gravity (cyan vectors) are overlaid on the 850 $\mu$m Stokes I map. (b) The distribution of the offset between the position angles of the intensity gradients and local gravity, i.e. ΔθIG, LG = |(θIG − θLG)| over the contours of 850 $\mu$m Stokes I map. The contour levels are same as in Fig. 2.
Figure A1.

(a) The orientations of the intensity gradients (red segments) and local gravity (cyan vectors) are overlaid on the 850 |$\mu$|m Stokes I map. (b) The distribution of the offset between the position angles of the intensity gradients and local gravity, i.e. ΔθIG, LG = |(θIG − θLG)| over the contours of 850 |$\mu$|m Stokes I map. The contour levels are same as in Fig. 2.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.