Gravity-driven Magnetic Field at ∼1000 au Scales in High-mass Star Formation

High-mass stars dominate the energy input and chemical enrichment of galaxies. Among all the ingredients that have to be considered in their formation, the magnetic field is by far the least explored. Indeed, it is still debated how the magnetic energy compares to other energies in play, namely, turbulence, gravity, and rotation.

Observations of linearly polarized dust emission are currently the best available tool to infer the magnetic field in molecular clouds and denser regions associated with star formation (e.g., Hull & Zhang 2019). Dust polarization observations of high-mass star-forming regions suggest that the magnetic field appears to be dynamically important during the collapse and fragmentation of parsec-scale molecular clumps and the mass assembly of dense cores at scales of 0.01–0.1 pc (Zhang et al. 2014; Hull & Zhang 2019; Cortes et al. 2019). However, it remains unclear whether at smaller scales (core–disk interface, ∼1000 au) the star formation process is magnetically dominated (Girart et al. 2009; Beltrán et al. 2019; Beuther et al. 2020; Cortes et al. 2021; Fernández-López et al. 2021). So far, since polarization observations are scarce at these small scales, it is difficult to evaluate the overall importance of the magnetic field in the high-mass star formation process.

Located at a parallax distance of 2.34 kpc (Xu et al. 2011) with a bolometric luminosity of 1.3 × 10⁴ L_⊙ (Sridharan et al. 2002), the high-mass star-forming region IRAS 18089−1732 is an ideal laboratory to assess the importance of the magnetic field with respect to turbulence, gravity, and rotation. Earlier studies at arcsec resolution show that IRAS 18089−1732 has a deeply embedded hot core (Beuther et al. 2004a, 2004b) and a disk-like rotating structure roughly perpendicular to a molecular outflow. A line-of-sight magnetic field strength of 8.4 and 5.5 mG has been estimated from measurements of Zeeman splitting of the 6.7 GHz CH₃OH maser line using the 100 m Effelsberg telescope (Vlemmings 2008) and the MultiElement Radio Linked Interferometer Network (MERLIN; Dall'Olio et al. 2017), respectively. Early observations of IRAS 18089−1732 with the Submillimeter Array (SMA) show the detection of at least a few independent polarization measurements in dust continuum emission (Beuther et al. 2010). Taking advantage of the superlative capabilities of the Atacama Large Millimeter/submillimeter Array (ALMA), we have observed IRAS 18089−1732 in polarized dust continuum and molecular line emission to better understand the role of the magnetic field in the formation of high-mass stars. This target was observed as part of the Magnetic fields in Massive star-forming Regions (MagMaR) survey that in total contains 30 sources. Details on the survey and source selection will be given in P. Sanhueza et al. (2021, in preparation). Early results on two high-mass star-forming regions, G5.89-0.39 and NGC 6334I(N), are presented in Fernández-López et al. (2021) and Cortes et al. (2021), respectively.

ALMA polarization observations (Project ID: 2017.1.00101.S; PI: Sanhueza) of IRAS 18089−1732 were taken on 2018 September 25. A total of 47 antennas of the 12 m array were used, covering baselines from 15 to 1400 m that resulted in an angular resolution of ∼03. The data set consists of full polarization observations in Band 6 (∼250.486 GHz; 1.2 mm). The correlator setup includes three spectral windows of width 1875 MHz, with a spectral resolution of 1.95 MHz (∼2.4 km s⁻¹), and two spectral windows of width 234 MHz, with a spectral resolution of 0.488 MHz (0.56 km s⁻¹).

Linearly polarized dust continuum emission is detected in the inner ∼8'' of the observed field, the inner one-third of the primary beam (24''), with polarization angles having less than 1% errors. We note, however, that polarization angles on angular scales up to the width of the primary beam have only few percent errors (Hull et al. 2020). Line contamination was removed from the continuum (Stokes I) image following the procedure described in Olguin et al. (2021). Stokes I was self-calibrated in phase and amplitude, while Stokes Q and U were only self-calibrated in phase. Self-calibration solutions were then applied to the spectral cubes.

The continuum imaging was done by independently cleaning each Stokes parameter using the CASA task tclean with Briggs weighting and robust parameter of 1. The resulting images have an angular resolution of 027 × 034 and sensitivities of 175 μJy beam⁻¹ for Stokes I and 31.4 μJy beam⁻¹ for both Stokes Q and U. The polarized intensity image was debiased following Vaillancourt (2006). The peak of the polarized dust emission is 1.4 mJy beam⁻¹. The mean (median) polarization fraction is 5% (4%).

The H¹³CO⁺ line emission was imaged using the automatic masking procedure yclean from Contreras et al. (2018). The CASA task tclean with Briggs weighting and robust parameter of 1 was used, resulting in a noise level of 3.2 mJy beam⁻¹ per 0.56 km s⁻¹channel.

The quasar J1924–2914 was used for the calibration of flux, bandpass, and polarization. The quasar J1832–2039 was used for phase calibration. Data calibration and imaging were performed using CASA 5.1.1 and 5.5.0, respectively.

The ALMA observations at 1.2 mm with a spatial resolution of 700 au (03) allow us to observe the internal structure of IRAS 18089−1732 in great detail (Figure 1). The dust continuum emission peaks at the position of the previously reported hot core (Beuther et al. 2004a, 2004b; Zapata et al. 2006). The overall emission is asymmetric, with an extension toward the northwest of the brightest 1.2 mm peak. In addition, there are spiral-like streamers associated with the central hot core having a counterclockwise orientation (mimicking a whirlpool). The magnetic field projected in the plane of the sky also traces spiral-like features connected to the central hot core, roughly following the dusty spirals and making more evident the whirlpool shape (Figure 1).

Zoom In Zoom Out Reset image size

Along with the continuum emission, there are myriad molecular line transitions detected with different excitation conditions (i.e., tracing different temperature and density regimes). Among these, here we focus on H¹³CO⁺ (J = 3−2) and leave the tracers of the inner hot core for a forthcoming work (V. Chen et al. 2021, in preparation). The emission from the H¹³CO⁺ (J = 3−2) line shows a spatial distribution coincident with that of the dust emission, as can be seen in the integrated intensity map (Figure 2(a)). This molecule is a good tracer of the relatively cold (E_u = 25 K), dense material in the spiral-like and filamentary structures found in IRAS 18089−1732. The H¹³CO⁺ spectra show two clear velocity components separated by ∼4 km s⁻¹ that can be traced continuously from the regions in which they overlap (presenting double-line profiles) to regions where only a single component is detected. To properly trace the velocity field of both components, a simultaneous Gaussian fitting of both velocity components was performed. Figure 2(b) shows the velocity field traced by the blueshifted gas component (with respect to the systemic cloud velocity v_LSR ≈ 33.0 km s⁻¹), while Figure 2(c) shows the velocity field traced by the redshifted gas component. Near the center of the hot core, the H¹³CO⁺ profile becomes complex, exhibiting extremely broad line widths, absorption features, and blending with hot core lines, precluding a successful Gaussian fitting.

Zoom In Zoom Out Reset image size

A prominent spiral-like feature in the blueshifted gas component that follows the magnetic field morphology is seen in Figure 2(b). Although there is significant substructure in the velocity pattern along this spiral-like filament, it has an overall velocity difference of ∼1.4 km s⁻¹ over a length of 7.3'' (17000 au), resulting in a 17 km s⁻¹ pc⁻¹ velocity gradient. The redshifted gas component shows two filaments that coincide with the area emitting polarized emission (Figure 2(c)), one extending to the north and the other one to the south of the central hot core. The northern filament starts at 33 (7700 au) from the hot core with blueshifted velocities of ∼34.1 km s⁻¹ and then connects to the red side of the rotating central hot core at ∼36.1 km s⁻¹ (velocity gradient of 53.6 km s⁻¹ pc⁻¹). The southern filament extends over 29 (6800 au) with extreme velocities of ∼34.2 and ∼36 km s⁻¹, resulting in a velocity gradient of 54.6 km s⁻¹ pc⁻¹. In the position–position–velocity (PPV) space displayed in Figure 3, the distribution of the gas and the main structures, spiral and northern/southern filaments, can be seen. Figure 3(a) simultaneously displays both the blueshifted and redshifted velocity components and is comparable to displaying Figures 2(b) and (c) in a single image. Figures 3(b) and (c) show the three main structures from different angles. Figure 3(d) corresponds to a position–velocity (PV) diagram. Figure 3(d) shows, at the positions indicated by the black arrows, how the velocity increases closer to the central hot core, especially in the southern and northern filaments. Such an increase in velocity is typical of gas accelerating close to the central object and is a sign of infall (Tobin et al. 2012).

Zoom In Zoom Out Reset image size

4.1. Dust Continuum Emission

4.2. Line Emission and Accretion Flows

4.3. Modeling the Magnetic Field

To model the polarization pattern observed in IRAS 18089−1732, we used the DustPol module contained in the ARTIST package (Padovani et al. 2012) and followed the same procedure described in Beltrán et al. (2019). We then perform synthetic ALMA observations of the models and compare the resulting polarization position angles with the observations. DustPol creates a set of FITS images related to the Stokes parameters (I, Q, U) that is directly used as an input for the CASA simobserve/simanalyze tasks, which use the same antenna configuration of the observing runs.

The model used for the magnetic field configuration is an axisymmetric singular toroid threaded by an hourglass-shaped poloidal field (Li & Shu 1996; Padovani & Galli 2011) with an added toroidal component to mimic the effect of rotation as in Padovani et al. (2013). More details on the model are given in Appendix C. The main purpose of the modeling is to determine basic properties of the magnetic field (mass-to-magnetic-flux ratio, toroidal-to-poloidal ratio, inclination) to be compared with observations, and to provide the configuration of the mean field needed to perform the Davis–Chandrasekhar–Fermi (DCF; Davis 1951; Chandrasekhar & Fermi 1953) analysis of the polarization angle residuals.

We performed a χ² test for all the combinations of the parameters λ, b₀, and i and found that the set λ = 8.38, ${b}_{0}\,={0.30}_{-0.18}^{+0.28}$, and $i=-{5}_{-15}^{+25}$° gives the lowest reduced χ² value (${\bar{\chi }}^{2}=4.13$). The magnetic field is dominated by a poloidal component, but with an important contribution from the toroidal component, whose magnitude is 30% of the poloidal component. This significant contribution from the toroidal component indicates that rotation is affecting the magnetic field. Figure 4, shows the comparison between the observed and modeled polarization angles. The excellent agreement is shown quantitatively in the inset of the same figure. The latter illustrates the distribution of the polarization angle residuals, Δψ = ψ_obs − ψ_mod, defined as the difference between the observed (ψ_obs) and modeled (ψ_mod) polarization angles, whose Gaussian fit gives a mean value of 〈Δψ〉 = 115 and a standard deviation σ_ψ = 1809. The latter is needed for the DCF method. The λ value of 8.38 of the best model characterizes the mass inside a flux tube, while observationally the mass is estimated assuming a spherical source. To correct for this, λ must be divided by a correction factor of 2.32 (Li & Shu 1996), leading to an effective λ of 3.61. Despite the relative simplicity of the model, this value of λ is very close to that estimated from observations (see Section 4.4). It should be stressed that cloud models with smaller values of λ cannot be completely ruled out. For example, a model with λ = 2.66 also provides an acceptable fit to the data, resulting in ${b}_{0}={0.40}_{-0.21}^{+0.50}$ and $i={0}_{-22}^{+22}$° with a slightly worse value of ${\bar{\chi }}^{2}=4.65$ (for smaller values of λ, the ${\bar{\chi }}^{2}$ becomes progressively larger). We note that the magnetic field strength determined in the following section would not change adopting the λ = 2.66 model because its σ_ψ = 1852 practically has the same value (1809) as in the preferred λ = 8.38 model.

Zoom In Zoom Out Reset image size

4.4. Analysis of the Energy Balance

ALMA observations of the high-mass star-forming region IRAS 18089−1732 have revealed that the dense molecular envelope surrounding the high-mass star has a complex spiral pattern at the 0.003–0.1 pc scales. This spiral-like morphology is seen in the gas and dust (traced by the H¹³CO⁺ (J = 3−2) line and the 1.2 mm continuum emission, respectively), as well as in the magnetic field (traced by the linearly polarized 1.2 mm dust emission). At the observed size scales, gravitational infall clearly dominates over the support from the magnetic field, turbulence, and rotation, resulting in the feeding of the inner dense, hot circumstellar disk-like structure with a high accretion rate of 0.9–2.5 × 10⁻⁴ M_⊙ yr⁻¹. We show that high-mass star formation can occur in weakly magnetized environments and that gravity is shaping the immediate surrounding around the high-mass star. The spiral magnetic field indicates that angular momentum is high enough to twist the field lines, as supported by the model and the energy analysis. With these observations and consistent with previous works in other high-mass star-forming regions (e.g., Koch et al. 2014, 2018), we suggest that the importance of the magnetic field in the process of high-mass star formation depends on the size scales traced and the evolutionary stage of the observed region.

P.S. and B.W. were partially supported by a Grant-in-Aid for Scientific Research (KAKENHI number 18H01259) of the Japan Society for the Promotion of Science (JSPS). J.M.G. is supported by the Spanish grant AYA2017-84390-C2-R (AEI/FEDER, UE). C.L.H.H. acknowledges the support of the NAOJ Fellowship and JSPS KAKENHI grants 18K13586 and 20K14527. J.M.J.'s research was conducted in part at the SOFIA Science Center, which is operated by the Universities Space Research Association under contract NNA17BF53C with the National Aeronautics and Space Administration. K.T. was supported by JSPS KAKENHI grant No. 20H05645. 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. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00101.S. 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.

Facility: ALMA. -

Software: CASA (v5.1.1, 5.5; McMullin et al. 2007).

The total gas mass can be calculated from the dust emission, in the optically thin limit, as

$$\begin{eqnarray}&&M={\mathbb{R}}\,\displaystyle \frac{{F}_{\nu }{D}^{2}}{{\kappa }_{\nu }{B}_{\nu }(T)}\,,\end{eqnarray} \tag{ A1 }$$

where F_ν is the source flux density, ${\mathbb{R}}$ is the gas-to-dust mass ratio, D is the distance to the source, κ_ν is the dust opacity per gram of dust, and B_ν is the Planck function at the dust temperature T. Assuming a dust-to-gas mass ratio of 100, dust opacity of 1.03 cm² g⁻¹ (interpolated to 1.2 mm assuming β = 1.6; Ossenkopf & Henning 1994), and a temperature of 30 K (Lu et al. 2014), the computed total mass is 75 M_⊙. The number density, n(H₂) = M/(Volume × ${\mu }_{{{\rm{H}}}_{2}}\,{m}_{{\rm{H}}}$) with ${\mu }_{{{\rm{H}}}_{2}}$ the molecular weight per hydrogen molecule and m_H the hydrogen mass, and the surface density, Σ = M/(π r²), can be calculated assuming a spherical core. We define an effective radius (r_eff) determined by the area (A) emitting above 10σ as r_eff = (A/π)^1/2, resulting in r_eff = 236 (equal to 0.027 pc or 5500 au). Assuming ${\mu }_{{{\rm{H}}}_{2}}=2.8$, we obtain an n(H₂) of 1.3 × 10⁷ cm⁻³ and Σ = 6.9 g cm⁻², values characteristic of cores forming high-mass stars. The free-fall time is determined from

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}=\sqrt{\displaystyle \frac{3\pi }{32G\rho }}\,,\end{eqnarray} \tag{ A2 }$$

where G is the gravitational constant and ρ is the mass density, n(H₂) × ${\mu }_{{{\rm{H}}}_{2}}\,{m}_{{\rm{H}}}$. With ρ equal to 6.3 × 10⁻¹⁷ g cm⁻³, the free-fall time is 8.4 × 10³ yr.

We estimated the total gas mass in the spiral-like feature and filaments from the H¹³CO⁺ emission following a standard procedure (Sanhueza et al. 2012) as follows. First, the column density for a linear, rigid rotor in the optically thin regime, assuming a filling factor of unity, can be calculated from

$$\begin{eqnarray}\begin{array}{rcl}N & = & \displaystyle \frac{3k}{8{\pi }^{3}{B}_{\mathrm{rot}}{\mu }^{2}}\displaystyle \frac{({T}_{\mathrm{ex}}+{{hB}}_{\mathrm{rot}}/3{k}_{{\rm{B}}})}{(J+1)}\displaystyle \frac{\exp ({E}_{{}_{J}}/{{kT}}_{\mathrm{ex}})}{[1-\exp (-h\nu /{{kT}}_{\mathrm{ex}})]}\\ & & \times \,\displaystyle \frac{1}{[J({T}_{\mathrm{ex}})-J({T}_{\mathrm{bg}})]}\displaystyle \int {T}_{{\rm{b}}}\,{dv}\,,\end{array}\end{eqnarray} \tag{ B1 }$$

where k_B is the Boltzmann constant, h is the Planck constant, T_ex is the excitation temperature (30 K; Lu et al. 2014), ν is the transition frequency (260.255339 GHz), μ is the permanent dipole moment of the molecule (3.89 D), J is the rotational quantum number of the lower state, E_J = hB_rot J(J + 1) is the energy in the level J, B_rot is the rotational constant of the molecule (43.377302 GHz), T_b is the brightness temperature, T_bg the background temperature, and J(T) is defined as

$$\begin{eqnarray}&&J(T)=\displaystyle \frac{h\nu }{k}\displaystyle \frac{1}{{e}^{h\nu /{kT}}-1}\,.\end{eqnarray} \tag{ B2 }$$

The column density is then converted into mass using

$$\begin{eqnarray}&&M={\left[\displaystyle \frac{{{\rm{H}}}^{13}{\mathrm{CO}}^{+}}{{{\rm{H}}}_{2}}\right]}^{-1}{\mu }_{{{\rm{H}}}_{2}}A\sum N({{\rm{H}}}^{13}{\mathrm{CO}}^{+})\,,\end{eqnarray} \tag{ B3 }$$

where A is the size of the emitting area in an individual position (pixel of 0.05''), [H¹³CO⁺/H₂] is the H¹³CO⁺ to molecular hydrogen abundance ratio, and the sum is over all the observed positions. We have assumed an [H¹³CO⁺/H₂] abundance ratio typical of high-mass star-forming regions equal to 1.28 × 10⁻¹⁰ (Hoq et al. 2013).

The idealized model used represents a cloud in magnetostatic equilibrium with density and magnetic field strength decreasing with radius as r⁻² and r⁻¹, respectively. The field lines are spirals wrapping around nested hourglass-shaped magnetic flux tubes, and have a kink in the equatorial plane where the toroidal component changes sign. The parameters of the model are (i) the mass-to-magnetic-flux ratio normalized to its critical value, λ, given by

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{(M/{{\rm{\Phi }}}_{B})}{{\left(M/{{\rm{\Phi }}}_{B}\right)}_{\mathrm{cr}}}=2\pi {G}^{1/2}\displaystyle \frac{M}{{{\rm{\Phi }}}_{B}}\,,\end{eqnarray} \tag{ C1 }$$

where G is the gravitational constant, M is the core mass, and Φ_B is the magnetic flux; (ii) the ratio of the strength of the toroidal and poloidal magnetic field components in the midplane, b₀; and (iii) the inclination of the toroid with respect to the plane of the sky, i (i = 0 for face-on). The latter is assumed positive if the magnetic field in the north side is directed toward the observer. In principle the orientation of the projection of the magnetic axis on the plane of the sky, φ, should also be considered as a further parameter, but as we will show, the best-fit model indicates that the equatorial plane of the cloud is practically in the plane of the sky (i = −5°), so that the model turns out to be insensitive to φ.

We considered four values for λ. The extreme cases include a nearly spherical density profile with a weak magnetic field (λ = 16.2) and a flat density profile with a strong magnetic field (λ = 1.63). The two intermediate cases correspond to λ = 2.66 and 8.38. For the magnetic field configuration, we have considered the range of b₀ from the pure poloidal case (b₀ = 0) to the case where the strength of the toroidal component is twice the poloidal one (b₀ = 2). Finally, we have taken inclinations with respect to the plane of the sky between −90° and 90°.

The lowest χ² value, 4.13, is found for the set λ = 8.38, ${b}_{0}={0.30}_{-0.18}^{+0.28}$, and $i=-{5}_{-15}^{+25}$° (uncertainties on b₀ and i have been estimated using the method of Lampton et al. 1976). The best-fit model can be seen in Figure 4. For completeness, Figure 5 shows the map of the ${\bar{\chi }}^{2}$ distribution as a function of the explored range of b₀ and i. Superimposed on this map are the isocontours of the average value of the polarization angle residuals, 〈Δψ〉 (left panel), and the skewness (also known as moment 3) parameter of their distribution (right panel). The mass-to-magnetic-flux ratio that corresponds to the minimum ${\bar{\chi }}^{2}$ also gives the distribution of Δψ with the smallest degree of skewness (equal to 0.12, namely, the distribution is only slightly positively skewed, which means that the peak of the distribution is shifted toward negative values as seen in the inset of Figure 4).

Zoom In Zoom Out Reset image size

The magnitude of the magnetic field was estimated using the DCF method, following the procedure of Beltrán et al. (2019)

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}=\xi \displaystyle \frac{{\sigma }_{\mathrm{los}}}{\delta {\psi }_{\mathrm{int}}}\sqrt{4\pi \rho }\,,\end{eqnarray} \tag{ D1 }$$

where ξ = 0.5 is a correction factor derived from turbulent cloud simulations (Ostriker et al. 2001) and σ_los is the line-of-sight velocity dispersion. The DCF method assumes that σ_los is from turbulence, and indeed, the observed velocity dispersion is 2.6 times the thermal line width. The average velocity dispersion (σ_los) for both H¹³CO⁺ velocity components above the 10σ area, in the continuum image, is 0.75 km s⁻¹ (both blueshifted and redshifted components have practically the same line width). The density used in the estimation of the magnetic field strength is 6.3 × 10⁻¹⁷ g cm⁻³. δ ψ_int is the intrinsic angle dispersion given by δ ψ_int = ${({\sigma }_{\psi }^{2}-\delta {\psi }_{\mathrm{obs}}^{2})}^{1/2}$. σ_ψ of 181 is the standard deviation of the polarization angle residuals obtained from the magnetic field modeling (see inset in Figure 4, top), and δ ψ_obs of 60 is the mean angle uncertainty in the whole debiased polarization angle image. The derived δ ψ_int of 171 results in a magnetic field strength of B_pos = 3.5 mG. The Alfvén speed, given by ${\sigma }_{A}=B/\sqrt{4\pi \rho }$, is 1.26 km s⁻¹.

We calculate the mass-to-magnetic-flux ratio normalized to the critical value, λ, using Equation (C1) and the fact that Φ_B = BπR². Adopting R = r_eff, we obtain a value of λ equal to 3.2.

The ratio of the turbulent to magnetic energy, β_t , is given by

$$\begin{eqnarray}&&{\beta }_{t}=3{\left(\displaystyle \frac{{\sigma }_{\mathrm{los}}}{{\sigma }_{A}}\right)}^{2}\,.\end{eqnarray} \tag{ D2 }$$

The calculated β_t is 1.07, which indicates that turbulence and the magnetic field play a comparable role in the energy budget of the system.

We assess the importance of the rotational energy with respect to gravity and the magnetic field as follows. First, we take the ratio between the gravitational energy, E_G,

$$\begin{eqnarray}&&{E}_{{\rm{G}}}=-\displaystyle \frac{{{GM}}^{2}}{R}\left(\displaystyle \frac{3-\alpha }{5-2\alpha }\right)\,,\end{eqnarray} \tag{ D3 }$$

and the rotational energy, E_rot,

$$\begin{eqnarray}&&{E}_{\mathrm{rot}}=\displaystyle \frac{1}{2}I{{\rm{\Omega }}}^{2}=\displaystyle \frac{1}{3}{{Mv}}_{\mathrm{rot}}^{2}\left(\displaystyle \frac{3-\alpha }{5-\alpha }\right)\,,\end{eqnarray} \tag{ D4 }$$

where the moment of inertia (I) of a sphere has been assumed, the angular velocity is given by Ω = v_rot/R, v_rot is the rotational velocity, and α corresponds to the exponent in a density profile of the form ρ(R) ∝ R^−α (where α = 0 implies a uniform density profile and α = 2 a centrally peaked density profile; Belloche 2013). The maximum velocity gradient observed is of 2 km s⁻¹. If we assume as the upper limit that this gradient is purely produced by rotation, we obtain a E_rot/E_G ratio of 0.11 for a uniform density profile (α = 0) and 0.037 for a centrally peak density profile (α = 2).

To evaluate the relative importance between rotation and the magnetic field, we follow two approaches: the rotation-to-magnetic-field-energy ratio and the prescription from Machida et al. (2005). First, the energy in the magnetic field (E_B) is given by

$$\begin{eqnarray}&&{E}_{{\rm{B}}}=\displaystyle \frac{1}{8\pi }{B}^{2}V=\displaystyle \frac{1}{6}{B}^{2}{R}^{3}\,,\end{eqnarray} \tag{ D5 }$$

where V is the volume, assumed here to be a sphere. Taking into account the same consideration from above (velocity gradient entirely produced by rotation), the E_rot/E_B ratio is 1.0 for a uniform density distribution and 0.56 for a centrally peaked density distribution. These results suggest that the magnetic field is slightly more important than rotation, depending on the exact density profile.

Second, for assessing whether the magnetic field dominates over rotation during the collapse of clouds, Machida et al. (2005) suggest that the ratio of the angular velocity to the magnetic flux density can be used. According to the following relation, if the observed ratio is

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{{\rm{\Omega }}}{B}\right)}_{\mathrm{obs}}\gt 0.39\displaystyle \frac{\sqrt{G}}{{\sigma }_{\mathrm{th}}}\\ =\,1.69\times {10}^{-7}\times {\left(\displaystyle \frac{{\sigma }_{\mathrm{th}}}{0.19\,\mathrm{km}\ \,{{\rm{s}}}^{-1}}\right)}^{-1}{\mathrm{yr}}^{-1}\,\mu {{\rm{G}}}^{-1}\,,\end{array}\end{eqnarray} \tag{ D6 }$$

rotation dominates over the magnetic field; otherwise, the magnetic field dominates over rotation. Evaluating the thermal velocity dispersion (isothermal sound speed), ${\sigma }_{\mathrm{th}}={\left(\tfrac{{k}_{B}T}{\mu {m}_{H}}\right)}^{1/2}$, for a temperature of 30 K and using the mean molecular weight per free particle μ = 2.37, we obtain the right-hand side of Equation (D6) equal to 9.8 × 10⁻⁸ yr⁻¹ μG⁻¹. Evaluating the observed ratio, we obtain 2.2 × 10⁻⁸ yr⁻¹ μG⁻¹. The observed ratio is lower than the right-hand side of Equation (D6) by a factor 4.5, suggesting that from this analysis the magnetic field dominates over rotation. Based on both methods employed to assess the importance of rotation with respect to the magnetic field, we conclude that the magnetic field is slightly more important than rotation.

The dynamical state of cores is generally evaluated by using the virial theorem. The ratio between the virial mass, M_vir, and the total mass defines the virial parameter, α_vir. A virial parameter of unity implies equilibrium, α_vir < 1 implies gravitational collapse, and α_vir > 1 means the core is expanding and it will disperse. In its simplest form, the virial analysis only includes gravity and the kinetic energy (turbulence and thermal energy) as follows:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{vir}}=\displaystyle \frac{{M}_{\mathrm{vir}}}{M}=3\left(\displaystyle \frac{5-2\alpha }{3-\alpha }\right)\displaystyle \frac{R{\sigma }_{\mathrm{los}}^{2}}{{GM}}\,,\end{eqnarray} \tag{ D7 }$$

resulting in α_vir = 0.24 for a uniform density profile (α = 0) and 0.14 for a centrally peaked density profile (α = 2). Both low α_vir values, in the absence of other energies, indicate that the core is collapsing dominated by gravity.

The virial parameter considering the magnetic field, is written as follows (e.g., Sanhueza et al. 2017):

$$\begin{eqnarray}&&{\alpha }_{\mathrm{vir},{\rm{B}}}=\displaystyle \frac{{M}_{\mathrm{vir}}}{M}=3\left(\displaystyle \frac{5-2\alpha }{3-\alpha }\right)\displaystyle \frac{R}{{GM}}\left({\sigma }_{\mathrm{los}}^{2}+\displaystyle \frac{{\sigma }_{A}^{2}}{6}\right),\end{eqnarray} \tag{ D8 }$$

which results in α_vir,B of 0.35 for a uniform density profile and 0.21 for a centrally peaked density profile. Including the magnetic field, gravity continues to be dominant and the core is collapsing.

We derive the contribution of rotation in the virial equation as below. The virial theorem can be written as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}I}{{{dt}}^{2}} & = & 2{E}_{{\rm{k}}}+{E}_{{\rm{G}}}+{E}_{{\rm{B}}}\,,\\ & = & 2({E}_{{\rm{n}}-\mathrm{rot}}+{E}_{\mathrm{rot}})+{E}_{{\rm{G}}}+{E}_{{\rm{B}}}\,,\end{array}\end{eqnarray} \tag{ D9 }$$

in which E_k is the kinetic energy that can be separated in a rotational part (E_rot) and nonrotational part (E_n−rot; McKee & Zweibel 1992). The nonrotational part includes thermal and turbulent energies as follows:

$$\begin{eqnarray}&&{E}_{{\rm{n}}-\mathrm{rot}}=\displaystyle \frac{3}{2}M{\sigma }_{\mathrm{los}}^{2}\,,\end{eqnarray} \tag{ D10 }$$

where the observed velocity dispersion is ${\sigma }_{\mathrm{los}}^{2}={\sigma }_{\mathrm{th}}^{2}+{\sigma }_{\mathrm{tur}}^{2}$ and σ_tur is the turbulent component. E_n−rot is frequently assumed as E_k when rotation is neglected. In virial equilibrium, the moment of inertia does not vary over time and the left-hand side of Equation (D9) becomes zero. Solving Equation (D9) for the mass, one can find the virial mass and the virial parameter can be calculated. If rotation and magnetic energies are ignored, Equation (D7) is recovered. If only rotation is ignored, Equation (D8) is found. Solving for all energies in play, we obtain

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{\mathrm{vir},{\rm{B}},\mathrm{rot}}=3\left(\displaystyle \frac{5-2\alpha }{3-\alpha }\right)\displaystyle \frac{R}{{GM}}\\ \,\times \,\,\left({\sigma }_{\mathrm{los}}^{2}+\displaystyle \frac{{\sigma }_{A}^{2}}{6}+\displaystyle \frac{2}{9}\left(\displaystyle \frac{3-\alpha }{5-\alpha }\right){v}_{\mathrm{rot}}^{2}\right).\end{array}\end{eqnarray} \tag{ D11 }$$

For a uniform density profile and for a centrally peaked density distribution, we obtain virial parameters of 0.57 and 0.28, respectively. This implies that at the scales we probe with our ALMA observations, even when both the magnetic field and rotation are considered, gravity is still the dominant dynamical force (independently of the density profile assumed).