Revealing Gravitational Collapse in the Serpens G3–G6 Molecular Cloud Using Velocity Gradients

We obtained a new ¹³CO (J = 10) fully sampled map of Serpens using the IRAM 30 m telescope (Carter et al. 2012). The observations were obtained in July 2020 using 16 hr of telescope time under average summer weather conditions (6 mm median water vapor). We covered a field of view (FoV) of 10' × 40'. The ¹³CO (J = 10) emission was observed using the EMIR receiver and the VESPA spectrometer using a bandwidth of 60 MHz at 0.092 MHz resolution (≈0.212 km s⁻¹). The half-power beamwidth (HPBW) at 110.201354 GHz is ≈235.

We used the on-the-fly scanning strategy with a dump time of 0.7 s and a scanning speed of 11''s^–1 to ensure a sampling of three dumps per beam along the scanning direction, with the scanning direction reversed after each raster line (i.e., zigzag scanning mode). We covered the full FoV (≈0.11 square degrees) with 15 rectangular tiles (of 15 OTF scans each) along the decl. direction and 20 rectangular tiles (of 5 OTF scans each) along the R.A. direction, followed by a calibration measurement. The reference position (REF) was observed for 10 s after each raster line following the pattern REFOTFOTFOTF–OTF–OTF–REF along the decl. direction and the pattern REFOTFOTFOTF–REF along the R.A. direction. Each tile is of approximately 8'' × 40'' size. In total, we employed about 13 minutes per tile.

Data reduction was carried out using the GILDAS1/CLASS software. ⁵ The data were first calibrated to the T_A scale and were then corrected for atmospheric absorption and spillover losses using the chopper-wheel method (Penzias & Burrus 1973). A polynomic baseline of second order was subtracted from each spectrum, avoiding velocities with molecular emission. The spectra were then gridded into a data cube through convolution with a Gaussian kernel of FWHM ∼1/3 of the IRAM 30 m telescope beamwidth at the rest line frequency. The typical (1σ) rms noise level achieved in the map is 0.33 K per 212 m s⁻¹ velocity channel. A large table of the individual spectra was made, and the spectra were finally combined to obtain a regularly gridded position–position–velocity data cube, setting the pixel size to 5''. Here we convert the measured antenna temperature T_A to brightness temperature T_b through T_b = (F_eff/B_eff)T_A, where F_eff = 0.95 is the forward efficiency of the IRAM 30 m telescope and B_eff = 0.79 is the main beam efficiency at 110.201354 GHz (Pety et al. 2017).

The ¹²CO (2–1) and ¹³CO (2–1) emission lines were observed with the Heinrich Hertz Submillimeter Telescope (Burleigh et al. 2013), while the H₂ column density data were obtained from the Herschel Gould Belt Survey (André et al. 2010). Each line was measured with 256 filters of 0.25 MHz bandwidth, giving a total spectral coverage of 41 km s⁻¹ at a resolution of 0.33 km s⁻¹.

The angular resolution of the emission lines is 38'' (0.04 pc) with a sensitivity of 0.12 K rms noise per pixel in one spectral channel (Burleigh et al. 2013). The radial velocity of the bulk of the emission ranges from about −1 to +18 km s⁻¹ for ¹²CO (2–1) and from +2 to +13 km s⁻¹ for ¹³CO (2–1) (Burleigh et al. 2013). We select the emissions within these ranges for our analysis.

To trace the magnetic field orientation in the plane of the sky (POS), we use the Planck 353 GHz polarized dust signal data from the Planck 3rd Public Data Release (DR3) 2018 of the High Frequency Instrument (Planck Collaboration et al. 2020). ⁶

The Planck observations define the polarization angle ϕ and polarization fraction p through Stokes parameter maps I, Q, and U:

$$\begin{eqnarray}\begin{array}{rcl}\phi & = & \displaystyle \frac{1}{2}\arctan (-U,Q)\\ p & = & \sqrt{{Q}^{2}+{U}^{2}}/I,\end{array}\end{eqnarray} \tag{ 1 }$$

where −U converts the angle from HEALPix convention to IAU convention, and the two-argument function $\arctan$ is used to account for the π periodicity. To increase the signal-to-noise ratio, we smooth all maps from nominal angular resolution $5^{\prime}$ up to a resolution of $10^{\prime}$ using a Gaussian kernel. The magnetic field angle is inferred from ϕ_B = ϕ + π/2.

VGT (González-Casanova & Lazarian 2017; Yuen & Lazarian 2017a; Hu et al. 2018; Lazarian & Yuen 2018a) is the main analysis tool in the work. It is theoretically rooted in the advanced MHD turbulence theory (Goldreich & Sridhar 1995, henceforth GS95), including the concept of fast turbulent reconnection theory (Lazarian & Vishniac 1999, henceforth LV99). These theories explained that turbulent eddies are anisotropic (see GS95) so that their semimajor axis is elongating along the local magnetic fields (see LV99). This was numerically demonstrated by Cho & Vishniac (2000) and Maron & Goldreich (2001). In particular, LV99 derived the anisotropy relation for the eddies in the local reference frame:

$$\begin{eqnarray}&&{l}_{\parallel }\simeq {L}_{\mathrm{inj}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{\tfrac{2}{3}}{M}_{{\rm{A}}}^{-4/3},\end{eqnarray} \tag{ 2 }$$

where l_⊥ and l_∥ are the perpendicular and parallel sizes of eddies with respect to the local magnetic field. M_A is the Alfvén Mach number, and L_inj is the injection scale of turbulence. The corresponding scaling of velocity fluctuation v_l at scale l is then

$$\begin{eqnarray}&&{v}_{l}\simeq {v}_{\mathrm{inj}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{\tfrac{1}{3}}{M}_{{\rm{A}}}^{1/3},\end{eqnarray} \tag{ 3 }$$

where v_inj is the injection velocity. Explicitly, as the anisotropic relation indicates l_⊥ ≪ l_∥, the velocity gradients scale as (Yuen & Lazarian 2020)

$$\begin{eqnarray}&&{\rm{\nabla }}{v}_{l}\propto \displaystyle \frac{{v}_{l}}{{l}_{\perp }}\simeq \displaystyle \frac{{v}_{\mathrm{inj}}}{{L}_{\mathrm{inj}}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{-\tfrac{2}{3}}{M}_{{\rm{A}}}^{\tfrac{1}{3}},\end{eqnarray} \tag{ 4 }$$

which means that the smallest resolved eddies induce the largest gradients. The gradients are perpendicular to the local magnetic fields.

As we discussed above, the velocity fluctuation, i.e., the velocity variance, is related to the eddy's size along the line of sight (LOS). Large velocity variance eddies give the most significant contribution to the velocity gradients. To get the largest velocity variance, Hu et al. (2018) proposed the principal component analysis (PCA) to preprocess the spectroscopic PPV cubes.

PCA treats the observed brightness temperature T_b(x, y, v) in the PPV cube as the probability density function of three random variables x, y, v. By splitting the PPV cube into n_v velocity channels along the LOS, the covariance matrix S(v_i , v_j ) and its corresponding eigenvalue equation are defined as

$$\begin{eqnarray}\begin{array}{rcl}S({v}_{i},{v}_{j}) & \propto & \displaystyle \int {{dxdyT}}_{{\rm{b}}}(x,y,{v}_{i}){T}_{{\rm{b}}}(x,y,{v}_{j})\\ & & -\displaystyle \int {{dxdyT}}_{{\rm{b}}}(x,y,{v}_{i})\displaystyle \int {{dxdyT}}_{{\rm{b}}}(x,y,{v}_{j}),\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\boldsymbol{S}}\cdot {\boldsymbol{u}}=\lambda {\boldsymbol{u}},\end{eqnarray} \tag{ 6 }$$

where i, j = 1, 2, ..., n_v and λ are the eigenvalues associated with the eigenvector u . The equation gives a total of n_v eigenvalues. Each eigenvalue λ_i is proportional to the squared velocity variance v_i ² along the LOS. The velocity channel corresponding to the largest eigenvalue, therefore, has the most significant contribution to the velocity gradients.

In the frame of VGT, the PPV cube can be preprocessed by PCA in two different ways. The first one is treating the eigenvalues λ as the weighting coefficients for each velocity channel (Hu et al. 2018). By integrating the weighted velocity channels along the LOS, we obtained the velocity centroid map C(x, y):

$$\begin{eqnarray}&&C(x,y)=\displaystyle \frac{\int {{dvT}}_{{\rm{b}}}(x,y,v)\cdot v\cdot \lambda (v)}{\int {{dvT}}_{{\rm{b}}}(x,y,v)}.\end{eqnarray} \tag{ 7 }$$

Large eigenvalues enhance the significance of their corresponding velocity channels, while small eigenvalues give a suppression.

Alternatively, instead of weighting velocity channels, the PPV cube can be projected into a new orthogonal basis constructed by n_v eigenvectors (Hu et al. 2020e). Its corresponding eigenvalue constrains the length of each eigenvector. Because the new axes are oriented along the direction of maximum variance, the PPV cube on a new basis exhibits the largest velocity variance. The projection is operated by the weighting channel T_b(x, y, v_j ) with the corresponding eigenvector element u_ij , in which the corresponding eigenchannel I_i (x, y) is

$$\begin{eqnarray}&&{I}_{i}(x,y)={\displaystyle \sum _{j}}^{{n}_{v}}{u}_{{ij}}\cdot {T}_{{\rm{b}}}(x,y,{v}_{j}).\end{eqnarray} \tag{ 8 }$$

This step produces a total of n_v eigenchannels in the eigenvector space. Generally, both the weighting and projecting approaches are figuring out the maximum velocity variance, which is the most important in velocity gradients' calculation. In this work, we adopt the projection approach to construct the pseudo-Stokes parameters; see the discussion below.

Note here that the channel number n_v should be sufficiently large so that the channel width Δv satisfies ${\rm{\Delta }}v\lt \sqrt{\delta ({v}^{2})}$, where $\sqrt{\delta ({v}^{2})}$ is the 1D velocity dispersion. We call the channel, which satisfies the criteria, a thin velocity channel. Due to the velocity caustic effect, which comes from the nonlinear mapping from real space to PPV space, the thin velocity channel's intensity fluctuation is mainly induced by velocity fluctuations instead of density fluctuations (Lazarian & Pogosyan 2000). The thin velocity channels, therefore, record velocity information, which is used to calculate velocity gradients.

The PPV cubes are preprocessed by PCA to extract the most crucial velocity components resulting in n_v eigenchannels. Each eigenchannel is convolved with 3 × 3 Sobel kernels ⁷ G_x and G_y to calculate pixelized gradient map ${\psi }_{g}^{i}(x,y)$:

$$\begin{eqnarray}\begin{array}{rcl}{\bigtriangleup }_{x}{I}_{i}(x,y) & = & {G}_{x}\ast {I}_{i}(x,y)\\ {\bigtriangleup }_{y}{I}_{i}(x,y) & = & {G}_{y}\ast {I}_{i}(x,y)\\ {\psi }_{g}^{i}(x,y) & = & {\tan }^{-1}\left(\displaystyle \frac{{\bigtriangleup }_{y}{I}_{i}(x,y)}{{\bigtriangleup }_{x}{I}_{i}(x,y)}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where △_x I_i (x, y) and △_y I_i (x, y) are the x and y components of the gradient, respectively, and ∗ denotes the convolution.

Note that the turbulent eddy discussed above is a statistical concept. A single gradient, therefore, does not necessarily correlate with the magnetic field. The orthogonal relative orientation between velocity gradients and the magnetic field appears only when the sampling is statistically sufficient. The statistical sampling procedure is proposed by Yuen & Lazarian (2017a), i.e., the subblock averaging method. First, the subblock averaging method takes all gradient orientations within a subblock of interest and then plots the corresponding histogram. A Gaussian fitting is then applied to the histogram. The Gaussian distribution's expectation value is the statistically most probable gradient's orientation within that subblock. Incidentally, the subblock averaging also partially suppresses the rms noise in the spectroscopic data.

The selection of subblock size, which controls the number of sampling points, is crucial. Hu et al. (2020c) later improved the subblock averaging method to be adaptive. The subblock centers were selected continuously, located at the position (i, j), i, j = 1, 2, 3, etc. All gradients pixels within the rectangular boundaries [i − d/2, i+d/2] and [j − d/2, j+d/2] are selected to do averaging, where d is the subblock size. We vary the subblock size and check its corresponding fitting errors within the 95% confidence level. When the fitting error reaches its minimum value, the corresponding subblock size is the optimal selection. We refer to this procedure as the adaptive subblock (ASB) averaging method (Hu et al. 2020c).

By repeating the gradient's calculation and ASB for each eigenchannel, we obtain a total of n_v eigengradient maps ${\psi }_{{gs}}^{i}(x,y)$ with i = 1, 2, ..., n_v . In analogy to the Stokes parameters of polarization, the pseudo Q_g and U_g of gradient-inferred magnetic fields are defined as

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{g}(x,y) & = & \displaystyle \sum _{i=1}^{{nv}}{I}_{i}(x,y)\cos (2{\psi }_{{gs}}^{i}(x,y))\\ {U}_{g}(x,y) & = & \displaystyle \sum _{i=1}^{{nv}}{I}_{i}(x,y)\sin (2{\psi }_{{gs}}^{i}(x,y))\\ {\psi }_{g} & = & \displaystyle \frac{1}{2}{\tan }^{-1}\left(\displaystyle \frac{{U}_{g}}{{Q}_{g}}\right).\end{array}\end{eqnarray} \tag{ 10 }$$

The pseudo polarization angle ψ_g is then defined correspondingly. Similar to the Planck polarization, ψ_B = ψ_g + π/2 gives the POS magnetic field orientation.

The relative alignment between magnetic fields orientation inferred from Planck polarization ϕ_B and velocity gradient ψ_B is quantified by the alignment measure (AM; González-Casanova & Lazarian 2017):

$$\begin{eqnarray}&&\mathrm{AM}=2\left(\langle {\cos }^{2}{\theta }_{r}\rangle -\displaystyle \frac{1}{2}\right),\end{eqnarray} \tag{ 11 }$$

where θ_r = ∣ϕ_B − ψ_B∣ and 〈〉 denotes the average within a region of interest. The value of AM spans from −1 to 1. AM = 1 mean ϕ_B and ψ_B are parallel, while AM = −1 indicates ϕ_B and ψ_B are perpendicular. The standard deviation divided by the sample size's square root gives the uncertainty σ_AM.

The double-peak histogram achieves the detection of the gravitationally collapsing region. Unlike the gradients ψ_B in diffuse regions, the gravitational collapse flips the gradients by 90°, ψ_B + π/2. The histogram of the gradients' angle can be used to extract this change. In the diffuse region, the histogram exhibits a single Gaussian peak located at ψ_B. The single peak of the histogram becomes ψ_B + π/2 in gravitationally collapsing regions. However, in transitional regions, which cover, for instance, half diffuse material and half gravitationally collapsing, the histogram is therefore expected to show two peak values of ψ_B and ψ_B + π/2. We call this type of histogram the double-peak histogram (DPH).

The DPH works as follows. Every single pixel of ψ_B is defined as the center of a subblock ⁸ to draw the histogram. An envelope, which is a smooth curve outlining the extremes, is adopted for the histogram to suppress noise and the effect from insufficient bins. Any term whose histogram weight is less than the mean weight value of the envelope is masked. After masking, we work out the peak value of each consecutive envelope. Once more than one peak value appear and the maximum difference is within the range 90° ± σ_g , where σ_g is the total standard deviation, the center of this second subblock is labeled as the boundary of a gravitationally collapsing region. The details of the algorithm are presented in Hu et al. (2020c).

Figure 1 presents the H₂ column density map for the Serpens G3–G6 clump. The H₂ column density structures are elongated and filamentary. Here we plot the distribution of young stellar objects (YSOs) within this clump. The YSOs are identified by Harvey et al. (2007). The evolutionary stage of a YSO can be classified as (in order of youngest to oldest): Class I, flat, Class II, and Class III (Lada 1987; Andre & Montmerle 1994; Greene et al. 1994). We find Class I and flat YSOs concentrate on the south dense clump core and the northeast filamentary tail. Class II YSOs are mainly located at the north filamentary tail, while Class III only occupies a small fraction. YSOs' high surface density at the dense clump core and the north tail indicates these two regions are actively forming stars.

Zoom In Zoom Out Reset image size

The ¹³CO (1–0) emission line observed with the IRAM 30 m telescope zooms into the south dense clump with a resolution of 235 (≈0.025 pc; see Figure 2). The radial velocity of the emission ranges from about 2.6 to 12.0 km s⁻¹. The spectral line appears to have two apparent peaks: 2.22 K at 7.65 km s⁻¹ and 0.58 K at 4.9 km s⁻¹. After fitting a double Gaussian profile to the integrated spectral line shown in Figure 2, we obtain the 1D velocity dispersion σ_v,1D = 1.1 ± 0.2 km s⁻¹ for the dominating emission feature. Figure 2 also shows the integrated intensity maps for the south dense clump core of Serpens G3–G6. The integration of ¹³CO (1–0) considers pixels where the brightness temperature is larger than 0.9 K, which is about three times the rms noise level. The ¹³CO (1–0) emission lines cover a wider 060 × 015 area. The clump's ¹³CO (1–0) is still filamentary, spanning from north to south.

Zoom In Zoom Out Reset image size

Figure 3 shows representative velocity channels from the ¹³CO (1–0) emission-line data, averaged over 636 m s⁻¹. The intensity structures seen in the velocity channels do not appear filamentary. In particular, two distinct dense structures appear at 7.0 and 7.6 km s⁻¹. These structures suggest that the velocity caustic effect is significant, i.e., velocity fluctuations dominate the thin velocity channels (Lazarian & Pogosyan 2000).

Zoom In Zoom Out Reset image size

The column density PDFs are widely used to study turbulence and self-gravity in ISM. The PDFs in gravitationally collapsing isothermal supersonic turbulence follow a hybrid of log-normal distribution P_LN(s) in the low-intensity range and power-law distribution P_PL(s) in the high-intensity range (Robertson & Kravtsov 2008; Ballesteros-Paredes et al. 2011; Collins et al. 2012; Burkhart 2018; Körtgen et al. 2019):

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{LN}}(s) & = & \displaystyle \frac{D}{\sqrt{2\pi {\sigma }_{s}^{2}}}{e}^{-\tfrac{{\left(s-{s}_{0}\right)}^{2}}{2{\sigma }_{s}^{2}}},s\lt {S}_{t}\\ {P}_{\mathrm{PL}}(s) & = & {{DCe}}^{{ks}},s\gt {S}_{t},\end{array}\end{eqnarray} \tag{ 12 }$$

where $s=\mathrm{ln}(N/{N}_{0})$ is the logarithm of the normalized column density N. S_t is transitional density between the P_LN and P_LN. ${s}_{0}=-\tfrac{1}{2}{\sigma }_{s}^{2}$ is the mean logarithmic density. The standard deviation σ_s for magnetized turbulence is related to the sonic Mach number M_S, driving parameter b, and compressibility β (Molina et al. 2012):

$$\begin{eqnarray}&&{\sigma }_{s}^{2}=\mathrm{log}\left(1+{b}^{2}{{M}_{{\rm{S}}}}^{2}\displaystyle \frac{\beta }{\beta +1}\right).\end{eqnarray} \tag{ 13 }$$

The turbulence-driving parameter b = 1/3 for purely solenoidal driving, while b = 1 for purely compressive driving. For a natural mixture of solenoidal and compressive driving, b is ≈0.4 (Federrath & Banerjee 2015; Federrath et al. 2016). Assuming P_NL(s) + P_PL(s) is normalized, continuous, and differentiable, the coefficient D and C are (Burkhart et al. 2017; Burkhart 2018)

$$\begin{eqnarray}\begin{array}{rcl}D & = & \left(\displaystyle \frac{{{Ce}}^{{S}_{t}k}}{-k}+\displaystyle \frac{1}{2}+\displaystyle \frac{1}{2}\mathrm{erf}\left(\displaystyle \frac{2{S}_{t}+{\sigma }_{s}^{2}}{s\sqrt{2}{\sigma }_{s}}\right)\right)\\ C & = & \displaystyle \frac{{e}^{0.5(k+1)k{\sigma }_{s}^{2}}}{{\sigma }_{s}\sqrt{2\pi }}.\end{array}\end{eqnarray} \tag{ 14 }$$

In Figure 4, we plot the H₂ column density PDFs for the entire Serpens G3–G6 clump. The PDFs are log-normal with dispersion σ_s ≈ 0.50 ± 0.04 until S_t ≈ 0.88, which indicates the gas is gravitationally collapsing when its density is larger than $N={e}^{{S}_{t}}{N}_{0}\approx 1.02\times {10}^{22}$ cm⁻². This column density threshold reveals that the south dense clump is gravitationally collapsing. The characteristic slope k of that power-law is −2.40 ± 0.17. Also, we zoom into the south dense region (see the red box in Figure 4), which corresponding to the measured ¹³CO (1–0). The south dense clump's PDFs are still the combination of P_LN and P_PL. The transition density S_t ≈ 0.66 and the characteristic slope k = −1.95 ± 0.26 identified the same gravitationally collapsing regions as before. The characteristic slope is related to the cloud mean free-fall time, the magnetic fields, and the efficiency of feedback (Federrath & Klessen 2013; Girichidis et al. 2014; Burkhart et al. 2017; Guszejnov et al. 2018). For instance, k = −2.40 ± 0.17 and k = −1.95 ± 0.26 correspond to radial density distributions ρ(r) ∝ r^{−1.83±0.06} and ρ(r) ∝ r^{−2.03±0.14}, respectively. ⁹ The power-law slopes of the radial density distributions suggest that all the high-density gas have collapsed into isothermal cores ρ(r) ∝ r⁻² (Shu 1977) and the star formation efficiency is in the range of 0%–20% (Federrath & Klessen 2013).

Zoom In Zoom Out Reset image size

In Figure 5, we present the results of the VGT analysis using the ¹²CO (2–1) and ¹³CO (2–1) emission-line data sets. We use the emission within the velocity ranges [+1, +13] km s⁻¹ for calculation following the recipe presented in Section 3. Pixels with brightness temperature less than three times the rms noise level (0.4 K) are blanked out. We average the gradients over 20 × 20 pixels and smooth the gradients map ψ_g (see Equation (10)) with a Gaussian filter with kernel width 1, i.e., over 5 pixels.

Zoom In Zoom Out Reset image size

We also compare the magnetic field inferred from the Planck 353 GHz polarized dust signal data (FWHM $\approx 10^{\prime}$). We regrid the Planck polarization further to achieve the same pixel size as emission lines. The polarization vector is also averaged over 20 × 20 pixels to match the gradients map. We find that the resulting gradients of ¹²CO (2–1) have good agreement with the Planck polarization, showing AM = 0.69 ± 0.04. Several anti-alignment vectors appear in the south dense region and image boundary. It is likely that the dust polarization and the ¹²CO (2–1) emission probe different spatial regions. The optically thick tracer ¹²CO samples the outskirt diffuse region of the cloud with volume density n ≈ 10² cm⁻³, while dust polarization likely traces denser regions. The theory of radiative torque (RAT) alignment (see Lazarian & Hoang 2007; Andersson et al. 2015 for a review) predicts that dust grains can remain aligned at high densities, especially in the presence of embedded stars. This is in agreement with numerical simulations (Bethell et al. 2007; Seifried et al. 2019). Therefore, we expect that for n ≥ 10³ cm⁻³, grains are aligned in the regions that we study. In addition, the low resolution ($10^{\prime}$) of the Planck data also contributes to the misalignment. The VGT measurements from ¹³CO (2–1) give less alignment (AM = 0.41 ± 0.07) with the Planck polarization. In particular, the gradients become perpendicular to the magnetic field in the south dense region. As the south dense region is identified to be gravitationally collapsing by the PDFs (see Section 4.2), this change in the gradients' direction suggests the presence of a gravitational collapse (see Section 3). Because the overall magnetic fields horizontally cross the clump, while the gradients are nearly vertical, the resolution effects do not change this conclusion.

We plot the histograms of the relative alignment between rotated gradients and the magnetic fields in Figure 6. The histogram based on the ¹²CO (2–1) data is more concentrated around 0, which means that two vectors are parallel. The distribution based on the ¹³CO (2–1) data, however, spreads all the way to π/2, which means that two vectors are perpendicular. Previous studies reported that the VGT results of the ¹³CO emission are more accurate than the ones derived using the ¹²CO emission when comparing with dust polarization (Hu et al. 2019a; Alina et al. 2020) in nongravitationally collapsing clouds. Here we find that the VGT of ¹²CO gives better alignment. It is likely because the gravitational collapse in the Serpens G3–G6 south clump happens in dense gas n ≥ 10³ cm⁻³. The ¹²CO emission samples more diffuse regions so that its gradients are less affected by self-gravity.

Zoom In Zoom Out Reset image size

The comparison of VGT and Planck polarization directly reveals the region undergoing gravitational collapse. However, insufficient resolution in polarization measurements limits this approach in small-scale studies. For instance, our high-resolution ¹³CO (1–0) can measure the pixelized gradient up to 5'', which is 120 times higher than Planck's resolution, $10^{\prime}$. Nevertheless, the double-peak histogram (see Section 3) extends the VGT to identify the gravitationally collapsing region independent of polarization measurements.

In Figure 7, we present the gradients map and the identified gravitationally collapsing regions using the ¹³CO (1–0) emission line. In dense areas, we find that the gradients flip their direction by 90° compared with surrounding low-intensity regions. The DPH finds three separate gravitationally collapsing regions. The largest gravitationally collapsing region covers the one identified by the PDFs (see Figure 4) and covers parts of low-intensity regions. It is likely the VGT is sensitive to gravity-induced inflows, which also span to low-intensity regions, while the PDFs are detecting the already-formed cores. Also, the DPH covers partially diffuse regions so that the actual gravitationally collapsing area is slightly overestimated. Nevertheless, both the VGT and the PDFs reveal that the central clump is gravitationally collapsing.

Zoom In Zoom Out Reset image size

Gravitational collapse usually shows a signature of infalling motions. Explicitly, an infall signature in a spectral line presents itself in the form of a red–blue asymmetry, generally with a diminished redshifted component (Walker et al. 1994; Myers et al. 1996). To search for the signature, we integrated the ¹³CO (J = 1–0) emission in a velocity range of +3.0 to +7.0 km s⁻¹ for blueshifted gas, as the velocity of bulk motion is 7.65 km s⁻¹ (see Section 4). Redshifted gas is integrated in the velocity range of +8.3 to +12.0 km s⁻¹. In Figure 8, the blueshifted and redshifted gases are overlaid in the integrated intensity map of ¹³CO (J = 1–0). We can see the clump is dominated by blueshifted gas. In the high-intensity region, blueshifted and redshifted gases are overlapped, which indicates an infalling motion along the LOS.

Zoom In Zoom Out Reset image size

Through both PDFs and VGT, we confirm that the Serpens G3–G6 south clump is gravitational collapsing. Here we determine the dynamics of the clump, focusing on the energy balance. The measured and derived physical parameters are listed in Table 1.

Table 1. Physical Parameters of the Serpens G3–G6 Clump

Physical Parameter	Symbol/Definition	Value (uncertainty)	Reference
Area	A	2.14 (0.20) pc2	Measured
1D velocity dispersion	σv,1D	1.10 (0.20) km s−1	Measured
Polarization dispersion	σp	8.54 (0.72) deg	Measured
H2 column density	N0	9.18 (0.43) × 1021 cm−2	Measured
Effective diameter	$L=2\sqrt{A/\pi }$	1.65 (0.08) pc	Derived
Effective radius	R = L/2	0.83 (0.04) pc	Derived
H2 volume number density	n0 = N0/L	1800.2 (117.62) cm−3	Derived
Mass of an H atom	mH	1.67 × 10−24 g	Ref.1
Mean molecular weight	${\mu }_{{{\rm{H}}}_{2}}$	2.8	Ref.1
Volume mass density	${\rho }_{0}={n}_{0}{\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}$	8.42 (0.55) × 10−21 g cm−3	Derived
Mass	$M={N}_{0}{\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}A$	440.8 (45.79) M⊙	Derived
Magnetic field strength	$B=f\sqrt{4\pi {\rho }_{0}}{\sigma }_{v,1{\rm{D}}}/{\sigma }_{p}$	119.88 (24.35) μG	Derived
Kinetic energy	${E}_{K}=3M{\sigma }_{v,1{\rm{D}}}^{2}/2$	1.59 (0.60) × 1046 erg	Derived
Gravitational energy	EG = − 3GM2/(5R)	−1.21 (0.26) × 1046 erg	Derived
Magnetic energy	EB = B2 R3/6	3.97 (1.71) × 1046 erg	Derived
Viral parameter	αvir = ∣2Ek /EG ∣	2.63 (1.14)	Derived
Free-fall time	${t}_{\mathrm{ff}}=\sqrt{3\pi /(32G{\rho }_{0})}$	0.73 (0.02) Myr	Derived
Sound speed (isothermal)	${c}_{s}=\sqrt{{k}_{{\rm{B}}}T/{\mu }_{p}{m}_{{\rm{H}}}}$	188 m s−1	Ref.1 (μp = 2.33)
Alfvén speed	${v}_{{\rm{A}}}=B/\sqrt{4\pi {\rho }_{0}}$	3.69 (0.75) km s−1	Derived
3D sonic Mach number	${M}_{{\rm{S}}}=\sqrt{3}{\sigma }_{v,1{\rm{D}}}/{c}_{s}$	10.1 (1.80)	Derived
3D Alfvén Mach number	${M}_{{\rm{A}}}=\sqrt{3}{\sigma }_{v,1{\rm{D}}}/{v}_{{\rm{A}}}$	0.52 (0.14)	Derived
Compressibility	$\beta =2{({M}_{{\rm{A}}}/{M}_{{\rm{S}}})}^{2}$	0.005 (0.001)	Derived

Note. All physical parameters are derived for pixels that fall within the 28.68 K km s⁻¹ (mean intensity value) ¹³CO (1–0) intensity contours drawn in Figure 2. All uncertainties consider the error propagation among each physical parameter. The gas temperature T is assumed to be 10 K (Draine & Lazarian 1998). The calculation of gravitational energy and magnetic energy assumes that the cloud is spherical and has a uniform density. Reference: (1) Kauffmann et al. (2008).

Download table as:  ASCIITypeset image

The area A of the south clump is measured with ¹³CO (1–0). We include all pixels where their integrated brightness temperature is above the mean value (see Figure 2). Adopting 415 ± 25 pc as the distance (Dzib et al. 2010), we find A ≈ 2.14 (0.2) pc. The value in brackets indicates the uncertainty, which is the average of the upper bound and lower bound. Assuming a simple spherical geometry, we have the effective LOS distance L ≈ 1.65 pc. The 1D velocity dispersion σ_v,1D ≈ 1.1 (0.2) km s⁻¹ is measured from the ¹³CO (1–0) emission line (see Figure 2). σ_v,1D contains the contribution from both turbulence velocity and shear velocity. The polarization dispersion σ_p ≈ 8.54 (0.72) deg is obtained from the Planck 353 GHz polarization; see Figure 9. The mean H₂ column density N₀ ≈ 9.18(0.43) × 10²¹ cm⁻².

Zoom In Zoom Out Reset image size

From these physical parameters, we derive the total mass M = 440.8 (45.79) M_⊙ and volume mass density ρ₀ = 8.42(0.55) × 10⁻²¹ g cm⁻³. Adopting f = 0.5, the total magnetic field strength B is calculated from the Davis–Chandrasekhar–Fermi method (Davis 1951; Chandrasekhar & Fermi 1953), giving $B=f\sqrt{4\pi {\rho }_{0}}{\sigma }_{v,1{\rm{D}}}/{\sigma }_{p}=119.88(24.35)$ μG. Assuming a spherically homogeneous cloud, the corresponding total kinetic energy E_K = 1.59(0.60) × 10⁴⁶ erg, total gravitational energy E_G = − 1.21(0.26) × 10⁴⁶ erg, and total magnetic field energy E_B = 3.96(1.71) × 10⁴⁶ erg. The ratio of kinetic energy and gravitational energy ∣E_K /E_G ∣ ≈ 1.31. In particular, the magnetic field energy to kinetic energy and gravitational energy ratios are ∣E_B/E_K ∣ ≈ 2.50 and ∣E_B/E_G ∣ ≈ 3.28, respectively. The role of the magnetic field in the Serpens G3–G6 south clump is much more significant than turbulence and self-gravity.

Note that the kinetic energy and magnetic field energy may be overestimated, as the dispersion σ_v,1D considers both turbulence velocity and shear velocity. The ratio E_B/E_K after simplification is equivalent to ${E}_{{\rm{B}}}/{E}_{K}=(2{f}^{2})/(9{\sigma }_{p}^{2})$, which is independent of σ_v,1D. Similarly, we have $| {E}_{{\rm{B}}}/{E}_{G}| =(10{f}^{2}{\sigma }_{v,1{\rm{D}}}^{2})/(9\pi {\sigma }_{p}^{2}{{GN}}_{0}{\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}L)$. The overestimated σ_v,1D leads to a stronger magnetic field and also an overestimated ∣E_B/E_G ∣. Nevertheless, here we have the ratio ¹⁰ ∣E_B/E_G ∣ ≈ 3.28; it is unlikely that the velocity dispersion is overestimated by a factor of 2 at least so that ∣E_B/E_G ∣ ≈ 1. Therefore, we expect the estimated magnetic field to be indeed stronger than turbulence and self-gravity. The significant magnetic field energy suggests that the gravitational collapse can happen in a strong magnetic field environment, as numerically suggested by Hu et al. (2020c). ¹¹

In addition, assuming the gas temperature T ∼ 10 K (Draine & Lazarian 1998), we have the isothermal sound speed c_s = 188 m s⁻¹ and Alfvén speed v_A = 3.69(0.75) km s⁻¹. The corresponding sonic Mach number M_S and Alfvén Mach number M_A are therefore 10.1 (1.80) and 0.52 (0.14), respectively. Also, we have the compressibility β = 0.005 (0.001). Recall that we have the measured PDF dispersion σ_s = 0.55 ± 0.05 (see Figure 4). The turbulence-driving parameter b can be derived from Equation (13):

$$\begin{eqnarray}&&{b}^{2}=\displaystyle \frac{({e}^{{\sigma }_{s}^{2}}-1)}{{M}_{{\rm{S}}}^{2}}\displaystyle \frac{\beta +1}{\beta };\end{eqnarray} \tag{ 15 }$$

using our derived parameters, we get b ≈ 0.84, which suggests the turbulence in the Serpens G3–G6 south clump is mainly driven by compressive force. The study in Federrath & Klessen (2012) shows that the supernova-driven turbulence is more effective in producing compressive motions and is expected to be more prominent in regions of enhanced stellar feedback. Therefore, supernovae in this region may contribute to the compressive driving force here. On the other hand, solenoidal motions are more prominent in the quiescent areas with low star formation activity. Intense star formation (i.e., gravitational collapse) can also contribute, although our analysis shows that the gravity-induced motions do not dominate over large volumes. The latter point requires further studies.

The Planck measurement gives an overall view of the magnetic field. However, its resolution limits our scope to smaller scales. The dispersion of polarization measures only large-scale magnetic fields, so that its value may be underestimated. The smooth field lines on a large scale suggest this underestimation is not significant.

In addition to the DCF method, we also estimated the M_A directly from velocity gradients using the new approach proposed in Lazarian et al. (2018b). It was shown that the properties of velocity gradients over the subblock are a function of M_A. In particular, the dispersion of the velocity gradients' orientation was shown to exhibit the power-law relation with M_A, and a different power-law relation was obtained for the so-called "Top-to-Bottom" ratio of the thin channel velocity gradients' (VChGs) distribution:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{A}}} & \approx & 1.6{\left({T}_{v}/{B}_{v}\right)}^{\tfrac{1}{-0.60\pm 0.13}},{M}_{{\rm{A}}}\leqslant 1\\ {M}_{{\rm{A}}} & \approx & 7.0{\left({T}_{v}/{B}_{v}\right)}^{\tfrac{1}{-0.21\pm 0.02}},{M}_{{\rm{A}}}\gt 1,\end{array}\end{eqnarray} \tag{ 16 }$$

where T_v denotes the maximum value of the fitted histogram of the velocity gradient's orientation, while B_v is the minimum value. The analytical justification of these relations is provided in Lazarian et al. (2020). Here we, however, use the empirically obtained dependencies (Hu et al. 2019b). This new way of obtaining the Alfvén Mach number provides us an M_A distribution map for the Serpens G3–G6 south clump, as shown in Figure 10. We find the median value of M_A is 0.62 for the south clump. We also repeat the analysis for the ¹³CO (2–1) emission line and find the median value M_A ≈ 0.53.

Zoom In Zoom Out Reset image size

We note that the difference between evaluating M_A using the new gradient approach in Lazarian et al. (2018b) and the traditional DCF method is that the new technique gets the value of M_A over an individual subblock. Therefore, we get not a single value of M_A, but a distribution of magnetization of M_A over the cloud image. Naturally, the new way of measuring magnetization is much more informative compared to measuring the dispersion of the projected magnetic field over the cloud image. This difference was demonstrated earlier in Hu et al. (2019c), where for a set of molecular clouds the polarization provided the mean value of magnetization and the distribution of velocity gradients provided the detailed maps of magnetization with the averaged value in good agreement with that obtained using polarization.

In the present study, we also have a similar situation. The average value of M_A estimated by the new VGT approach is close to the value that is obtained by measuring the directions of polarization. The latter is M_A = 0.52 ± 0.1. This correspondence is important as the technique in Lazarian et al. (2018b) and the DCF-type measurements of M_A are very different both in terms of the information employed and how this information is processed. This increases our confidence in our results.

Having M_A in hand, one can use a new technique of measuring magnetic strength in Lazarian et al. (2020). The technique is termed MM2 there, as it uses the values of two Mach numbers, the sonic one M_S and the Alfvén one M_A. Using the relation derived in Lazarian et al. (2020), one can evaluate the POS magnetic field as

$$\begin{eqnarray}&&B={\rm{\Omega }}{c}_{s}\sqrt{4\pi {\rho }_{0}}{M}_{{\rm{S}}}\ {M}_{{\rm{A}}}^{-1},\end{eqnarray} \tag{ 17 }$$

where Ω is a geometrical factor. By adopting Ω = 1 (i.e., the magnetic field perpendicular to the LOS), M_A = 0.62 (i.e., measured by VGT), M_S = 10.1, ρ₀ = 8.42 × 10⁻²¹ g cm⁻³, and c_s = 188 m s⁻¹, we get B ≈ 100 μG, which is close to the one (≈120 μG) derived from the DCF method.

Measuring the magnetic field in molecular clouds is generally difficult. The polarization measurement is one practical approach to trace the magnetic fields. It is rooted in the theory of RAT alignment (Lazarian & Hoang 2007; Andersson et al. 2015), which predicts that the polarized dust thermal emission is perpendicular to the magnetic field and the polarized starlight is parallel to the magnetic fields. Through the measured dispersion of the dust polarization directions and the information on spectral broadening, one can estimate the magnetic field through the DCF method (Falceta-Gonçalves et al. 2008; Cho & Yoo 2016). However, the polarization only gives the integrated magnetic field measurement along the LOS.

To achieve the local magnetic field measurement in molecular clouds, several techniques, for instance, the correlation function analysis (CFA; Lazarian et al. 2002; Esquivel & Lazarian 2011), the structure function analysis (SFA; Hu et al. 2021b; Xu & Hu 2021; Hu & Lazarian 2021), and the more recent VGT (González-Casanova & Lazarian 2017; Yuen & Lazarian 2017a; Hu et al. 2018; Lazarian & Yuen 2018a), have been proposed. These techniques are based on the anisotropy of MHD turbulence, i.e., turbulent eddies are elongating along their local magnetic field direction (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999). In particular, for VGT, the velocity gradients of eddies are perpendicular to their local magnetic fields so that the gradients reveal the direction of the magnetic fields. Compared with dust polarization, the VGT, used in this work, exhibits several advantages in tracing magnetic fields. First, it employs spectroscopic data to measure the molecular cloud's local magnetic fields so that it gets rid of the contamination from the foreground (Lazarian & Yuen 2018a; Hu et al. 2020e; Lu et al. 2020). Using spectroscopic data and VGT, one can achieve a higher resolution of the resulting magnetic field map. Also, VGT can be applied to multiple emission lines to reveal how the magnetic field changes with the variation of density (Hu et al. 2019a; Hu & Lazarian 2021).

Observationally identifying gravitationally collapsing regions is notoriously challenging. The column density PDFs provide one possible solution. In isothermal molecular clouds, the PDFs follow a log-normal distribution for nongravitationally collapsing gas. The self-gravity, however, introduces a power-law tail to the high-density range. This power-law tail statistically reveals the critical density threshold above which the gas becomes gravitational collapsing. In this work, we find the gravitationally collapsing gas corresponds to density ${N}_{{{\rm{H}}}_{2}}\approx 1.02\times {10}^{22}$ cm⁻² in the Serpens G3–G6 clump.

VGT provides an alternative way to identify gravitational collapse. In turbulence-dominated diffuse media, the velocity gradients are perpendicular to the magnetic fields (Lazarian & Yuen 2018a). In the presence of gravitational collapse, the gravity-induced velocity acceleration changes the velocity gradients by 90° being parallel to the magnetic fields (Yuen & Lazarian 2017b; Hu et al. 2020c). This change reveals the gravitationally collapsing regions. For instance, VGT confirms that the G3–G6 south clump is gravitationally collapsing here. In addition, utilizing different emission lines, VGT can tell us the volume density range in which the collapse occurs. As shown in Figure 5, the change of gradients appears only in the ¹³CO map, which means the collapse happens at volume density n ≥ 10³ cm⁻³.

The application of VGT is not limited to nearby molecular clouds. Multiple molecular clouds and molecular filaments can be easily separated from spectroscopic PPV cubes. A complete survey of gravitationally collapsing clouds and filaments is achievable using abundant spectroscopic data sets, for instance, the JCMT (Liu et al. 2019), GAS (Kauffmann et al. 2017), COMPLETE (Ridge et al. 2006), FCRAO (Young et al. 1995), ThrUMMS (Barnes et al. 2015), CHaMP (Yonekura et al. 2005), and MALT90 (Foster et al. 2011) surveys.

In addition to velocity gradients, column density gradients provide an alternative way to study the magnetic fields. For instance, Soler et al. (2013) used the histogram of relative orientation (HRO) to characterize the relative orientation of column density gradients and magnetic fields. The HRO technique relies on the polarimetry measurement to get magnetic field information.

The HRO should not be confused with the intensity gradient technique (IGT) which is the outshoot of the VGT (Hu et al. 2019b). The IGT can identify the direction of the magnetic field in subsonic turbulence ¹² and identify shocks in supersonic turbulence. Different from HRO, IGT is independent of polarimetry and keeps the spatial information of density gradients by employing the subblock averaging method, which was first implemented in VGT (Yuen & Lazarian 2017a).

Here we make a comparison of HRO and IGT using the Herschel and Planck data. A numerical comparison is presented in Hu et al. (2019b). For HRO, the density gradients were directly calculated for each pixel following the recipe proposed in Soler et al. (2013). The gradients were then segmented according to their corresponding column density value. After that, we draw the histogram of the relative angle θ between segmented gradients and the magnetic field inferred from Planck polarization. The histogram shape parameter ζ = A_c − A_e , where A_c is the area under the central region of the HRO curve $\left(\tfrac{3}{8}\pi \lt \theta \lt \tfrac{5}{8}\pi \right)$, and A_e is the area in the extremes of the HRO ($0\lt \theta \lt \tfrac{1}{8}\pi$ and $\tfrac{7}{8}\pi \lt \theta \lt \pi$). As for IGT, the pixelized gradient map was further processed by the adaptive subblock averaging method (see Section 3). A similar segmentation was also implemented in IGT, and we use the AM (Equation (11)) to quantify the performance of IGT. Note here that we rotate the density gradients by 90° for IGT so that both positive AM and ζ indicate that the nonrotated gradients are perpendicular to the magnetic field while negative AM and ζ represent a parallel relative orientation.

Figure 11 presents the results of HRO and IGT for the Serpens south clump. We find both ζ and AM are positive when ${N}_{{{\rm{H}}}_{2}}\lt 2.8\times {10}^{21}$ cm⁻². ζ drops to approximately zero in the range of $2.8\times {10}^{21}\lt {{\rm{N}}}_{{{\rm{H}}}_{2}}\lt 5.1\times {10}^{21}$ cm⁻² and then drops to negative value further. The transition density 5.1 × 10²¹ cm⁻² is smaller than the 1.01 × 10²² cm⁻² given by the PDFs. The behavior of AM is similar to ζ, but the change is more dramatic. It is contributed by the subblock averaging method, which enhances the statistically most important components so that density gradients can be used to trace the magnetic field in the diffuse region. The transition from positive AM and ζ to negative AM and ζ is likely caused by self-gravity, as gravitational collapse also flips the density gradients by 90° (Hu et al. 2020c). However, this change can also be caused by shocks, which make it more difficult to use the density gradients in identifying gravitationally collapsing regions than velocity gradients.

Zoom In Zoom Out Reset image size

By analogy to the HRO for column density gradients, Hu et al. (2019b) introduce it to the study of thin channel velocity gradients, called V-HRO. Here we use the velocity gradients calculated from the ¹³CO (1–0) emission line, Herschel column density data, and Planck polarization data to implement the V-HRO. The velocity gradients are calculated per pixel without the subblock averaging method applied and are also segmented based on their corresponding column density value. We find that at the low-density range (${N}_{{{\rm{H}}}_{2}}\lt 4.8\times {10}^{21}$ cm⁻²), the ζ of V-HRO is positive and is larger than the ζ of HRO, which indicates that velocity gradients (unrotated) have better orthogonal alignment with the magnetic fields. At the high-density range, the ζ of V-HRO drops dramatically to be negative and a flip in the velocity gradients' direction occurs. As this change of direction is insensitive to shocks and can be amplified by the subblock averaging method, we therefore can identify gravitationally collapsing regions through velocity gradients (Yuen & Lazarian 2017b; Hu et al. 2020a, Hu et al. 2020c).

Star formation is the subject of hot debate in modern astrophysics. The present two classes of star formation theory differ in the role played by magnetic fields. Mestel & Spitzer (1956) first proposed the strong magnetic field model that clouds, which are mainly regulated by the magnetic fields, are formed with subcritical masses. The magnetic support can be removed if mass could move across field lines. This removal is mainly conducted by neutral gas through the process of ambipolar diffusion (Mestel 1966; Spitzer 1968). The ambipolar diffusion increases the mass to flux ratio. As a consequence, cloud envelopes are subcritical, and cloud cores are supercritical. It, therefore, leads to the low efficiency of star formation. This explanation assumes flux-freezing conditions of the magnetic field in the absence of turbulence. However, the turbulence is ubiquitous in the ISM (Armstrong et al. 1995; Chepurnov & Lazarian 2010; Xu & Zhang 2017; Zhang et al. 2020), and in the presence of turbulence, flux freezing no longer holds, as the process of fast turbulent magnetic reconnection breaks the magnetic freezing (Lazarian & Vishniac 1999). Based on the theory of turbulent reconnection, Lazarian (2005, 2014) and Lazarian et al. (2012) identified reconnection diffusion (RD) as an essential process of removing magnetic flux. Consequently, the total magnetic field may be larger than its viral value, while the region can still be contracting. The RD theory predictions were confirmed in Santos-Lima et al. (2010, 2012, 2013) for molecular clouds and circumstellar disk settings. ¹³ It increases the mass to flux ratio and allows the magnetic field to equalize inside and outside the cloud, decreasing magnetic support. Unlike the ambipolar diffusion, which depends on the ionization, the reconnection diffusion is only related to the turbulent eddies' scale and the turbulent velocities (Lazarian 2005). The weak magnetic field model (Mac Low & Klessen 2004; Padoan & Nordlund 1999) suggests that clouds are formed at the intersection of turbulent supersonic flows. This model requires supercritical clouds so that gravitational collapse can be triggered without removing magnetic support.

Our analysis reveals that the Serpens G3–G6 south clump is supersonic (M_S = 10.1) and sub-Alfvénic (M_A = 0.52). The ratio λ between the actual and critical mass-to-magnetic flux ratios M/Φ is (Crutcher 2004) $\lambda =7.6\times {10}^{-21}{N}_{{{\rm{H}}}_{2}}/B\approx 0.58$. λ < 1 suggests the cloud is subcritical. Also, the total magnetic field energy in the Serpens G3–G6 south clump significantly surpasses the sum of total kinetic energy and gravitational energy (see Table 1), suggesting the magnetic field is relatively strong in this gravitationally collapsing clump.

The typical strong magnetic field model ignores the turbulence effect. Consequently, it is difficult to accumulate enough gas along the field to overcome the magnetic field (Mestel & Spitzer 1956). However, turbulent compression can contribute to material accumulation. Because the strong magnetic fields provide pressures perpendicular to the field lines, the compression preferentially follows the magnetic fields and accumulates material along the perpendicular direction. Once the collapse is triggered, the inflows also predominantly move along the field lines. Consequently, the gravitationally bound cloud will be thin oblate spheroids, as we see in Figure 2. Our results suggest compressive turbulence in the subcritical Serpens G3–G6 south clump. It supports the turbulent, strong magnetic field model.

Another observation of the Serpens south region's magnetic field suggests that magnetic supercriticality sets in at visual extinctions larger than 21 mag, which also indicates that gravitational collapse can be triggered in a strongly magnetized environment (Pillai et al. 2020).

In this work, we use the DCF method to estimate the total magnetic field strength. The DCF method assumes that (i) ISM turbulence is an isotropic superposition of linear small-amplitude Alfvén waves, (ii) turbulence is homogeneous in the region studied, (iii) the compressibility and density variations of the media are negligible, and (iv) the variations of the magnetic field direction and the velocity fluctuations arise from the same region in space. These assumptions could raise uncertainties in our estimation. Nevertheless, the direct measurement from Zeeman splitting reveals that the LOS magnetic field strength is B ≈ 80 μG for H i volume density ≈ 4 × 10³ cm⁻³ (Crutcher 2012). As the magnetic field is amplified in a gravitationally collapsing region, we expect that the uncertainties in our estimated POS magnetic field strength B ≈ 120 μG from the DCF method and B ≈ 100 μG from the MM2 method are not significant.

The measured velocity dispersion includes the contribution from both turbulent velocity and shear velocity. In our calculation, we did not separate the two components. As discussed in Section 4.5, the ratio between kinetic energy and magnetic energy is independent of the velocity dispersion. The uncertainty here mainly comes from the polarization measurement. The velocity dispersion gives contributions to the ratio of magnetic energy and gravitational energy. However, the overestimated velocity dispersion is unlikely to compensate for the difference so that the ratio is less than 1.

Also, we assume the gas temperature T ≃ 10 K. This assumption introduces uncertainties to the calculation of sound speed and sonic Mach number. The gas temperature can be roughly estimated from the optically thick ¹²CO (2–1) emission line (Pineda et al. 2010; Kong et al. 2015):

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\nu }(T) & = & \displaystyle \frac{h\nu /{k}_{{\rm{B}}}}{\exp (h\nu /({k}_{{\rm{B}}}T))-1}\\ T & = & \displaystyle \frac{h{\nu }_{12}}{{k}_{{\rm{B}}}}{\left[\mathrm{log}(1+\displaystyle \frac{h{\nu }_{12}/{k}_{{\rm{B}}}}{{T}_{b}^{12}+{J}_{\nu }({T}_{\mathrm{bg}})})\right]}^{-1}\\ & = & 11.06{\left[\mathrm{log}(1+\displaystyle \frac{11.06}{{T}_{b}^{12}+0.194})\right]}^{-1},\end{array}\end{eqnarray} \tag{ 18 }$$

where h is the Planck constant, ν₁₂ is the emission frequency of ¹²CO (2–1), k_B is the Boltzmann constant, T_b ¹² is the peak intensity of ¹²CO (2–) in units of kelvin, J_ν (T) is the effective radiation temperature, and T_bg = 2.725 K is the temperature of the cosmic microwave background radiation. By using a range of observed T_b values, we estimate that T ranges from 3.7 K to 14.3 K, with a median value of 7.1 K. A species at transition level J = 1–0 is brighter than the one at J = 2–1, so it would correspond to a higher temperature. We, therefore, expect that the assumption T = 10 K does not change significantly our conclusions.