# Magnetic Fields in Massive Star-forming Regions (MagMaR): Unveiling an Hourglass Magnetic Field in G333.46-0.16 Using ALMA

Piyali Saha, Patricio Sanhueza, Marco Padovani, Josep Girart, Paulo Cortés, Kaho Morii, Junhao Liu, Á. Sánchez-Monge, Daniele Galli, Shantanu Basu, et al.

► To cite this version:

Piyali Saha, Patricio Sanhueza, Marco Padovani, Josep Girart, Paulo Cortés, et al. Magnetic Fields in Massive Star-forming Regions (MagMaR): Unveiling an Hourglass Magnetic Field in G333.46-0.16 Using ALMA. The Astrophysical Journal Letters, 2024, 972 (1), pp.L6. 10.3847/2041-8213/ad660c hal-05346073

HAL Id: hal-05346073

https://cnrs.hal.science/hal-05346073v1

Submitted on 6 Nov 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Magnetic Fields in Massive Star-forming Regions (MagMaR): Unveiling an Hourglass Magnetic Field in G333.46–0.16 Using ALMA

Piyali Saha¹, Patricio Sanhueza¹,², Marco Padovani³, Josep M. Girart⁴,⁵, Paulo C. Cortés⁶,⁷, Kaho Morii¹,⁸, Junhao Liu¹, Á. Sánchez-Monge⁴,⁵, Daniele Galli³, Shantanu Basu⁹, Patrick M. Koch¹⁰, Maria T. Beltrán³, Shanghuo Li¹¹, Henrik Beuther¹¹, Ian W. Stephens¹², Fumitaka Nakamura¹,²,⁸, Qizhou Zhang¹³, Wenyu Jiao¹⁴,¹⁵, M. Fernández-López¹⁶, Jihye Hwang¹⁷, Eun Jung Chung¹⁷, Kate Pattle¹⁸, Luis A. Zapata¹⁹, Fengwei Xu¹⁴,¹⁵,²⁰, Fernando A. Olguin²¹, Ji-hyun Kang¹⁷, Janik Karoly¹⁸, Chi-Yan Law³, Jia-Wei Wang²¹, Timea Csengeri²², Xing Lu²³, Yu Cheng¹, Jongsoo Kim¹⁷, Spandan Choudhury¹⁷, Huei-Ru Vivien Chen²¹, and Charles L. H. Hull¹

¹National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan; s.piyali16@gmail.com, piyali.saha@nao.ac.jp

² Astronomical Science Program, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

³ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

⁴Institut de Ciències de l'Espai (ICE-CSIC), Campus UAB, Can Magrans S/N, E-08193 Cerdanyola del Vallés, Catalonia, Spain

⁵Institut d'Estudis Espacials de Catalunya (IEEC), Esteve Terradas 1, PMT-UPC, E-08860 Castelldefels, Catalonia, Spain

⁶Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile

⁷National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

⁸Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

⁹Department of Physics and Astronomy, University of Western Ontario, London, ON N6A 3K7, Canada

¹⁰Academia Sinica, Institute of Astronomy and Astrophysics, Taipei, Taiwan

¹¹Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

¹²Department of Earth, Environment, and Physics, Worcester State University, Worcester, MA 01602, USA

¹³Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

¹⁴Department of Astronomy, School of Physics, Peking University, Beijing, 100871, People's Republic of China

¹⁵Kavli Institute for Astronomy and Astrophysics, Peking University, Haidian District, Beijing 100871, People's Republic of China

¹⁶Instituto Argentino de Radioastronomía (CCT-La Plata, CONICET; UNLP; CICPBA), C.C. No. 5, 1894, Villa Elisa, Buenos Aires, Argentina

¹⁷Korea Astronomy and Space Science Institute, 776 Daedokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea

¹⁸Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

¹⁹Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, P.O. Box 3-72, 58090, Morelia, Michoacán, Mexico

²⁰I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany

²¹Institute of Astronomy and Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan

²²Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France

²³Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China

Received 2024 May 1; revised 2024 July 20; accepted 2024 July 22; published 2024 August 22

## Abstract

The contribution of the magnetic field to the formation of high-mass stars is poorly understood. We report the high angular resolution (~0′′3, 870 au) map of the magnetic field projected on the plane of the sky ($B_{POS}$) toward the high-mass star-forming region G333.46–0.16 (G333), obtained with the Atacama Large Millimeter/submillimeter Array at 1.2 mm as part of the Magnetic fields in Massive star-forming Regions survey. The $B_{POS}$ morphology found in this region is consistent with a canonical “hourglass” with an embedded flattened envelope in a perpendicular direction, which suggests a dynamically important field. This region is fragmented into two protostars that appear to be gravitationally bound in a stable binary system with a separation of ~1740 au. Interestingly, by analyzing $H^{13}CO^+$ ($J = 3–2$) line emission, we find no velocity gradient over the extent of the continuum, which is consistent with a strong field. We model the $B_{POS}$, obtaining a marginally supercritical mass-to-flux ratio of 1.43, suggesting an initially strongly magnetized environment. Based on the Davis–Chandrasekhar–Fermi method, the magnetic field strength toward G333 is estimated to be 5.7 mG. The absence of strong rotation and outflows toward the central region of G333 suggests strong magnetic braking, consistent with a highly magnetized environment. Our study shows that despite being a strong regulator, the magnetic energy fails to prevent the process of fragmentation, as revealed by the formation of the two protostars in the central region.

Unified Astronomy Thesaurus concepts: Dust continuum emission (412); Magnetic fields (994); Polarimetry (1278); Star formation (1569)

## 1. Introduction

Star formation is a complex process controlled by several factors, among which gravity, turbulence, and magnetic field play a key role. Magnetic fields are believed to oppose the

gravitational contraction and fragmentation of dense cores, thus delaying the formation of protostars (e.g., Mouschovias & Ciolek 1999; McKee & Ostriker 2007; Li et al. 2017; Palau et al. 2021; Hwang et al. 2022). However, magnetic fields can also help to channel cloud material toward overdense regions, thus acting as a catalyst in the formation of protostars (e.g., Soler & Hennebelle 2017; Hennebelle & Inutsuka 2019).

In a strongly magnetic environment, the morphology of the magnetic field remains preserved up to scales of dense molecular

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
cores of 0.01–0.1 pc (e.g., Qiu et al. 2014; Li et al. 2015; Cortés et al. 2021). Another possible scenario shows the dragging of the frozen-in magnetic field along with the gas material toward the dense core by the gravitational collapse. This phenomenon creates a pinching effect in the magnetic field lines, leading to an “hourglass”-like appearance (e.g., Girart et al. 2006, 2009; Rao et al. 2009; Tang et al. 2009; Stephens et al. 2013; Hull et al. 2014; Qiu et al. 2014; Koch et al. 2018; Maury et al. 2018; Beltrán et al. 2019; Kwon et al. 2019; Cortés et al. 2021; Huang et al. 2024). This specific structure may be partially due to projection and line-of-sight integration effects that lead to the observed dust polarization. It is of utmost importance to characterize hourglass patterns when observed, in order to better understand the initial conditions of star formation.

Despite being a significant regulator of star formation, the role of the magnetic field during the birth of massive stars is still poorly understood. To make progress, the survey Magnetic fields in Massive star-forming Regions (MagMaR) has been carried out. In MagMaR, a total of 30 high-mass star-forming regions were observed at 1.2 mm with the Atacama Large Millimeter/ submillimeter Array (ALMA). A few detailed characteristics of some targets, such as G5.89–0.39, IRAS 18089–1732, NGC 6334I(N), and NGC 6334I have been discussed by Fernández-López et al. (2021), Sanhueza et al. (2021), Cortés et al. (2021), and Cortés et al. (2024), respectively.

Out of the 30 targets, the magnetic field projected on the plane of the sky ($B_{POS}$) toward the high-mass star-forming region G333.46–0.16 (hereafter G333) shows the most extended hourglass morphology on a few thousand au scale (assuming a distance of 2.9 kpc; Lin et al. 2019). This target was studied as a part of the ATLASGAL survey of massive clumps by Csengeri et al. (2014) and Lin et al. (2019). The bolometric luminosity of this target is $4.4 \times 10^3 L_\odot$ (Lin et al. 2019). Based on its spectral energy distribution, Lin et al. (2019) estimated a dust temperature and mass of 25.2 K and 282 $M_\odot$, respectively, for this region. Here, we aim to assess the magnetic properties and investigate the importance of the magnetic field, gravity, and turbulence in the star-forming region G333.

## 2. Observations

ALMA polarimetric observations toward G333 were carried out on 2018 September 27 using ALMA Band 6 (1.2 mm) as a part of the MagMaR project (project ID: 2017.1.00101.S and 2018.1.00105.S; PI: P. Sanhueza). The maximum recoverable scale (MRS) was ~4".5. During the observing runs, J1650–5044 was used as phase calibrator, while J1427–4206 was used as a calibrator for bandpass, flux, and polarization. Details of the observational setup are discussed by Cortés et al. (2021), Fernández-López et al. (2021), and Sanhueza et al. (2021).

Linearly polarized dust continuum emission is identified within the the inner one-third area of the primary beam (24"). Bright emission lines have been eliminated from the Stokes *I* continuum following the method elaborated by Olguin et al. (2021). We used CASA 6.5.2 to perform self-calibration and imaging. Self-calibration of the Stokes *I* was performed with three iterations in phase with a final solution interval of 10 s. The self-calibration solutions were applied to the respective spectral cubes. The imaging of each Stokes parameter was performed separately employing the CASA task *tclean* using Briggs weighting with a robust parameter of 1. The resulting Stokes *I* image has an angular resolution of 0".318 × 0".290 with a position angle of 49°.4 (922 × 841 au). The sensitivities of the images are 160 µJy beam⁻¹

for the final Stokes *I* and 29 µJy beam⁻¹ for the final Stokes *Q* and *U*. The source is bright in continuum Stokes *I* emission, making the image dynamic range limited. However, the polarized emission (Stokes *Q* and *U*) is significantly weaker compared to the continuum Stokes *I*; therefore, it is unlikely to be dynamic-range-limited, making it possible to reach the thermal noise. The debiased linear polarized intensity, polarization fraction, and polarization angle images were constructed following Vaillancourt (2006).

The image of H¹³CO⁺ (*J* = 3–2) line emission, also included in the spectral setup, was obtained by the automatic masking method *yclean* (Contreras et al. 2018). The CASA task *tclean* was performed using Briggs weighting with a robust parameter of 1, leading to a noise level of 2.9 mJy beam⁻¹ (0.61 K) per 0.56 km s⁻¹ channel.

## 3. Results

The ALMA 1.2 mm continuum emission with ~900 au spatial resolution is shown in Figure 1(a) and reveals details of the internal structure of G333. The dust continuum emission shows a flattened structure that is elongated in the northwest-southeast direction. The center of this flattened structure is fragmented into two more condensations with a separation of ~1740 au along the major axis of the elongation. Both the condensations show the presence of hot molecular core emission lines (Taniguchi et al. 2023), suggesting that they are potential protostars. The peak flux at the position of the brightest component, located toward the southeast (MM1), is 53.5 mJy beam⁻¹, while the fainter one located toward the northwest (MM2) has a peak flux of 38.8 mJy beam⁻¹.

The direction of the $B_{POS}$ is inferred by assuming the dust grains are aligned with respect to the magnetic fields (i.e., rotating the polarization segments by 90°; Cudlip et al. 1982; Hildebrand et al. 1984; Hildebrand 1988; Lazarian 2000; Andersson et al. 2015). The polarized emission detected in G333 suggests an hourglass-like geometry of the magnetic field aligned with the symmetry axis of the hourglass almost parallel to the minor axis of the flattened envelope. For our analysis, we focus on the circular area as marked in Figure 1(a) that harbors the hourglass magnetic field. We chose this area to avoid the additional distortion of $B_{POS}$ produced by other surrounding cores. The circle is centered at $\alpha = 16:21:20.183$ and $\delta = -50:09:46.662$, with a radius of 1".5 (4350 au). As the diameter of the area of analysis (3") is less than the MRS, the extended emission should not be significantly affected by filtering.

Figure 2(a) presents the integrated intensity map of H¹³CO⁺ (*J* = 3–2) line, which is a cold dense gas tracer (upper energy level, $E_u$, of 25 K), toward G333. Figure 2(b) shows the H¹³CO⁺ distribution of the intensity-weighted velocity structure (moment 1 map) toward G333 with respect to the systemic velocity ($v_{lsr}$ ~ −43 km s⁻¹). No clear signature of a large-scale velocity gradient is detected with a spectral resolution of 0.56 km s⁻¹, suggesting a quiescent environment at ~10,000 au scale. However, the velocity field over MM1 and MM2 is relatively more blueshifted than their immediate vicinities, which could be a sign of infall, as suggested by Estalella et al. (2019) and Olguin et al. (2021).

## 4. Discussion

### 4.1. Dust Continuum Emission and Velocity Dispersion

The flux density enclosed in the circular region is 493 mJy. The temperature of the region marked with a circle as in
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) The $B_{POS}$ geometry (black line segments) toward G333, obtained after rotating the polarization segments by 90°, overplotted on ALMA 1.2 mm dust continuum emission. The $B_{POS}$ segments have arbitrary length and are plotted above the $3\sigma$ ($\sigma = 29~\mu\text{Jy beam}^{-1}$) for Stokes Q and U level. Contours correspond to dust continuum emission 5, 10, 50, 100, 190, 230, 300, and 340 times $\sigma$ (=160 $\mu\text{Jy beam}^{-1}$). The positions of MM1 and MM2 are marked as plus symbols. The circle in magenta represents the area of analysis. The scale bar and the beam size of 880 au are displayed on the bottom right and bottom left sides, respectively. The line segments are drawn following the Nyquist sampling (every 3 pixels each of size $0''05$). (b) Magnified view of the area of analysis. Symbols are the same as in panel (a).</figcaption></figure>

Figure 1 is determined using the CH$_3$CN ($J=14-13$) rotational transitions. We primarily detected the four brightest K components ($K=0, 1, 2, 3$) of the transition toward this area. We avoided the region very close to MM1 and MM2 to refrain from line contamination by radiation from MM1 and MM2. In Figure 3, we present the average brightness distribution of the CH$_3$CN ($J=14-13$) transitions in black. The blue dashed line represents the fitted spectrum that is obtained by fitting the observed spectrum with XCLASS (Möller et al. 2017) resulting in a temperature of 50 K. Using this temperature, and assuming optically thin dust emission and a spherical geometry, the total gas mass is estimated using

$$M = \Lambda \frac{F_{\nu} D^2}{\kappa_{\nu} B_{\nu}(T)}, \qquad (1)$$

where $\Lambda$ is the gas-to-dust mass ratio, $F_\nu$ is the flux density of the source, D is the distance to the target, $\kappa_\nu$ is the dust opacity per gram of dust, and $B_\nu$ is the Planck function with a dust temperature T. The gas-to-dust ratio is assumed to be 100:1. The dust opacity is assumed as $1.03 \text{ cm}^2 \text{ g}^{-1}$ (interpolated to 1.2 mm; Ossenkopf & Henning 1994). Using a dust temperature of 50 K, the total mass is estimated as $23 M_\odot$.

We estimate the mass density $\rho$ as $3.9 \times 10^{-17} \text{ gm cm}^{-3}$ using a volume ($V = \frac{4}{3}\pi R^3$, where $R = 1''5$). The number density $n(H_2)$ is $\rho/\mu_{H_2}m_H$, where $\mu_{H_2}$ (=2.86) is the mean molecular weight per hydrogen molecule (Kirk et al. 2013) and $m_H$ is the atomic mass of hydrogen. The $n(H_2)$ enclosed in the area of analysis is estimated as $8.4 \times 10^6 \text{ cm}^{-3}$.

Using the dendogram technique (Rosolowsky et al. 2008) adopted in the astrodendro Python package,<fnref n="24" /> the flux densities

of MM1 and MM2 are estimated as 59.5 mJy and 25.8 mJy, respectively. Assuming a similar temperature, MM1 is 2.3 times more massive than MM2. Table 1 lists the dendogram results for both protostars.

Considering a bolometric luminosity of $4.4 \times 10^3 L_\odot$ (Lin et al. 2019) and assuming that a single main-sequence star is responsible for the bolometric luminosity, the stellar mass would then be ~6 $M_\odot$ (Mottram et al. 2011). However, there are two protostars at the center. At larger scales, based on the ALMA 7 m array observations with an angular resolution of $3''7$ (Csengeri et al. 2017), the mass of the structure containing MM1 and MM2 is determined as $120 M_\odot$, suggesting that most of the clump mass ($282 M_\odot$) is concentrated in the central region of the clump. Thus, these protostars are potential candidates for becoming high-mass stars by accreting gas material from the available mass reservoir.

Using H$^{13}$CO$^+$ line emission, we derive the dispersion in the turbulent velocity along the line of sight, $\sigma v_{los}$ (= $\sqrt{\sigma_{obs}^2 - \sigma_{th}^2}$). The $\sigma_{obs}$ and $\sigma_{th}$ are the observed and thermal velocity dispersions of the H$^{13}$CO$^+$ line, respectively. The $\sigma_{th}$ can be expressed as $\sqrt{k_B T/\mu m_H}$, where $k_B$ is the Boltzmann constant and $\mu$ is the molecular weight of H$^{13}$CO$^+$ (30). Assuming a temperature of 50 K, $\sigma_{th}$ and $\sigma v_{los}$ are estimated as 0.12 and 1.21 km s$^{-1}$, respectively.

## 4.2. Modeling of the Hourglass-shaped Magnetic Field

In order to model the magnetic field pattern, we used the DustPol module included in the ARTIST package (Padovani et al. 2012), which is based on the Line Modeling Engine radiative transfer code (Brinch & Hogerheijde 2010). DustPol generates the synthetic Stokes I, Q, and U images in FITS format, which are used as an input for the simobserve and simanalyze tasks of CASA, considering

<footnote n="24">²⁴ http://www.dendrograms.org/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) The moment zero or integrated intensity map of H$^{13}$CO$^{+}$ line emission. (b) The moment 1 or intensity-weighted velocity map of the same. The white line segments show the B$_{POS}$ geometry, and the contours outline the dust continuum emission (same as Figure 1). The positions of MM1 and MM2 are marked as plus symbols. The scale bar and the spatial resolution are displayed on the bottom right and bottom left sides, respectively, in all panels.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The spectra of CH$_{3}$CN (J = 14–13) toward G333. The black line shows the observed spectrum, and the dashed blue line shows the fitted spectrum.</figcaption></figure>

the same antenna configuration of the observing runs in ALMA. In this way, polarization position angle maps are created and compared with the observations. The physical modeling of the 3D magnetic field is done by combining an axisymmetric singular toroid threaded by a poloidal field (Li & Shu 1996; Padovani & Galli 2011). Based on Padovani et al. (2013), we added a toroidal force-free component of the magnetic field to represent the effects of rotation.

<table><caption>**Table 1**

Results Obtained from Astrodendro<sup>a</sup></caption><thead><tr><th>R.A.<br>(ICRS)</th><th>Decl.<br>(ICRS)</th><th>FWHM<br>(arcsec × arcsec)</th><th>Peak<br>Intensity<br>(mJy<br>beam<sup>-1</sup>)</th><th>Flux<br>Density<br>(mJy)</th></tr></thead><tbody><tr><td rowspan="2">16:21:20.20</td><td rowspan="2">−50:09:46.91</td><td>MM1<br>0.31 × 0.21</td><td rowspan="2">53.5</td><td rowspan="2">59.5</td></tr><tr><td>MM2<br>0.23 × 0.14</td></tr><tr><td>16:21:20.17</td><td>−50:09:46.41</td><td></td><td>38.8</td><td>25.8</td></tr></tbody></table>

**Note.**

<sup>a</sup> F$_{min}$ = 5σ: the threshold above which the structures are defined. δ = 1σ: the minimum significance for separation of the structures. S$_{min}$ = half of the beam size: the minimum area to be allocated in the individual structure.

The model used in this work has four free parameters: (1) the mass-to-flux ratio normalized to its critical value, λ, defined as

$$ \lambda = \frac{(M/\Phi_B)}{(M/\Phi_B)_{cr}} = 2\pi\sqrt{G} \frac{M}{\Phi_B}, \quad (2) $$

where G is the gravitational constant, Φ$_{B}$ is the magnetic flux (=πR$^{2}$B), and M is the mass of the core; (2) the ratio of the strengths of the toroidal and the poloidal components of the magnetic field in the midplane of the source, b$_{0}$; (3) the orientation of the magnetic axis projected on the plane of the sky, ϕ, starting from north and increasing eastward; and (4) the inclination with respect to the plane of the sky, i, with an assumption of being positive/negative when the magnetic field in the northern part directs toward/away from us. A χ$^{2}$ test is performed on the polarization angle residuals (Δψ) of the polarization angle ψ, obtained from the difference between the observed (ψ$_{obs}$) and the modeled (ψ$_{mod}$) polarization angles (Δψ = ψ$_{obs}$ − ψ$_{mod}$), enclosed within the circular area for each combination of the four free parameters. The χ$^{2}$ is estimated accounting for regions lying above the 3σ level (σ = 29 µJy beam$^{-1}$) of the Stokes Q and U images. A constant temperature (50 K) has been adopted in the modeling, as no temperature map of this region is available currently. On the other hand, the region close to MM1 and MM2 is expected to be relatively hot and have more weight in estimating the Stokes parameters.

The best-fit model provides the minimum reduced χ$^{2}$ value of 9.39 with the combination of model parameters: λ = 1.63, b$_{0}$ = −0.1 ± 0.1, i = 45°+13°−6°, and ϕ = 40°+1°−4°. The errors are estimated following Lampton et al. (1976). A low value of b$_{0}$ in the best-fit model is consistent with no apparent rotation toward G333 as found in H$^{13}$CO$^{+}$. The value of λ = 1.63 obtained from the model (see Equation (2)) is relative to the mass enclosed in a flux tube, but observationally, the mass is derived considering a spherical system. Therefore, the estimated λ is transformed to an effective λ of 1.43 (Li & Shu 1996), which is marginally supercritical.

In Figure 4(a), we show a comparison between the observed and modeled magnetic field geometry, limited to the area within the circle shown in Figure 1(a). Although the hourglass model reproduces quite well the observed magnetic field geometry, deviations can be seen around MM1, mainly toward its southern and eastern regions. Since MM1 is relatively more massive than MM2, it is possible that the gravity of MM1 might have dragged and distorted the magnetic field lines more effectively than MM2. In Figure 4(b), we show the histogram of the Δψ. A Gaussian fit
<figure><img src="image_7.png" /><figcaption>Figure 4. (a) Comparison between the observed (green segments) and modeled (yellow segments) magnetic field orientations overplotted on the polarized intensity map (gray scale). The black contours are the same as Figure 1. (b) Histogram of the polarization angle residuals (Δψ) along with the Gaussian fitting.</figcaption></figure>

to the whole distribution of Δψ provides a mean value of $\langle \Delta\psi \rangle = -3.32$ and a standard deviation of $\sigma_{\psi} = 19.48$.

### 4.3. Magnetic Field Strength

We estimate $B_{POS}$ toward G333 using the well-known Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953),<fnref n="25" /> which states

$$B_{POS} = Q \sqrt{4\pi\rho} \frac{\sigma v_{los}}{\delta\psi_{int}}, \quad (3)$$

where $Q$ is a correction factor, 0.5,<fnref n="26" /> adopted from simulations in turbulent clouds (Ostriker et al. 2001); $\rho$ is the mass density of the cloud; $\sigma v_{los}$ is the dispersion in the turbulent velocity

along the line of sight; and $\delta\psi_{int}$ is the intrinsic dispersion in the polarization angle $\psi$. The $\delta\psi_{int}$ can be obtained from $\sqrt{\sigma_{\psi}^2 - \delta\psi_{obs}^2}$. The average uncertainty in the observed $\psi$ within the circle of analysis is $\delta\psi_{obs} = 5.4$. The $\delta\psi_{int}$ obtained is 18.7. Consequently, $B_{POS}$ is estimated as 4.0 mG. Adopting $i = 45^\circ$ from our model, we also estimate the total magnetic field strength $B = B_{POS}/\cos(i) = 5.7$ mG. With a density of $n(H_2) = 8.4 \times 10^6$ cm$^{-3}$, this strength is consistent as found in the compilation of B-field strength with density by Liu et al. (2022a) using the DCF method.

The DCF method assumes the deviation of the uniform magnetic field by turbulent gas motions without considering gravity, which is most likely responsible for the hourglass shape. In this work, gravity is taken into account while using the DustPol model, and by considering the Δψ, we eliminate the gravitational effect, and the $\delta\psi$ is supposed to be the interplay between the turbulence and the magnetic field.

While computing $\lambda$ using Equation (2), we assumed that all the mass is in the envelope, while the protostellar mass should also be added to the envelope mass. Therefore, the estimated $\lambda$ should be considered as a lower limit. As mentioned in Section 4.1, based on the luminosity, an additional stellar mass of ~6 $M_\odot$ could be added to the envelope mass, making the total mass 29 $M_\odot$, which leads to a mass-to-flux ratio of 1.22. The fact that the dust emission may be optically thick at the very central bright fragmented region leads to an underestimate of the mass of the circular region that cannot be ignored. However, the estimated value of $\lambda$ is similar to the one obtained from the modeling (1.43).

We compute the Alfvén speed, $v_A = B/\sqrt{4\pi\rho} = 2.58$ km s$^{-1}$, through which we estimate the turbulent to magnetic energy, $\beta_t = 3(\sigma v_{los}/v_A)^2 = 0.66$. This indicates that the magnetic energy dominates moderately over the turbulent energy in this region. Thus, the magnetic field toward G333, in spite of being quite strong, could not prevent the fragmentation and eventually star formation.

Based on nonideal magnetohydrodynamic simulations of an initially subcritical core by Machida & Basu (2020), the collapse happens after the core reaches a marginally supercritical mass-to-flux ratio, leaving behind a considerable amount of magnetic flux, which helps to transport out the angular momentum through strong magnetic braking. As a result, the protostar may harbor no disk or a very small one, with a very weak outflow signature. In G333, we do not find any evidence of rotation in the envelope region (Figure 2(b)) or in any of the protostars at the current angular resolution. Even in higher angular resolution (~0.05°) observations of the same region obtained as part of the survey Digging into the Interior of Hot Cores with ALMA (DIHCA; Olguin et al. 2022, 2023) survey (P. Sanhueza, private communication), no clear signature of rotation was found in any of the protostars, making the argument of strong magnetic braking more evident. Additionally, we found no significant outflow signature coming from the two protostars in G333 (see Appendix A). This is consistent with the overall physical scenario of very weak outflow and no disk when starting with strongly magnetized initial conditions (Machida & Basu 2020). This leads to the possibility that star formation in G333 started in an initially subcritical or transcritical cloud.

<footnote n="25">²⁵ Also see Liu et al. (2022a, 2022b).</footnote>

<footnote n="26">²⁶ In Ostriker et al. (2001), the size scale adopted is ~8 pc, whereas Liu et al. (2021) estimated $Q$ to be ~0.28 at a size scale of 1–0.2 pc, which is relatively closer to the size scale of the area of analysis. Adopting the latter, the $B_{POS}$ is estimated as ~2.3 mG.</footnote>
## 4.4. Analysis of the Energy Balance

The current dynamical state toward G333 can be inferred from the virial parameter ($\alpha_{vir}$), the ratio of the virial mass ($M_{vir}$) to the total mass ($M$) of the system with considerations of both the turbulent and magnetic supports. While $\alpha_{vir} \sim 1$ indicates the equilibrium state, $\alpha_{vir} > 1$ implies expansion, and, conversely, $\alpha_{vir} < 1$ signifies contraction under gravity. We estimate $\alpha_{vir}$ as 1.3 for a centrally peaked density profile (see Appendix B). This suggests that G333 is currently close to a quasi-equilibrium state or undergoing quasi-static evolution, which is further supported by the quiescent environment of G333 as implied by the velocity field traced by the H$^{13}$CO$^{+}$ line emission (Figure 2(b)). A strongly magnetized environment exerts adequate support against rapid gravitational collapse, although the magnetic field could not hinder the fragmentation. This is consistent with a second stage of fragmentation within an initial massive core that forms out of a transcritical cloud (Bailey & Basu 2014).

## 4.5. Flattened Envelope Harboring the Hourglass Magnetic Field

Theoretical studies suggest that for a relatively strong magnetic field environment ($\lambda \le 3$), the core material is channeled primarily along the magnetic field lines by the gravitational pull, building a flattened pancake-like structure perpendicular to the field lines and developing a small toroidal component of the magnetic field (Allen et al. 2003; Price & Bate 2007). As shown in Figure 1, the flattened structure where the fragmentation took place agrees well with the strong magnetic field case, which is further supported by a low $\lambda$ obtained from our best-fit model as well as the observations. Additionally, a small contribution of the toroidal component of the field ($\sim 10\%$ of the poloidal component) obtained from the best-fit model supports this phenomenon. This is consistent with some previous studies (e.g., Qiu et al. 2014; Beltrán et al. 2019; Kwon et al. 2019), which suggest a major contribution of the magnetic field in the formation of high-mass protostars.

## 4.6. Stability Analysis of the Central Protostars

To investigate whether the central protostars MM1 and MM2 are gravitationally bound, we followed the methodology adopted in Pineda et al. (2015) and Li et al. (2024). We estimate the gravitational potential energy ($V_i$) by

$$V_i = -\sum_{i \neq j} \frac{Gm_i m_j}{r_{ij}}, \qquad (4)$$

where $m_i$ and $m_j$ are the masses of objects i and j, respectively, and $r_{ij}$ is the separation between them. The kinetic energy ($E_i$) is given by

$$E_i = \frac{1}{2} m_i (v_i - v_{\text{com}})^2, \qquad (5)$$

where $v_i$ is the line-of-sight velocity of the object i and $v_{\text{com}}$ is the velocity of the center of mass of the system. We estimate $v_{\text{com}}$ by

$$v_{\text{com}} = \frac{\sum_k m_k v_k}{\sum_k m_k}. \qquad (6)$$

The velocities of MM1 and MM2 are obtained from Taniguchi et al. (2023) as -44.06 and -43.55 km s$^{-1}$,

respectively. Assuming that the measured velocity difference between MM1 and MM2 is in 1D, the full velocity difference in 3D is $\Delta v_{3D} = \sqrt{3}(v_i - v_{\text{com}})$ along the line of sight. Also, the total separation between MM1 and MM2 is estimated by multiplying the measured projected separation (1740 au) with $4/\pi$, i.e., 2215 au, assuming a random orientation between them.

We measured the masses of MM1 and MM2 from the flux density obtained from the astrodendro task to estimate the $E_i/V_i$ for each of them. The $E_i/V_i$ for MM1 and MM2 are estimated as 0.08 and 0.17, respectively. As both MM1 and MM2 show $E_i/V_i < 1$, they could be gravitationally bound in a binary system at the present stage.

## 5. Summary and Conclusions

The ALMA high angular resolution ($\sim 0''3$) observations of linearly polarized 1.2 mm dust emission toward the high-mass star-forming region G333 revealed an hourglass-shaped pattern of $B_{POS}$ at a scale of $\sim 6000$ au. The hourglass shape is found to be more pinched toward MM1 than MM2, which might be because of the larger mass of MM1 compared to MM2. The protostars MM1 and MM2 are found to be gravitationally bound in a binary system at present with a separation of 1740 au.

The H$^{13}$CO$^{+}$ line emission shows no strong velocity gradient that could hint at any rotation toward G333. Also, none of the protostars show any clear signature of outflows, suggesting a strong magnetic braking.

The hourglass-shaped $B_{POS}$ is well fitted with an axisymmetric semianalytical magnetostatic model. The best-fit model is primarily governed by a poloidal component tangled with a small toroidal component ($\sim 10\%$ of the poloidal one). The best-fit results also provide $\lambda = 1.43$, $i = 45^{\circ} \pm 13^{\circ}$, and $\phi = 40^{\circ} \pm 4^{\circ}$.

Based on the DCF relation, the total magnetic field strength is estimated as 5.7 mG. We also estimate the mass-to-flux ratio using the dispersion of polarization angle obtained from the model fitting, which results in 1.22, consistent with the same obtained from the best-fit model ($\lambda = 1.43$).

Our analysis of energy balance in this area suggests a quasi-equilibrium state, which is consistent with the current dynamic environment inferred by H$^{13}$CO$^{+}$ line emission. Our speculation based on this work is that the magnetic field toward G333 is strong enough to stabilize the environment toward the center. A low magnitude of turbulent-to-magnetic energy ratio indicates a suppression of the magnetic field over the turbulence in this region. However, both are ultimately overwhelmed by the gravity, responsible for the fragmentation and eventually star formation. This leads to a possible scenario that the star formation in G333 might have started in an initially subcritical core, when the central peak density area reached the marginally supercritical mass-to-flux ratio.

## Acknowledgments

We thank the anonymous referee for providing constructive suggestions that enriched the draft significantly. P.S. was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Nos. JP22H01271 and JP24K17100) of the Japan Society for the Promotion of Science (JSPS). P.S. was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Nos. JP22H01271 and JP23H01221) of the JSPS. J.M.G. acknowledges support by grant PID2020-117710GB-I00 (MCI-AEI-FEDER, UE). This work is also partially supported by the program Unidad
de Excelencia María de Maeztu CEX2020-001058-M. P.C.C. was supported by the NAOJ Research Coordination Committee, NINS (NAOJ-RCC-2202-0401). J.L. is partially supported by a Grant-in-Aid for Scientific Research (KAKENHI No. JP23H01221) of JSPS. A.S.-M. acknowledges support from the RyC2021-032892-I grant funded by MCIN/AEI/10.13039/501100011033 and by the European Union “Next GenerationEU”/PRTR, as well as the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and support from the PID2020-117710GB-I00 (MCI-AEI-FEDER, UE). M.T.B. acknowledges financial support through the INAF Large Grant *The role of MAGnetic fields in MAssive star formation* (MAGMA). M.F.L. acknowledges support from the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement No. 101086388 (LACEGAL). M.F.L. also acknowledges the warmth and hospitality of the ICE-UB group of star formation. K.P. is a Royal Society University Research Fellow, supported by grant No. URF\R1\211322. X.L. acknowledges support from the National Key R&D Program of China (No. 2022YFA1603101), the Natural Science Foundation of Shanghai (No. 23ZR1482100), the National Natural Science Foundation of China (NSFC) through grant Nos. 12273090 and 12322305, and the Chinese Academy of Sciences (CAS) “Light of West China” Program (No. zbzgzdsys-202212). T.Cs. has received financial support from the French State in the framework of the IdEx Université de Bordeaux Investments for the Future Program. Y.C. was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI No. JP24K17103) of the JSPS. This Letter makes use of the following ALMA data: ADS/JAO. ALMA#2017.1.00101.S. Data analysis was (in part) carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research made use of astrodendro, a Python package to compute dendrograms of astronomical data (http://www.dendrograms.org/).

Facility: ALMA

Software: CASA (McMullin et al. 2007; CASA Team et al. 2022).

# Appendix A

## Molecular Outflows Identified around G333

Several molecular outflows are identified around G333 using the ¹²CO (J = 3–2) line emission observed by ALMA (project ID: 2013.1.00960.S; PI: T. Csengeri). Using the CASA task *imregrid*, we match the phase center of the ¹²CO line image and our 1.2 mm dust continuum image.

Figure 5 shows the ¹²CO molecular outflows overplotted on the dust continuum image. The ¹²CO line emission is integrated from −106.1 to −65.3 km s⁻¹ for the blueshifted lobe and −22.1–0.7 km s⁻¹ for the redshifted lobe. The sensitivities of the integrated images of the blueshifted and redshifted ¹²CO line emission are 0.14 and 0.10 Jy beam⁻¹, respectively. The outflow lobes in the north and south of the region enclosed by the hourglass-shaped magnetic field (shown as a cyan colored circle in Figure 5) are highly collimated. Interestingly, it is noticeable that although there is a significant amount of redshifted ¹²CO line emission within the circle (which seems to

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Molecular outflows identified in ¹²CO (J = 3–2) overplotted on our observed 1.2 mm dust continuum emission. The plus symbols mark the positions of MM1 and MM2. The circle in cyan shows the region of analysis. The blueshifted and redshifted outflows are shown in blue and red contours, respectively. The blue contours are drawn with levels 3, 5, 7, 9, 12, and 18 times 0.13 Jy beam⁻¹. The red contours are drawn with levels 3, 5, 7, 9, 11, 15, 20, and 25 times 0.10 Jy beam⁻¹. The possible sources driving the outflows are marked using green crosses. The red and blue arrows show the directions of the redshifted and blueshifted outflows, respectively.</figcaption></figure>

be the redshifted lobe of the nearby protostar located south of the circular area of analysis), none of the outflow lobes are driven by any of the protostars but likely the low-mass protostars residing in the outskirts of the dense gas. There are some dust condensations that likely could host these protostars, marked with crosses in Figure 5.

The absence of outflows in any of the protostars makes it difficult to predict the angular momentum direction. In a classical scenario of magnetized star formation, the angular momentum is efficiently lost by magnetic braking (Mouschovias & Paleologou 1979). Our speculation is that magnetic braking is relatively strong in the case of G333, because of which we do not find any observational signature of strong angular momentum.

# Appendix B
## Virial Parameter

The virial theorem can be expressed as follows:

$$ \frac{1}{2} \frac{d^2 I}{dt^2} = 2E_k + E_G + E_B, \quad (B1) $$

where *I* is the moment of inertia and $E_k$, $E_G$, and $E_B$ are the kinetic, gravitational, and magnetic energies, respectively. We consider a spherical cloud having a radial density profile $\rho \propto r^{-\alpha}$, where $\alpha = 2$ implies that the density profile is centrally peaked (Belloche 2013). The gravitational energy ($E_G$) can be expressed by

$$ E_G = - \frac{(3 - \alpha) GM^2}{(5 - 2\alpha) R}. \quad (B2) $$

The magnetic energy ($E_B$) is given by

$$ E_B = \frac{B^2 V}{8\pi}. \quad (B3) $$
The kinetic energy ($E_k$) is derived by

$$E_k = \frac{3}{2} M \sigma_{\text{tot}}^2. \qquad (B4)$$

where $\sigma_{\text{tot}}$ is the total gas velocity dispersion defined by $\sigma_{\text{tot}} = \sqrt{\sigma_{\text{th}}^2 + \sigma v_{\text{los}}^2}$. The $\sigma_{\text{th}}$ is the thermal velocity dispersion ($\sigma_{\text{th}} = \sqrt{k_B T / \mu_p m_H}$), where $\mu_p = 2.37$ is the mean molecular weight per free particle. The rotational energy ($E_r$) is expressed as

$$E_r = \frac{1}{3} M v_{\text{rot}}^2 \frac{(3 - \alpha)}{(5 - \alpha)}. \qquad (B5)$$

We do not find any evidence of rotation with a spectral resolution of 0.56 km s⁻¹ in the area of analysis. If we assume this resolution as the upper limit of the velocity due to rotation ($v_{\text{rot}}$), we estimate the ratio $E_r/E_G$ to be 0.006 for a centrally peaked density profile. With the same assumption of rotation, a ratio of $E_r/E_B$ is estimated to be 0.013 for the same profile. Similarly, we estimate the ratio of $E_r/E_k$ as 0.014. As the contribution of the rotational energy is negligible compared to the other energies, we do not account for the upper limit of the rotational term in the estimation of the virial parameter.

Including both the kinetic and magnetic energies, a system can be in stable phase, when $2E_k + E_B + E_G < 0$. The ratio between the virial mass $M_{\text{vir}}(=M_{k+B})$ and the total mass $M$ can be defined as the virial parameter $\alpha_{\text{vir}}$. Following Liu et al. (2020), $\alpha_{\text{vir}}$ is given as

$$\alpha_{\text{vir}} = \frac{M_{k+B}}{M} = \frac{1}{M} \left[ \sqrt{M_B^2 + \left(\frac{M_k^2}{2}\right) + \frac{M_k}{2}} \right], \qquad (B6)$$

where $M_{k+B}$ is the total virial mass considering the kinetic motion and magnetic field. The kinetic virial mass $M_k$ is defined as

$$M_k = 3 \left( \frac{5 - 2\alpha}{3 - \alpha} \right) \frac{\sigma_{\text{tot}}^2 R}{G}. \qquad (B7)$$

The virial mass accounting for the ordered magnetic field $M_B$ is expressed as

$$M_B = \frac{\pi R^2 B}{\sqrt{\frac{6(3 - \alpha)}{(5 - 2\alpha)} \pi^2 G}}. \qquad (B8)$$

We estimate the $\alpha_{\text{vir}}$ to be 1.3, which suggests an equilibrium state of G333, implying that the magnetic field provides significant support against gravity.

## ORCID iDs

Piyali Saha https://orcid.org/0000-0002-0028-1354

Patricio Sanhueza https://orcid.org/0000-0002-7125-7685

Marco Padovani https://orcid.org/0000-0003-2303-0096

Josep M. Girart https://orcid.org/0000-0002-3829-5591

Paulo C. Cortés https://orcid.org/0000-0002-3583-780X

Kaho Morii https://orcid.org/0000-0002-6752-6061

Junhao Liu https://orcid.org/0000-0002-4774-2998

Á. Sánchez-Monge https://orcid.org/0000-0002-3078-9482

Daniele Galli https://orcid.org/0000-0001-7706-6049

Shantanu Basu https://orcid.org/0000-0003-0855-350X

Patrick M. Koch https://orcid.org/0000-0003-2777-5861

Maria T. Beltrán https://orcid.org/0000-0003-3315-5626

Shanghuo Li https://orcid.org/0000-0003-1275-5251

Henrik Beuther https://orcid.org/0000-0002-1700-090X

Ian W. Stephens https://orcid.org/0000-0003-3017-4418

Fumitaka Nakamura https://orcid.org/0000-0001-5431-2294

Qizhou Zhang https://orcid.org/0000-0003-2384-6589

Wenyu Jiao https://orcid.org/0000-0001-9822-7817

M. Fernández-López https://orcid.org/0000-0001-5811-0454

Jihye Hwang https://orcid.org/0000-0001-7866-2686

Eun Jung Chung https://orcid.org/0000-0003-0014-1527

Kate Pattle https://orcid.org/0000-0002-8557-3582

Luis A. Zapata https://orcid.org/0000-0003-2343-7937

Fengwei Xu https://orcid.org/0000-0001-5950-1932

Fernando A. Olguin https://orcid.org/0000-0002-8250-6827

Ji-hyun Kang https://orcid.org/0000-0001-7379-6263

Janik Karoly https://orcid.org/0000-0001-5996-3600

Chi-Yan Law https://orcid.org/0000-0003-1964-970X

Jia-Wei Wang https://orcid.org/0000-0002-6668-974X

Timea Csengeri https://orcid.org/0000-0002-6018-1371

Xing Lu https://orcid.org/0000-0003-2619-9305

Yu Cheng https://orcid.org/0000-0002-8691-4588

Jongsoo Kim https://orcid.org/0000-0002-1229-0426

Spandan Choudhury https://orcid.org/0000-0002-7497-2713

Huei-Ru Vivien Chen https://orcid.org/0000-0002-9774-1846

Charles L. H. Hull https://orcid.org/0000-0002-8975-7573

## References

* Allen, A., Li, Z.-Y., & Shu, F. H. 2003, ApJ, 599, 363

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

* Bailey, N. D., & Basu, S. 2014, ApJ, 780, 40

* Belloche, A. 2013, in EAS Publications Series, Role and Mechanisms of Angular Momentum Transport During the Formation and Early Evolution of Stars, Evry Schatzman School 2012 62, ed. P. Hennebelle & C. Charbonnel (Les Ulis: EDP Sciences), 25

* Beltrán, M. T., Padovani, M., Girart, J. M., et al. 2019, A&A, 630, A54

* Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25

* CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

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

* Contreras, Y., Sanhueza, P., Jackson, J. M., et al. 2018, ApJ, 861, 14

* Cortés, P. C., Girart, J. M., Sanhueza, P., et al. 2024, arXiv:2406.14663

* Cortés, P. C., Sanhueza, P., Houde, M., et al. 2021, ApJ, 923, 204

* Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2017, A&A, 600, L10

* Csengeri, T., Urquhart, J. S., Schuller, F., et al. 2014, A&A, 565, A75

* Cudlip, W., Furniss, I., King, K. J., & Jennings, R. E. 1982, MNRAS, 200, 1169

* Davis, L. 1951, PhRv, 81, 890

* Estalella, R., Anglada, G., Díaz-Rodríguez, A. K., & Mayen-Gijon, J. M. 2019, A&A, 626, A84

* Fernández-López, M., Sanhueza, P., Zapata, L. A., et al. 2021, ApJ, 913, 29

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

* Girart, J. M., Rao, R., & Marrone, D. P. 2006, Sci, 313, 812

* Hennebelle, P., & Inutsuka, S.-i. 2019, FrASS, 6, 5

* Hildebrand, R. H. 1988, QJRAS, 29, 327

* Hildebrand, R. H., Dragovan, M., & Novak, G. 1984, ApJL, 284, L51

* Huang, B., Girart, J. M., Stephens, I. W., et al. 2024, ApJL, 963, L31

* Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13

* Hwang, J., Kim, J., Pattle, K., et al. 2022, ApJ, 941, 51

* Kirk, J. M., Ward-Thompson, D., Palmeirim, P., et al. 2013, MNRAS, 432, 1424

* Koch, P. M., Tang, Y.-W., Ho, P. T. P., et al. 2018, ApJ, 855, 39

* Kwon, W., Stephens, I. W., Tobin, J. J., et al. 2019, ApJ, 879, 25
Lampton, M., Margon, B., & Bowyer, S. 1976, ApJ, 208, 177

Lazarian, A. 2000, in ASP Conf. Ser. 215, Cosmic Evolution and Galaxy Formation: Structure, Interactions, and Feedback, ed. J. Franco et al. (San Francisco, CA: ASP), 69

Li, H.-B., Jiang, H., Fan, X., Gu, Q., & Zhang, Y. 2017, NatAs, 1, 0158

Li, H.-B., Yuen, K. H., Otto, F., et al. 2015, Natur, 520, 518

Li, S., Sanhueza, P., Beuther, H., et al. 2024, NatAs, 8, 472

Li, Z.-Y., & Shu, F. H. 1996, ApJ, 472, 211

Lin, Y., Csengeri, T., Wyrowski, F., et al. 2019, A&A, 631, A72

Liu, J., Qiu, K., et al. 2022a, ApJ, 925, 30

Liu, J., Zhang, Q., Commerçon, B., et al. 2021, ApJ, 919, 79

Liu, J., Zhang, Q., Qiu, K., et al. 2020, ApJ, 895, 142

Liu, J., Zhang, Q., & Qiu, K. 2022b, FrASS, 9, 943556

Machida, M. N., & Basu, S. 2020, MNRAS, 494, 827

Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS, 477, 2760

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Möller, T., Endres, C., & Schilke, P. 2017, A&A, 598, A7

Mottram, J. C., Hoare, M. G., Davies, B., et al. 2011, ApJL, 730, L33

Mouschovias, T. C., & Ciolek, G. E. 1999, in The Origin of Stars and Planetary Systems (Nato Science Series C), ed. C. J. Lada & N. D. Kylafis, Vol. 540 (Berlin: Springer), 305

Mouschovias, T. C., & Paleologou, E. V. 1979, ApJ, 230, 204

Olguin, F. A., Sanhueza, P., Chen, H.-R. V., et al. 2023, ApJL, 959, L31

Olguin, F. A., Sanhueza, P., Ginsburg, A., et al. 2022, ApJ, 929, 68

Olguin, F. A., Sanhueza, P., Guzmán, A. E., et al. 2021, ApJ, 909, 199

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

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

Padovani, M., Brinch, C., Girart, J. M., et al. 2012, A&A, 543, A16

Padovani, M., & Galli, D. 2011, A&A, 530, A109

Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114

Palau, A., Zhang, Q., Girart, J. M., et al. 2021, ApJ, 912, 159

Pineda, J. E., Offner, S. S. R., Parker, R. J., et al. 2015, Natur, 518, 213

Price, D. J., & Bate, M. R. 2007, MNRAS, 377, 77

Qiu, K., Zhang, Q., Menten, K. M., et al. 2014, ApJL, 794, L18

Rao, R., Girart, J. M., Marrone, D. P., Lai, S.-P., & Schnee, S. 2009, ApJ, 707, 921

Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338

Sanhueza, P., Girart, J. M., Padovani, M., et al. 2021, ApJL, 915, L10

Soler, J. D., & Hennebelle, P. 2017, A&A, 607, A2

Stephens, I. W., Looney, L. W., Kwon, W., et al. 2013, ApJL, 769, L15

Tang, Y.-W., Ho, P. T. P., Koch, P. M., et al. 2009, ApJ, 700, 251

Taniguchi, K., Sanhueza, P., Olguin, F. A., et al. 2023, ApJ, 950, 57

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