A Holistic Perspective on the Dynamics of G035.39-00.33: The Interplay between Gas and Magnetic Fields

The POS magnetic field is traced by polarized 850 μm continuum data obtained with the SCUBA-2/POL-2 instrument at the JCMT. The SCUBA-2/POL-2 observations of G035.39 (project code: M17BP050; PI: Tie Liu) were conducted from 2017 June to 2017 November using a version of the SCUBA-2 DAISY mapping mode (Holland et al. 2013) optimized for POL-2 observations (POL-2 DAISY mapping mode; Friberg et al. 2016).³⁷ In total, 70 scans were conducted. The beam size of the JCMT at 850 μm is 141. The POL-2 DAISY scan pattern uses a scan speed of 8'' s⁻¹ (compared to 155'' s⁻¹ for a SCUBA-2 DAISY scan pattern) and a fully sampled circular region with a diameter of 12', with a waveplate rotation speed of 2 Hz (Ward-Thompson et al. 2017). The full description of the SCUBA-2/POL-2 instrument and the POL-2 observational mode can be found in Friberg et al. (2016) and Ward-Thompson et al. (2017). Since only the central 3' diameter region has an approximately uniform coverage in the POL-2 DAISY observations, we obtained two adjacent maps to cover G035.39. The central pointings of the two maps are R.A.(J2000) = 18:57:07, decl.(J2000) = +02:11:30 and R.A.(J2000) = 18:57:10, decl.(J2000) = +02:08:00.

Data reduction is performed using a python script called pol2map written within the STARLINK/SMURF package (Chapin et al. 2013), which is specific for submillimeter data reduction (much of it specific to the JCMT). The default pixel size in SCUBA-2/POL-2 observations is 4'', but the final data are gridded to 8'' pixels in pol2map to improve sensitivity. The Stokes Q, U, and I data are all reduced with a filtering-out scale of 200''. The output polarization percentage values are debiased using the mean of their Q and U variances to remove statistical biasing in regions of low signal-to-noise (Kwon et al. 2018; Soam et al. 2018). The details of data reduction with pol2map can be found in Kwon et al. (2018) and Soam et al. (2018). The final co-added maps have an rms noise of ∼1.5 mJy beam⁻¹. The polarization angle θ is measured as θ = 0.5 arctan(U/Q). The angle increases from north toward east, following the IAU convention. Throughout this paper, the polarization orientations obtained are rotated by 90° to show the magnetic field orientation projected on the POS.

We use SCUBA-2 Stokes I 450 and 850 μm continuum data obtained from the legacy survey program "SCOPE" (Liu et al. 2018) and Herschel archival data from the Hi-GAL project (Molinari et al. 2010) to construct the pixel-by-pixel SEDs of the G035.39 field.

The SCUBA-2 observations were conducted on 2016 April 13 under better weather conditions than SCUBA-2/POL-2 observations in 2017. Therefore, the 450 μm data were also obtained. The beam sizes at 450 and 850 μm are 79 and 141, respectively. The pixel sizes are 2'' and 4'' at 450 and 850 μm, respectively. The rms levels at 450 and 850 μm are ∼60 and ∼10 mJy beam⁻¹, respectively.

We use the level 2.5 Herschel/SPIRE (250–500 μm) maps available in the Herschel Science Archive,³⁸ using extended source calibration. The resolutions of the original maps at 250, 350, and 500 μm are approximately 183, 249, and 363, respectively.

Large-scale C¹⁸O (1–0) and ¹³CO (1–0) are used to study the kinematics of the G035.39's natal molecular cloud. The C¹⁸O (1–0) and ¹³CO (1–0) mapping data are obtained from the legacy survey program "TRAO Observations of PGCCs (TOP)" (Liu et al. 2018). The observations were conducted on 2017 March 17. The map size is 30' × 30'. The center of those maps is R.A.(J2000) = 18:57:10, decl.(J2000) = +02:10:00. The FWHM beam size (θ_B) is 47''. The main-beam efficiency (η_B) is 51%. The system temperature during observations is 243 K. The OTF data were smoothed to 0.33 km s⁻¹ and the baseline removed with Gildas/CLASS. The rms level is 0.15 K in antenna temperature (${T}_{A}^{* }$) at a spectral resolution of 0.33 km s⁻¹.

Single-pointing observational data of the HCO⁺ (1–0), H¹³CO⁺ (1–0), and H₂CO (2_1,2–1_1,1) lines are used to investigate the dynamical status of a starless clump in G035.39. The data taken with the Korean VLBI Network (KVN) 21 m telescope (Kim et al. 2011) in Tamna station were obtained on 2017 November 26 in its single-dish mode. The rest frequencies of HCO⁺ (1–0), H¹³CO⁺ (1–0), and H₂CO (2_1,2–1_1,1) lines are 89.18852, 86.754288, and 140.83952 GHz, respectively. The pointing position is R.A.(J2000) = 18:57:11.38, decl.(J2000) = +02:07:27.9. The main-beam sizes at 86 and 140 GHz are 32'' and 23'', respectively. The main-beam efficiencies at 86 and 140 GHz are 44% and 36%, respectively. The data are reduced with Gildas/CLASS. The spectral resolution for both the HCO⁺ (1–0) and H¹³CO⁺ (1–0) lines is ∼0.11 km s⁻¹. The spectral resolution for H₂CO (2_1,2–1_1,1) is 0.07 km s⁻¹. The on-source times for the HCO⁺ (1–0), H¹³CO⁺ (1–0), and H₂CO (2_1,2–1_1,1) observations are 10, 15, and 25 minutes, respectively. The system temperatures during observations are 184, 188, and 171 K, respectively. The achieved rms levels in antenna temperatures are ∼0.05, ∼0.04, and ∼0.03 K, respectively.

We also use the NH₃ (1,1) line data from Sokolov et al. (2017). The GBT beam at NH₃ (1,1) line frequency is 32''. The details of the NH₃ (1,1) observations can be found in Sokolov et al. (2017).

Figure 1 shows the 850 μm Stokes I image. Besides the main filament (outlined with a red thick contour, hereafter denoted as F_M) located at the center of the image, which was identified in previous work (Kainulainen & Tan 2013), the deep SCUBA-2/POL-2 observations reveal several fainter adjacent elongated structures (F_W, F_SW, F_E, and F_NE; see Figure 2) connected to F_M. The skeletons of these elongated structures are identified by using the FILFINDER algorithm (Koch & Rosolowsky 2015) in 850 μm Stokes I emission above 3σ (1σ ∼ 1.5 mJy beam⁻¹). The skeletons are more easily identified in the high-contrast 850 μm Stokes I image than Herschel images because the extended diffuse emission is filtered out in SCUBA-2/POL-2 data. With larger filtering-out scale, more extended emission can be recovered in 850 μm Stokes I data (Liu et al. 2018). Contamination from extended emission in 850 μm continuum will reduce the contrast between the skeletons and the background emission. The FILFINDER algorithm adopting the techniques of mathematical morphology not only can identify the bright filaments but also can reliably extract a population of the faint filaments (Koch & Rosolowsky 2015). The gray thick curves in Figure 2 show the skeletons of the elongated structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

F_M has a length of ∼6.8 pc, measured from its skeleton. The longest elongated structure (outlined with a yellow thick contour in Figure 1; denoted as F_W) having a similar length (∼6.7 pc) to F_M is connected to the northern end of F_M. The mean intensities of F_M and F_W at 850 μm within the 10 mJy beam⁻¹ contours of the Stokes I image are ∼100 mJy beam⁻¹ and ∼24 mJy beam⁻¹, respectively, suggesting that F_W is about four times fainter than F_M.

The POS magnetic field orientations are shown in Figure 2. The field orientations are nearly perpendicular to the major axis of F_M at the middle ridge but tend to be parallel to its major axis at the lower-density tails. The field orientations of the elongated structures (F_SW, F_E, and F_NE) in their denser regions close to the junctions with F_M also tend to be perpendicular to their skeletons. In contrast, the field orientations associated with F_W are more parallel to its major axis. In this paper, we will mainly focus on F_M. More detailed analysis and modeling of magnetic field geometry in the whole G035.39 field will be presented in a forthcoming paper (Juvela et al. 2018b).

Panel (a) in Figure 3 shows the magnetic field orientations associated with only F_M. The magnetic field orientations are more disordered at the two ends and near the edges of the filament. In contrast, the magnetic field orientations become more ordered along the central spine of the filament. We average the orientations with a 16'' pixel boxcar filter as Pattle et al. (2017) did and present the averaged orientations overlaid on a centroid velocity image of NH₃ (1,1) from Sokolov et al. (2017) in panel (b) of Figure 3. We divide F_M into three regions ("N," "M," "S"), which show obvious differences in velocities and magnetic field geometries. "N" and "S" show redshifted and blueshifted line-of-sight velocities with respect to "M." In "N" and "S," the magnetic field orientations are more parallel to the filament skeletons. In contrast, the magnetic field orientations are more perpendicular to the filament skeletons in "M."

Zoom In Zoom Out Reset image size

To investigate how the ordered magnetic field orientations change along the filament, we average the magnetic field orientations with a 32'' pixel boxcar filter and calculate the angles (δθ) between the mean magnetic field orientations and the local orientations of the nearest skeletons. The boxcar filter has a size of 32'' (or 0.45 pc), similar to the radius of the filament. The local orientations of the skeletons were measured from their gradients, similar to what Koch et al. (2012) did. The δθ as a function of offset distance from the northern end of F_M is presented in panel (c) of Figure 3. The magnetic field orientations in the end regions ("N" and "S") are more parallel to the major axis of the filament with δθ ≤ 40°. In contrast, the magnetic field orientations in the central part ("M") are more perpendicular to the major axis of the filament with δθ ≥ 60°. From each end, δθ increases toward the middle part of F_M.

Figure 4 presents histograms of the magnetic field orientations in the whole "M" region (green) and in the subregions (red and blue). The subregions are associated with dense clumps identified in Section 3.3. The blue histogram shows the position angles in the clump "c3" region. The red histogram shows the position angles in the region covering the clumps "c5," "c6," and "c7." Only the orientations with Stokes I intensity larger than 100 mJy beam⁻¹ are included. From a Gaussian fit, we obtain a mean magnetic field orientation of ∼86° and an orientation dispersion (σ_θ) of ∼17° in the whole "M" region. The mean magnetic field orientations ($\overline{\theta }$) in different subregions of the "M" region are nearly the same, mostly perpendicular to the major axis of the filament.

Zoom In Zoom Out Reset image size

Panel (d) of Figure 3 shows how the magnetic field orientation variation (with a default pixel size of 8'') changes with Stokes I intensity. The magnetic field orientations (θ) are subtracted by a mean value ($\overline{\theta }=86^\circ$) in the densest part of "M." Although the uncertainties of orientations caused by noise have been considered in this plot, the orientation variations in low-density regions still show large dispersion, which can be improved by future higher-sensitivity observations. The orientation variation $| \theta -\overline{\theta }|$, however, is nearly constant toward higher Stokes I intensity (>200 mJy beam⁻¹), suggesting that the magnetic field orientations in the central part of the F_M are more perpendicular to the major axis. In contrast, the orientation variations $| \theta -\overline{\theta }|$ become larger toward less dense regions (with Stokes I intensity smaller than 100 mJy beam⁻¹), suggesting that magnetic field orientations in less central regions tend to be more parallel to the major axis.

Column density (N_H2) and dust temperature (T_d) maps of G035.39 are constructed from fitting the Herschel/SPIRE (250, 350, and 500 μm) and JCMT/SCUBA-2 (450 and 850 μm) data with a modified blackbody (MBB) function, assuming a dust emissivity spectral index (β) of 1.8 and the dust mass absorption coefficient κ = 0.1(ν/1 THz)^β cm² g⁻¹ (Juvela et al. 2018a). Only SCUBA-2 data at 450 μm when the signal is above 150 MJy sr⁻¹ and at 850 μm when the signal is above 30 MJy sr⁻¹ in the F_M are used. In such compact, bright regions, SCUBA-2 data are not much affected by spatial filtering. Outside such regions, the SCUBA-2 data are corrected with an offset obtained after comparing an image convolved to the 40'' resolution with a SPIRE-based prediction image at the same resolution.

The observations are fit with a model that consists of surface brightness maps at reference wavelengths and a color temperature map, all with a pixel size of 6''. The model is fit to the observations as a global optimization problem. This involves the convolution of the MBB predictions of the model to the resolution of each of the observed surface brightness maps. After optimization, the final model is convolved to a resolution of 15'', which is close to the resolution of the SCUBA-2 data. More details of the procedure are presented in a forthcoming paper (Juvela et al. 2018b). However, the results are found to be very close to those that would be obtained with the method described in the Appendix of Palmeirim et al. (2013), which similarly tries to maximize the resolution of the resulting column density maps. In our case, the dust temperature is constrained mainly by the SPIRE data, with the shortest wavelength (250 μm) being close to the peak of the spectrum of cold dust emission. However, Juvela et al. (2012) estimated that even with 7% surface brightness errors, the SPIRE data give a better than 1 K accuracy for the temperatures (for T ∼ 15 K), which corresponds only to ∼20% uncertainty in the column density. In the fits, we used the 4% and 10% error estimates for SPIRE and SCUBA-2 data, respectively.

The N_H2 and T_d maps of G035.39 are presented in Figure 5. Although the background emission is subtracted, the above method is still subject to the usual caveats with regard to the line-of-sight temperature variations (Juvela et al. 2018b). The mean dust temperature of 15 K that we derived, however, is very consistent with the mean dust temperature of 14 K (from a small median filter method) and mean kinetic temperature of 13 K derived by Sokolov et al. (2017).

Zoom In Zoom Out Reset image size

The main filament F_M has much higher column density and is colder than its surroundings. The N_H2 and T_d maps of F_M are also presented in the right panel of Figure 5. Here we focus on the highest column density part of the main filament where ${N}_{{{\rm{H}}}_{2}}$ > 7 × 10²¹ cm⁻², a column density threshold for core formation suggested in nearby clouds (André et al. 2014). The mean column density of the subregion is ∼1.8 × 10²² cm⁻². The length (L) and the projected area (A) of the ridge are ∼6.8 and ∼6.9 pc², respectively. Therefore, the total mass (M) of the subregion is calculated as

$$\begin{eqnarray}&&M={N}_{{{\rm{H}}}_{2}}\times A{\mu }_{g}{m}_{{\rm{H}}},\end{eqnarray} \tag{ 1 }$$

where μ_g = 2.8 is the molecular weight per hydrogen molecule and m_H is the mass of a hydrogen atom. The derived mass is ∼2800 M_⊙. Therefore, the line mass (mass per unit length) of the filament is M/L ∼ 410 M_⊙ pc⁻¹, which is within the range (223–635 M_⊙ pc⁻¹) given by Sokolov et al. (2017). Assuming that the filament has cylindrical geometry, the mean radius (r) of the circular end of the cylinder is

$$\begin{eqnarray}&&r=\displaystyle \frac{A}{2L}\approx 0.5\,\mathrm{pc},\end{eqnarray} \tag{ 2 }$$

where A is the projected area of the filament. Therefore, the volume (V) of the cylindrical filament is

$$\begin{eqnarray}&&V=\pi {r}^{2}\times L\approx 5.4\,{\mathrm{pc}}^{3},\end{eqnarray} \tag{ 3 }$$

and the mean volume density (N_H2) is

$$\begin{eqnarray}&&{n}_{{{\rm{H}}}_{2}}=\displaystyle \frac{M}{V{\mu }_{g}{m}_{{\rm{H}}}}\approx 7.3\times {10}^{3}\,{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 4 }$$

In contrast, the average column density and dust temperature of the fainter, western elongated structures F_W are 4.0 × 10²¹ cm⁻² and ∼21 K, respectively, which are calculated within the 10 mJy beam⁻¹ contours of the Stokes I intensity at 850 μm. It is noted, however, that the column density and dust temperature of F_W are less constrained than those of F_M because only Herschel data are used in SEDs fit toward F_W. The mass of F_W is ∼650 M_⊙. Considering the length (∼6.7 pc) of F_W, its line mass and volume density are ∼100 M_⊙ pc⁻¹ and 1.7 × 10³ cm⁻³, respectively. The total mass, line mass, column density, and volume density of F_W are all about four times smaller than those of F_M.

The main filament F_M shows a chain of clumps with nearly even spacing. As shown in the right panel of Figure 5, we extracted nine dense clumps from the SCUBA-2 850 μm image using the FELLWALKER (Berry 2015) source-extraction algorithm, which is a part of the Starlink CUPID package (Berry 2007). The core of the FELLWALKER algorithm is a gradient-tracing scheme that follows many different paths of steepest ascent in order to reach a significant summit, each of which is associated with a clump (Berry 2007). FELLWALKER is less dependent on specific parameter settings than other source-extraction algorithms (e.g., CLUMPFIND; Berry 2007).

The source-extraction process with FELLWALKER is the same as that used by the JCMT Plane Survey, and the details can be found in Moore et al. (2015) and Eden et al. (2017). A mask constructed above a threshold of 3σ in the signal-to-noise ratio (S/N) map is applied to the intensity map as input for the task CUPID:EXTRACTCLUMPS, which extracts the peak and integrated flux density values of the clumps. A further threshold for CUPID:FINDCLUMPS is the minimum number of contiguous pixels, which is set at 12, corresponding to the number of pixels expected to be found in an unresolved source with a peak S/N of 5σ, given a 14'' beam and 4'' pixels.

The coordinates, radii (R_eff), mean dust temperature (T_d), mean column density (${N}_{{{\rm{H}}}_{2}}$), and mean nonthermal velocity dispersion (σ_NT) of the clumps derived from NH₃ (1,1) line emission are presented in Table 1. The mean separation between clumps is ∼0.9 pc. The effective radii of clumps are ${R}_{\mathrm{eff}}=\sqrt{{ab}}$, where a and b are the sizes of the semimajor axis and semiminor axis of the clump from a 2D Gaussian fit, respectively. The beam-deconvolved effective radii (R_eff) of clumps range from 0.12 to 0.31 pc, with an average value of 0.23 pc. The dust temperatures of the clumps range from 13.2 to 17.0 K, with an average value of 14.6 K. Clump masses (M_clump) are derived as

$$\begin{eqnarray}&&{M}_{\mathrm{clump}}=\pi {R}_{\mathrm{eff}}^{2}{N}_{{{\rm{H}}}_{2}}{\mu }_{g}{m}_{{\rm{H}}}.\end{eqnarray} \tag{ 5 }$$

Assuming that the clumps have spherical geometry, their H₂ volume density (n_H2) is derived as

$$\begin{eqnarray}&&{n}_{{{\rm{H}}}_{2}}=\displaystyle \frac{M}{\tfrac{4}{3}\pi {R}_{\mathrm{eff}}^{3}{\mu }_{g}{m}_{{\rm{H}}}}.\end{eqnarray} \tag{ 6 }$$

The volume densities and masses of clumps are also presented in Table 1. The masses of clumps range from 16 to 219 M_⊙, with a mean value of ∼107 M_⊙. The volume densities of clumps range from 8 × 10³ cm⁻³ to 6.4 × 10⁴ cm⁻³, with a mean value of 3 × 10⁴ cm⁻³.

Table 1.  Parameters of Dense Clumps in the Main Filament

Clump	R.A.	Decl.	Reffa	Tdb	${N}_{{{\rm{H}}}_{2}}$ c	σNTd	${n}_{{{\rm{H}}}_{2}}$	Bclump	σA	${{ \mathcal M }}_{{ \mathcal A }}$	Mclump	Mvir	${M}_{\mathrm{vir}}^{B}$
No.	(J2000)	(J2000)	(pc)	(K)	(1022 cm−2)	(km s−1)	(104 cm−3)	(μG)	(km s−1)		(M⊙)	(M⊙)	(M⊙)
c1	18:57:03.84	+02:13:01.2	0.25	17.0	0.8	0.31	0.8	56	0.8	0.6	35	45	78
c2	18:57:06.48	+02:12:18.0	0.12	14.6	1.6	0.30	3.2	142	1.0	0.5	16	19	44
c3	18:57:07.92	+02:11:06.0	0.29	14.7	3.7	0.47	3.1	138	1.0	0.8	219	92	150
c4	18:57:08.40	+02:09:32.4	0.31	15.0	2.2	0.45	1.7	94	0.9	0.8	149	92	144
c5	18:57:08.88	+02:08:27.6	0.13	14.0	3.4	0.70	6.4	219	1.1	1.1	40	81	114
c6	18:57:06.48	+02:08:31.2	0.23	15.4	2.8	0.47	3.0	133	1.0	0.8	104	73	119
c7	18:57:08.40	+02:07:44.4	0.22	14.0	3.6	0.50	4.0	162	1.1	0.8	123	76	124
c8	18:57:11.52	+02:07:19.2	0.28	13.2	3.6	0.39	3.1	138	1.0	0.7	199	64	120
c9	18:57:11.76	+02:06:03.6	0.25	13.9	1.7	0.34	1.7	91	0.9	0.6	75	45	89

Notes.

^a ${R}_{\mathrm{eff}}=\sqrt{{ab}}$, where a and b are the sizes of the semimajor axis and semiminor axis of the clump from a 2D Gaussian fit, respectively. R_eff has been deconvolved with the beam. ^bClump-averaged dust temperature (T_d). ^cMean column density. ^dClump-averaged velocity dispersion derived from line widths of NH₃ (1, 1) emission (Sokolov et al. 2017).

Download table as:  ASCIITypeset image

Panel (a) of Figure 6 shows the averaged spectrum of the ¹³CO (1–0) line emission toward the main filament F_M (enclosed by the green contour in panels (c) to (g) of Figure 6). Multiple velocity components are seen in ¹³CO (1–0) line emission, suggesting that several clouds are overlapped along the line of sight. Those velocity components are well separated by at least 10 km s⁻¹. The four strongest components can be well fitted with Gaussian functions. Their peak velocities, line widths, and brightness temperature are listed in Table 2. The emission around 90 km s⁻¹ is weak and was not fitted.

Zoom In Zoom Out Reset image size

Panel (b) of Figure 6 shows the averaged spectrum of the ¹³CO (1–0) line emission toward the western elongated structures F_W (enclosed by the magenta contour in panels (c) to (g) of Figure 6). Five velocity components are seen in the ¹³CO (1–0) line emission. We fit four components with Gaussian functions and summarize the fitting parameters in Table 2. The emission (highlighted in red) between 50 and 70 km s⁻¹ contains a narrow component and a broad component. The broad component overlaps in velocity with the 45.4 km s⁻¹ component.

Panels (c)–(g) of Figure 6 present the integrated intensity maps of each velocity component in ¹³CO (1–0) line emission. The structures identified in 850 μm continuum emission are mainly associated with the velocity feature at 45 km s⁻¹ (see panel (e)). Previous molecular line studies have also suggested that the main filament F_M is dominated by emission around 45 km s⁻¹ (Sanhueza et al. 2012; Henshaw et al. 2013, 2014, 2017; Jiménez-Serra et al. 2014; Sokolov et al. 2017). The emission around 13 km s⁻¹ (see panel (c)) is very diffuse and is mainly distributed to the west of F_M. The emission around 27 km s⁻¹ (see panel (d)) is mainly located to the north of F_M and is partially overlapped with F_M. The emission around 93 km s⁻¹ (see panel (g)) is overlapped with F_M and F_W but shows a very large velocity difference with respect to the 45 km s⁻¹ component, suggesting that the 93 km s⁻¹ component is not really physically associated with F_M and F_W. Similar to the 45 km s⁻¹ component, the emission around 56 km s⁻¹ (see panel (f)) also shows filamentary structures. Its emission, however, is mainly distributed in the north and northwest of F_M and F_W. A long filament (marked by the yellow dashed line in panel (f)) in the 56 km s⁻¹ cloud overlaps with the southern part of F_M. Since the velocity difference between the 56 km s⁻¹ component and the 45 km s⁻¹ component is larger than 10 km s⁻¹, F_M and F_W may not be greatly affected by the 56 km s⁻¹ cloud. The 56 km s⁻¹ cloud, however, is adjacent to the west part (marked by the green dashed line in panel (e)) of the 45 km s⁻¹ cloud and shows broad emission therein (see panel (b)), which suggests that the 56 km s⁻¹ cloud may be interacting with the 45 km s⁻¹ cloud. The 45 km s⁻¹ cloud shows two elongated structures as shown by the dashed lines in panel (e). The main part, which contains F_M, has a length of at least ∼20 pc. The western part, which contains F_W, has a length of ∼17 pc.

Panel (a) of Figure 7 presents the averaged spectra of C¹⁸O (1–0) line emission in the F_M and F_W regions. The C¹⁸O (1–0) line emission also shows two velocity components at ∼45 km s⁻¹ (see panel (b)) and ∼56 km s⁻¹ (see panel (c)). The ∼56 km s⁻¹ component, however, is very weak and marginally detected toward F_M and F_W. Panels (b) and (c) present the integrated intensity maps of the two velocity components. F_M is clearly associated with the 45 km s⁻¹ emission but offset from the ∼56 km s⁻¹ emission. As shown in panel (d), the C¹⁸O (1–0) line emission around 45 km s⁻¹ has similar morphology to the 850 μm continuum emission. F_W is not obviously seen in C¹⁸O (1–0) emission maps. The C¹⁸O (1–0) emission signal is marginally detected as seen from the averaged spectrum in panel (a). From Gaussian fitting, the C¹⁸O (1–0) line luminosity of F_W is about five times smaller than that of F_M. F_W and F_M, however, show very similar velocity (∼44.9 ± 0.1 km s⁻¹ for both) and line widths (1.8 ±0.4 km s⁻¹ for F_W and 2.1 ± 0.1 km s⁻¹ for F_M) in C¹⁸O (1–0) emission, suggesting that F_W and F_M are kinematically and spatially connected.

Zoom In Zoom Out Reset image size

F_W has ¹³CO (1–0) and C¹⁸O (1–0) line widths similar to those of F_M. F_W, however, is about four times less dense than F_M, suggesting that turbulence in F_W likely plays a relatively more important role compared to F_M.

Table 2 summarizes the parameters of the averaged spectra of ¹³CO (1–0) and C¹⁸O (1–0) line emission toward F_M and F_W, including the centroid velocities, line widths, and peak brightness temperature from the Gaussian fitting. Their mean dust temperatures, column densities, and H₂ number densities are also listed. The column densities of ¹³CO (1–0) and C¹⁸O (1–0) lines are derived using the non-LTE radiation transfer code, RADEX (Van der Tak et al. 2007). The inputs for RADEX are kinetic temperature (T_k), H₂ number density (${n}_{{{\rm{H}}}_{2}}$), molecular line column density (N_line), and line width (Δv). We assume that the T_k is equal to the dust temperature. By fixing T_k, ${n}_{{{\rm{H}}}_{2}}$, and Δv, we can calculate the brightness temperature for a grid of N_line values and compare the model brightness temperature with observed values to find out the best N_line. The derived molecular line column densities and abundances are listed also in Table 2. The mean column densities of ¹³CO in F_M and F_W are ∼2.2 × 10¹⁶ cm⁻² and ∼8.5 × 10¹⁵ cm⁻², respectively. The mean column densities of C¹⁸O in F_M and F_W are ∼2.6 × 10¹⁵ cm⁻² and ∼4.5 × 10¹⁴ cm⁻², respectively. The ¹³CO abundance in F_M is about two times smaller than that in F_W, possibly due to larger opacity in ¹³CO (1–0) emission in F_M. In contrast, F_M and F_W have very similar C¹⁸O abundances. By comparing the observed C¹⁸O abundance (X_obs) with a canonical C¹⁸O abundance (X_cano ∼ 6.1 × 10⁻⁷), Jiménez-Serra et al. (2014) claimed high CO depletion factors (${f}_{{\rm{D}}}={X}_{\mathrm{cano}}/{X}_{\mathrm{obs}}\,\sim$ 5–12) in all three regions of F_M. We derive a mean CO depletion factor of ∼6 for both F_M and F_W if we adopt the same canonical C¹⁸O abundance. F_M is much denser than F_W, and thus CO may be expected to be more depleted in F_M. The C¹⁸O and ¹³CO abundances in both F_M and F_W, however, only vary by a factor of 1–2 in our observations, suggesting that CO in F_M may not be depleted as severely as suggested in previous studies (e.g., Jiménez-Serra et al. 2014). Observations in Jiménez-Serra et al. (2014) have much better spatial resolution than ours, and thus they arguably trace better the inner parts of the filament. Therefore, high CO depletion may occur in densest regions of the filament. The uncertainties of the adopted canonical abundance, however, may prevent an accurate determination of CO depletion levels because observed CO abundances vary cloud by cloud in the Galaxy (Liu et al. 2013c; Giannetti et al. 2014). Finally, the optical depths and excitation temperatures are listed in Table 2. The averaged C¹⁸O (1–0) and ¹³CO (1–0) lines are optically thin for both F_M and F_W. The excitation temperatures of C¹⁸O (1–0) and ¹³CO (1–0) in F_M are ∼15.0 and 14.1 K, respectively. Being similar to the dust temperature, these temperatures suggest that C¹⁸O and ¹³CO lines in F_M are nearly thermally excited and that local thermodynamical equilibrium conditions hold. The excitation temperatures of the C¹⁸O (1–0) and ¹³CO (1–0) line emissions in F_W, however, are ∼14.3 and 14.5 K, respectively, which are lower than the dust temperature (∼21 K), suggesting that C¹⁸O (1–0) and ¹³CO (1–0) lines in F_W are subthermally excited.

A large-scale, smooth velocity gradient of 0.4–0.8 km s⁻¹ pc⁻¹ in the northern part of the main filament (F_M) has been revealed in the ¹³CO, C¹⁸O, and N₂H⁺ line emissions (Henshaw et al. 2014; Jiménez-Serra et al. 2014). Recently, Sokolov et al. (2017) mapped the whole F_M in NH₃ lines and also found a smooth velocity gradient of ∼0.2 km s⁻¹ pc⁻¹ across the whole filament as shown in panel (b) of Figure 3. Several scenarios have been proposed to explain these gradients, including cloud rotation, gas accretion along the filaments, global gravitational collapse, and unresolved subfilament structures (Jiménez-Serra et al. 2014). A very promising scenario that could explain the presence of the smooth velocity gradient would involve the initial formation of substructures inside the turbulent molecular cloud, which interact with each other and may subsequently converge into each other as the cloud undergoes global gravitational collapse (Jiménez-Serra et al. 2014). From high-sensitivity and high spectral resolution molecular line (N₂H⁺ and C¹⁸O) observations, Henshaw et al. (2013) identify several velocity-coherent filaments inside F_M and argue that F_M formed via the collision of two relatively quiescent filaments moving at a relative velocity of ∼5 km s⁻¹.

Based on our large-scale ¹³CO (1–0) map, we argue here that the velocity gradients and the collision of filaments inside F_M are more likely caused by a large-scale (∼10 pc) cloud–cloud collision. F_M itself is also formed owing to the large-scale cloud–cloud collision.

From the channel maps of ¹³CO (1–0) line emission (see the Appendix), we identify two velocity-coherent clouds (G035.39-main and G035.39-west) whose spatial distributions are distinctly different. Panel (a) of Figure 8 presents the moment 1 map of ¹³CO (1–0) line emission toward the G035.39 clouds. Blueshifted emission (G035.39-main) and redshifted emission (G035.39-west) are clearly separated in the northern part of the map, indicating two clouds in collision. The southern part of G035.39-main is not affected by cloud–cloud collision and thus shows no clear velocity gradient. Hernandez & Tan (2015) argue that the large-scale velocity gradient in G035.39 is caused by cloud rotation. By inspecting the channel maps, however, the blueshifted emission gas and redshifted emission gas more likely belong to two different clouds. In addition, multiple velocity components in F_M have been seen in high spectral resolution observations (Henshaw et al. 2013), which cannot be well explained by the cloud rotation scenario.

Zoom In Zoom Out Reset image size

The integrated intensity maps of ¹³CO (1–0) for G035.39-west (46–47 km s⁻¹) and G035.39-main (41–43 km s⁻¹) are shown in panel (b) of Figure 8. In the colliding area, the cloud gas with redshifted velocities is well separated spatially from the cloud gas with blueshifted velocities. F_M is located in the interface layer, where the internal turbulence and the momentum exchange between the two colliding clouds may mix the gas distribution in both space and velocity and enhance the density therein. G035.39-west is curved as depicted by the yellow dashed line, suggesting that it is greatly compressed as it collides with G035.39-main. The widespread SiO emission in the northern part of F_M discovered by Jiménez-Serra et al. (2010) may be a sign of (large-scale) shocks from the resulting compression. The emission peak of G035.39-west is in the interface, suggesting that the majority of the gas of G035.39-west may have merged with the gas of G035.39-main. G035.39-west seems to have swung during collision, forming a long tail in the west (as seen in the red dashed box). In addition, the northern ("N") part of F_M appears to be more affected by the collision than the middle and southern parts ("M" and "S") because the emission peak of G035.39-west is located close to the north end of F_M. Moreover, the northern part of F_M is found to have more redshifted velocities than the southern part, as seen from the moment 1 map of ¹³CO (1–0) line emission. We suspect that the collision occurs from the northwest of F_M and slows down as it propagates to the south, causing a velocity gradient along F_M. Therefore, the previous findings of velocity gradients (Henshaw et al. 2014; Jiménez-Serra et al. 2014; Sokolov et al. 2017) and multiple velocity components (Henshaw et al. 2013) in F_M can be explained by the mixed gas distribution from the two larger-scale colliding clouds.

The relative velocity between G035.39-west and G035.39-main is ∼5 km s⁻¹, which is similar to that of the two colliding filaments suggested by Henshaw et al. (2013). Considering the projection effect, the collision velocity of the two clouds could be ∼10 km s⁻¹ (for an inclination angle of 45°), consistent with the collision speeds of GMCs in some simulations (Inoue & Fukui 2013; Wu et al. 2015, 2017). A cloud–cloud collision in G035.39 is also supported by [C ii] observations (Bisbas et al. 2018).

A schematic illustration of the cloud–cloud collision is shown in panels (a) and (b) of Figure 9. The smaller G035.39-west cloud collides with the northern part of G035.39-main in a northwest–southeast direction. The cloud–cloud collision enhances the density in the interface, where the massive filament F_M has formed. The dynamical effect of the cloud–cloud collision may perturb F_M and trigger its fragmentation. In contrast, the southern part (blue dashed box in panel (b) of Figure 8) of G035.39-main is not affected by the cloud–cloud collision, and its density is not enhanced, as seen from our ¹³COmap, as well as infrared extinction maps (Kainulainen & Tan 2013). No dense clump (or new stars) has been formed there yet. Therefore, cloud–cloud collision can enhance density and shorten the local freefall timescale for star formation.

Zoom In Zoom Out Reset image size

As discussed in Section 3.1, the magnetic field orientations are roughly perpendicular to the major axis of F_M along the central ridge and at the junctions with other filaments, while magnetic field orientations tend to be oblique in the lower-density surroundings (see panel (d) in Figure 3). To explore the large-scale field geometry further, we convolve the Stokes Q, U, and I maps of POL-2 data with a 2' beam and recalculate the polarization angles from the smoothed maps. The smoothed magnetic field orientations are overlaid onto the similar Stokes I image in Figure 10. The smoothed magnetic field orientations may indicate that the magnetic field is pinched around the middle and southern parts (marked by the blue dashed ellipse in Figure 10) of F_M. The pinched magnetic field can be associated with the accretion flow around and along the filament in a globally collapsing cloud (Klassen et al. 2017; Gomez et al. 2018; P. S. Li et al. 2018, in preparation). Therefore, the magnetic fields in these regions seem to have been dragged by the collapsing gas flow that forms those dense structures (Li et al. 2015b, 2018, in preparation; Klassen et al. 2017; Gomez et al. 2018).

Zoom In Zoom Out Reset image size

As shown in Figure 1 in Gomez et al. (2018), the magnetic field lines can be dragged by the accretion flow. In low-density regions away from a filament, the gas flow direction is perpendicular to the filament, and the dragged magnetic field must be mainly perpendicular to the filament as seen in the surroundings of other filamentary clouds (Chapman et al. 2011; Cox et al. 2016; Santos et al. 2016). In the low-density regions surrounding the filament spine, the magnetic field is affected by accretion flows both onto and along the filament, and thus the magnetic field lines must also develop a component parallel to the filament and appear oblique. Accretion flows along the filament can compress the magnetic field at the filament spine and increase the perpendicular component as seen in simulations. This picture, illustrated in panel (c) of Figure 9, can well explain the magnetic field geometry in the middle part of F_M, where the projected magnetic field lines are mainly perpendicular to the filament along its densest regions and are oblique in less dense regions (see panel (d) in Figure 3).

As shown in panel (b) of Figure 3, the magnetic field orientations in the northern end of F_M are nearly parallel to the filament. This pattern cannot be explained by gas flows along the dense filament because such flows would not increase the parallel component of the magnetic field, but would rather increase the perpendicular component, leading to a "U"-shaped magnetic geometry (Gomez et al. 2018; P. S. Li et al. 2018, in preparation). As marked by the green dashed box in Figure 10, the smoothed magnetic field orientations in a large area close to the northern end of the F_M are well aligned along the northwest–southeast direction. Indeed, such a nearly parallel pattern may be due to compression from the east (see panel (c) in Figure 9). Since the northern part of F_M is more affected by the cloud–cloud collision (see Section 4.1), the northern end of F_M may be elongated and compressed by the two colliding clouds. The magnetic field at this location is well aligned with the elongated filament (see Figure 9(c)).

4.2.1. Pinched Magnetic Field Surrounding Clump "c8"

The magnetic field surrounding clump "c8" is shown in panel (a) of Figure 11. The magnetic field orientations associated with "c8" are likely pinched. The pinched magnetic field may hint at gas inflows toward the center of "c8." The yellow dashed arrows in the panel indicate the possible gas flow directions. Interestingly, the magnetic field orientations are roughly parallel to the suggested gas flow directions, a behavior seen also in simulations (Klassen et al. 2017; Gomez et al. 2018; P. S. Li et al. 2018, in preparation). In addition, the magnetic field orientations close to the center become more perpendicular to the major axis of the local filament, suggesting that magnetic fields therein have been compressed by accretion flows and become more perpendicular as a result (Klassen et al. 2017; Gomez et al. 2018; P. S. Li et al. 2018, in preparation). Indeed, the magnetic fields surrounding "c8" hint at a "U"-shaped geometry caused by accretion flows (Gomez et al. 2018; P. S. Li et al. 2018, in preparation). The coarse resolution of our POL-2 observations, however, cannot resolve the field geometry of the clumps. In addition, there are only a handful of high-S/N orientations, too few to well constrain the field geometry. Future higher-sensitivity and higher-resolution polarization observations are needed to investigate the field geometry in greater detail.

Zoom In Zoom Out Reset image size

In panel (b) of Figure 11, we show the magnetic field surrounding an accretion core from a simulation of an IRDC (Li et al. 2015b, 2018, in preparation). This image is taken from an ideal MHD turbulence simulation driven at a thermal Mach number of 10 and an Alfvén Mach number of 1. The image is selected at half of the freefall time of the simulation. The core (bright yellow region) is accreting gas along the filament. The white dashed arrow shows the accretion direction. The magnetic field has been twisted along the accretion direction. Close to the densest part of the core, the magnetic field is compressed and is roughly perpendicular to the filament. The accretion significantly increases the perpendicular component of the magnetic field at the filament spine. Although "c8" is much more massive and larger than the simulated core, the magnetic field surrounding "c8" shows similar geometry to that associated with the simulated core, suggesting that the magnetic field surrounding "c8" has been similarly pinched by gas inflow along the filament.

Using ALMA, Henshaw et al. (2017) resolved clump "c3" into a network of narrow (∼0.028 ± 0.005 pc in width) subfilaments that contain 28 compact cores. Those cores may be still accreting gas along those subfilaments. Therefore, we suspect that the magnetic fields surrounding those cores are also pinched owing to the accretion, as seen in simulations (see panel (b) of Figure 11). Future polarization observations with ALMA will shed light on the magnetic field geometry surrounding these cores.

Most stability analyses of massive filaments suggest that the magnetic field is important, but a thorough analysis has been elusive given the difficulties in observing magnetic fields (e.g., Jackson et al. 2010; Contreras et al. 2016; Henshaw et al. 2016; Lu et al. 2018). In this section, we investigate the gravitational stability of the main filament F_M by taking into account thermal pressure, turbulence, and magnetic fields.

The critical line mass for the global gravitational stability of an isothermal filament supported by thermal pressure and turbulence is (Ostriker 1964; Pattle et al. 2017)

$$\begin{eqnarray}&&{\left(\displaystyle \frac{M}{L}\right)}_{\mathrm{crit}}=\displaystyle \frac{2{\sigma }_{3{\rm{D}}}^{2}}{G},\end{eqnarray} \tag{ 7 }$$

where G is the gravitational constant. Assuming that the velocity dispersion is isotropic, the three-dimensional velocity dispersion is

$$\begin{eqnarray}&&{\sigma }_{3{\rm{D}}}=\sqrt{3({\sigma }_{\mathrm{NT}}^{2}+{c}_{s}^{2})}.\end{eqnarray} \tag{ 8 }$$

We obtain a mean nonthermal velocity dispersion (σ_NT) of ∼0.4 km s⁻¹ from the line widths of the NH₃ (1,1) line (Sokolov et al. 2017). The 1D thermal velocity dispersion (or sound speed) is

$$\begin{eqnarray}&&{c}_{s}=\sqrt{\displaystyle \frac{{{kT}}_{k}}{\mu {m}_{{\rm{H}}}}}\end{eqnarray} \tag{ 9 }$$

where μ = 2.37 is the mean molecular weight per "free particle" (H₂ and He; the number of metal particles is negligible). c_s is ∼0.23 km s⁻¹ for a mean kinetic temperature (T_k) of 15 K. Here we assume that T_k equals T_d under the local thermodynamic equilibrium (LTE) conditions. Therefore, the σ_3D and ${(M/L)}_{\mathrm{crit}}$ are 0.80 km s⁻¹ and ∼296 M_⊙ pc⁻¹, respectively. The ${(M/L)}_{\mathrm{crit}}$ is smaller than the measured line mass (∼410 M_⊙ pc⁻¹), suggesting that the filament cannot be supported against collapse only by turbulent gas pressure in the absence of magnetic fields.

In contrast, the F_W filament is as turbulent as the F_M but is much less dense. The critical line mass for F_W is comparable to that of F_M, while its line mass is only ∼100 M_⊙ pc⁻¹. Therefore, turbulent motions in F_W can dominate over gravity in F_M and may further stretch F_W.

The criterion for filament stability with support from magnetic fields can be estimated as (Ostriker 1964; Fiege & Pudritz 2000; Pattle et al. 2017)

$$\begin{eqnarray}&&{\left(\displaystyle \frac{M}{L}\right)}_{\mathrm{crit},\mathrm{mag}}={\left(\displaystyle \frac{M}{L}\right)}_{\mathrm{crit}}{\left(1-\displaystyle \frac{{{ \mathcal M }}_{{ \mathcal B }}}{| { \mathcal W }| }\right)}^{-1},\end{eqnarray} \tag{ 10 }$$

where ${{ \mathcal M }}_{{ \mathcal B }}$ is the magnetic energy per unit length and $| { \mathcal W }| =\left(\tfrac{M}{L}\right)G\approx 4.6\times {10}^{26}$ erg cm⁻¹ is the gravitational energy per unit length. ${{ \mathcal M }}_{{ \mathcal B }}$ may be either positive or negative, depending on whether a poloidal or toroidal field, respectively, dominates the overall magnetic energy (Fiege & Pudritz 2000). If a poloidal field dominates, the factor ${\left(1,-,\tfrac{{{ \mathcal M }}_{{ \mathcal B }}}{| { \mathcal W }| }\right)}^{-1}$ is larger than 1 (Fiege & Pudritz 2000). Therefore, a poloidal field helps to support the cloud radially against self-gravity by increasing the critical mass per unit length (Fiege & Pudritz 2000). In contrast, the factor ${\left(1,-,\tfrac{{{ \mathcal M }}_{{ \mathcal B }}}{| { \mathcal W }| }\right)}^{-1}$ is smaller than 1 (Fiege & Pudritz 2000) for a toroidal field. A toroidal field works with gravity in squeezing the cloud by reducing the critical mass per unit length (Fiege & Pudritz 2000).

The total magnetic field strength (B_tot) can be estimated using the Davis–Chandrasekhar–Fermi (DCF) method (Davis & Greenstein 1951; Chandrasekhar & Fermi 1953):

$$\begin{eqnarray}&&{B}_{\mathrm{tot}}=1.3{B}_{\mathrm{pos}}=1.3\,Q^{\prime} \sqrt{4\pi \rho }\displaystyle \frac{{\sigma }_{\mathrm{NT}}}{{\sigma }_{\theta }},\end{eqnarray} \tag{ 11 }$$

where B_pos is the POS magnetic field strength, 1.3 is a factor considering projection effects, Q' is a factor of order unity accounting for variations in field strength on scales smaller than the beam (Crutcher et al. 2004), $\rho ={\mu }_{g}{m}_{{\rm{H}}}{n}_{{{\rm{H}}}_{2}}$ is the gas density, and σ_θ is the dispersion in polarization position angles. Here Q' is taken as 0.5 (Ostriker et al. 2001).

As seen in Figure 3, the magnetic field orientations in the middle part of F_M are well ordered and uniform, with their orientations roughly perpendicular to the major axis of the filament. In addition, the magnetic field orientations tend to be more perpendicular toward the denser regions (see panel (d) of Figure 3). In contrast, the magnetic field orientations in the northern and southern ends are more widely dispersed in direction. We only estimated σ_θ in the middle part of the main filament, where the Stokes I intensity at 850 μm is above 100 mJy beam⁻¹. As shown in panel (a) of Figure 4, from a Gaussian fit to the orientation angles, the measured σ_θ is ∼17°. We correct the angular dispersion σ_θ by subtracting the measurement uncertainty (δ_θ ∼ 9°) with $\sqrt{{\sigma }_{\theta }^{2}-{\delta }_{\theta }^{2}}$. The corrected σ_θ is ∼15°, which is smaller than the maximum value at which the standard DCF method can be safely applied (≤25°; Heitsch et al. 2001). Taking σ_θ as ∼15°, we obtain a POS magnetic field strength (B_pos) of 50 μG and hence a total magnetic field strength (B_tot) of 65 μG.

The estimated magnetic field strength should be treated as an upper limit because (1) σ_θ is estimated from the densest region of F_M with uniform magnetic field orientations. The magnetic field orientations in other parts of F_M are more widely dispersed and should have larger σ_θ. Therefore, σ_θ used in the above calculations with the DCF method is a lower limit. (2) σ_NT is estimated from the mean line width of NH₃ (1, 1) (Sokolov et al. 2017). The GBT beam at NH₃ (1, 1) line frequency is 32'', larger than the 141 beam of SCUBA-2/POL-2 at 850 μm. Therefore, some uncertainties remain in σ_NT for the estimation of magnetic field strength using the DCF method at this scale.

Crutcher et al. (2010) obtained an empirical relation for a maximum field strength (B_max) versus density from Zeeman observations:

$$\begin{eqnarray}&&{B}_{\max }\simeq 0.22{\left(\displaystyle \frac{{n}_{{\rm{H}}}}{{10}^{5}{\mathrm{cm}}^{-3}}\right)}^{0.65}\,\mathrm{mG}\,\,({n}_{{\rm{H}}}\gt 300\,{\mathrm{cm}}^{-3}).\end{eqnarray} \tag{ 12 }$$

For a mean n_H of ∼1.5 × 10⁴ cm⁻³ in F_M, the maximum field strength B_max estimated from the empirical relation is 64 μG, i.e., similar to the value (65 μG) derived from the DCF method. Therefore, 65 μG may represent an upper limit for the mean magnetic field strength of the main filament F_M.

Although lacking an accurate determination of magnetic field strength, we can use the upper limit of 65 μG to evaluate the importance of magnetic field in the gravitational stability of the main filament.

The corresponding Alfvénic velocity for a magnetic field strength of 65 μG is

$$\begin{eqnarray}&&{\sigma }_{A}=\displaystyle \frac{{B}_{\mathrm{tot}}}{\sqrt{4\pi \rho }},\end{eqnarray} \tag{ 13 }$$

where σ_A ≈ 1.0 km s⁻¹ for B_tot = 65 μG and a mean volume density of 7.3 × 10³ cm⁻³. The Alfvén Mach number is

$$\begin{eqnarray}&&{{ \mathcal M }}_{{ \mathcal A }}=\sqrt{3}{\sigma }_{\mathrm{NT}}/{\sigma }_{A}.\end{eqnarray} \tag{ 14 }$$

We derived ${{ \mathcal M }}_{{ \mathcal A }}\approx 0.7$, suggesting that the turbulent motions may be sub-Alfvénic in the main filament. The total magnetic energy (E_B) is (Pattle et al. 2017)

$$\begin{eqnarray}&&{E}_{B}=\displaystyle \frac{{B}_{\mathrm{tot}}^{2}V}{2{\mu }_{0}},\end{eqnarray} \tag{ 15 }$$

where μ₀ is the permeability of free space. Therefore, the magnetic energy per unit length is

$$\begin{eqnarray}&&| {{ \mathcal M }}_{{ \mathcal B }}| =\displaystyle \frac{{E}_{B}}{L}.\end{eqnarray} \tag{ 16 }$$

Considering the volume (V = 5.4 pc³) and length (L = 6.8 pc) of F_M, E_B and $| {{ \mathcal M }}_{{ \mathcal B }}|$ are ∼2.7 × 10⁴⁶ erg and ∼1.3 ×10²⁶ erg cm⁻¹, respectively. Therefore, if a poloidal field component dominates, it will increase the critical line mass by a factor of ${\left(1,-,\tfrac{1.3\times {10}^{26}}{4.6\times {10}^{26}}\right)}^{-1}\sim 1.39$. In this case, the critical line mass, taking into account the additional support from magnetic fields, will become ∼411 M_⊙ pc⁻¹, which is very similar to the measured value (∼410 M_⊙ pc⁻¹). If, however, a toroidal field component dominates, it will decrease the critical line mass by a factor of ${\left(1,-,\tfrac{-1.3\times {10}^{26}}{4.6\times {10}^{26}}\right)}^{-1}\sim 0.78$. If so, the critical line mass will become ∼231 M_⊙ pc⁻¹, much smaller than the measured value (∼410 M_⊙ pc⁻¹).

Judging from panels (a) and (b) in Figure 3, the northern part ("N") of F_M has magnetic field orientations parallel to the major axis, indicating that the magnetic field therein is likely poloidal. Therefore, the northern part may be stable with additional support from magnetic fields.

In contrast, the middle part of F_M is dominated by a magnetic field whose orientation is perpendicular to the major axis, resembling the projection of a toroidal magnetic field wrapping around the filament. If so, the middle part may become more unstable. Alternatively, the field can also just simply go straight through the filament, providing no support against gravitational collapse. Therefore, the middle part is likely unstable and may further fragment or collapse. Indeed, clump "c3" in the middle part already contains plenty of substructures (subfilaments and cores) detected in ALMA observations (Henshaw et al. 2017).

The southern part of F_M is more complicated. In contrast to the middle part, the magnetic field orientations in the southern part are more parallel to the major axis. We note, however, that the magnetic field orientations tend to be more perpendicular to the major axis in the densest region of the southern part. Therefore, the magnetic field in the southern part may contain comparable toroidal and poloidal components. The critical line mass considering magnetic field support will not deviate too much from that without magnetic field support. Hence, the southern part may be unstable and may fragment or collapse.

In this section, we investigate the gravitational stability of dense clumps from virial analysis by taking into account the support from thermal pressure, turbulence, and magnetic fields.

If we only consider support from thermal pressure and turbulence, the virial masses (M_vir) of the clumps, assuming a uniform density profile, are (Bertoldi & McKee 1992; Pillai et al. 2011; Sanhueza et al. 2017)

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}=\displaystyle \frac{5{R}_{\mathrm{eff}}}{G}({\sigma }_{\mathrm{NT}}^{2}+{c}_{s}^{2}).\end{eqnarray} \tag{ 17 }$$

The virial masses calculated are presented in Table 1. Three clumps ("c1," "c2," and "c5") have virial masses that are two to three times larger than their clump masses and hence may be gravitationally unbound, suggesting that "turbulent" gas motions in the clumps provide enough support against self-gravity. The other clumps have virial masses smaller than their clump masses, suggesting that they are bound and unstable without additional support from magnetic fields.

Henshaw et al. (2016) also suggested that the dense cores revealed in high-resolution interferometric observations are susceptible to gravitational collapse without additional support from magnetic fields. Therefore, it is important to evaluate the importance of magnetic fields in the gravitational stability of the dense clumps. Our SCUBA-2/POL-2 observations, however, do not resolve the magnetic fields surrounding those clumps, and thus we do not have estimation of their magnetic field strengths from observations. Instead, we estimate the magnetic field strengths with the empirical relation from Crutcher et al. (2010) and Li et al. (2015b).

Based on their MHD simulation results, Li et al. (2015b) suggested that the average field strength (B_clump) in molecular clumps in the interstellar medium is

$$\begin{eqnarray}&&{B}_{\mathrm{clump}}\simeq 42{\left(\displaystyle \frac{{n}_{{\rm{H}}}}{{10}^{4}{\mathrm{cm}}^{-3}}\right)}^{0.65}\,\mu {\rm{G}}.\end{eqnarray} \tag{ 18 }$$

Using Equation (18), we estimated the total magnetic field strength B_clump for clumps. The B_clump and Alfvénic speed σ_A of clumps are listed in Table 1. The B_clump values range from ∼56 to 219 μG. The mean Alfvén Mach number of clumps is ∼0.75, suggesting that the magnetic field may play a role as important as turbulence in supporting clumps against gravity. To investigate the gravitational stability of those dense clumps, we estimated the virial masses (${M}_{\mathrm{vir}}^{B}$) of the clumps considering thermal, turbulent, and magnetic pressures and assuming a uniform density profile (Bertoldi & McKee 1992; Pillai et al. 2011; Sanhueza et al. 2017):

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}^{B}=\displaystyle \frac{5{R}_{\mathrm{eff}}}{G}\left({\sigma }_{\mathrm{NT}}^{2}+{C}_{s}^{2}+\displaystyle \frac{{\sigma }_{A}^{2}}{6}\right).\end{eqnarray} \tag{ 19 }$$

${M}_{\mathrm{vir}}^{B}$ are also presented in Table 1. Three clumps ("c1," "c2," and "c5") have virial masses two to three times larger than their clump masses and hence may be gravitationally unbound, suggesting that "turbulent" gas motions and magnetic fields in them provide significant support against self-gravity. The most two massive clumps ("c3" and "c8"), however, have clump masses larger than its virial masses, suggesting that they will collapse and fragment. The other clumps have virial masses comparable to their clump masses, suggesting that they are close to virial equilibrium with additional support from magnetic fields.

Clump "c8" is particularly interesting because it is not visible at Herschel/PACS 70 and 160 μm bands and Spitzer/MIPS 24 μm band, indicating that it is very cold and maybe starless. The physical parameters (e.g., mass, density, size) of "c8" are similar to other Galactic massive starless clumps discovered in large surveys (e.g., Guzmán et al. 2015; Traficante et al. 2015; Contreras et al. 2017; Yuan et al. 2017). As noted earlier (Section 4.2.1), the magnetic field surrounding "c8" is pinched, hinting at gas inflow along the filament. The virial parameter (α_vir) of "c8" is α_vir = M_vir/M_clump ≤ 0.6 even if we consider additional support from magnetic fields, suggesting that "c8" is undergoing gravitational collapse.

Figure 11 shows evidence of the collapse of "c8" from the asymmetric "blue-skewed profiles" of optically thick lines (HCO⁺ (1–0) and H₂CO (2_1,2 – 1_1,1) from KVN observations. The systemic velocity of "c8" is 44.9 km s⁻¹, which is determined from Gaussian fitting to the single-peaked H¹³CO⁺ (1–0) line. In contrast, HCO⁺ (1–0) and H₂CO (2_1,2 – 1_1,1) show double-peaked emission, with the blueshifted peak stronger than the redshifted one, and typical "blue-skewed profiles" for infall signature (Zhou et al. 1993). Such a "blue-skewed profile" of optically thick lines is commonly seen in surveys toward massive clumps (Wu & Evans 2003; Fuller et al. 2005; Wu et al. 2007; Jin et al. 2016; Liu et al. 2016a), which can be interpreted as evidence for the global collapse of massive clumps (Liu et al. 2013a; Peretto et al. 2013). We highlight that, to our knowledge, "c8" could be the first discovered massive starless clump candidate exhibiting this characteristic infall profile.

We model the HCO⁺ (1–0) and H₂CO (2_1,2 – 1_1,1) lines using RATRAN following Peretto et al. (2013) and Yuan et al. (2018). For the modeling, a power-law density profile (ρ ∝ r^−1.5) is assumed, and the kinetic temperature is set to be 13 K. We have tried a grid of models by varying molecular abundances, infall velocities, and velocity dispersions. The resulting infall velocity inferred from the best models for HCO⁺ (1–0) is 0.32 ± 0.04 km s⁻¹, while the resulting infall velocity derived from the best models for H₂CO (2_1,2 – 1_1,1) is 0.20 ± 0.10 km s⁻¹. Though the values are arguably the same within the uncertainties, the infall velocity traced by H₂CO (2_1,2 – 1_1,1) is smaller than that traced by HCO⁺ (1–0). Since the effective excitation density (1.5 × 10⁵ cm⁻¹) of H₂CO (2_1,2 – 1_1,1) at 10 K is much larger than that (9.5 × 10² cm⁻¹) of HCO⁺ (1–0) (Shirley 2015), H₂CO (2_1,2 – 1_1,1) should trace denser, inner regions of the clump than HCO⁺ (1–0). Therefore, the gas inflow indicated by H₂CO (2_1,2 – 1_1,1) and HCO⁺ (1–0) seems to be decelerated from the outer part to the inner part. The decelerated inflow may be caused by the enhanced magnetic field strength near the clump center, which will help resist gravity. Considering the uncertainties in the infall velocities, future work is needed.

Assuming a power-law density profile (ρ ∝ r^−1.5), the mass enclosed in r_o is

$$\begin{eqnarray}&&M={\int }_{0}^{{r}_{o}}4\pi {r}^{2}{\rho }_{o}{\left(\displaystyle \frac{r}{{r}_{o}}\right)}^{-1.5}{dr}=\displaystyle \frac{4\pi }{1.5}{r}_{o}^{3}{\rho }_{o},\end{eqnarray} \tag{ 20 }$$

where r_o is the outer radius and ρ_o is the density at r_o. Therefore, the mass inflow rate at r_o can be estimated using

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=4\pi {r}_{o}^{2}{\rho }_{o}{v}_{\mathrm{in}}=1.5{{Mv}}_{\mathrm{in}}/{r}_{o}.\end{eqnarray} \tag{ 21 }$$

We take the total clump mass (∼200 M_⊙) for M and clump radius (0.28 pc) for r_o. We assume an infall velocity v_in of 0.32 km s⁻¹ obtained from the HCO⁺ (1–0) measurement because the mean volume density of "c8" is closer to the critical density of HCO⁺ (1–0). The inferred mass inflow rate (${\dot{M}}_{\mathrm{in}}$) is thus ∼4 × 10⁻⁴ M_⊙ yr⁻¹. This mass inflow rate is consistent with those measured in other high-mass star-forming clumps (Wu et al. 2009, 2014; Sanhueza et al. 2010; Liu et al. 2011a, 2011b, 2013a, 2013b, 2016b; Ren et al. 2012; Peretto et al. 2013; Qin et al. 2016; Yuan et al. 2018). The clump mass of "c8" also exceeds the empirical threshold ($M\gt 870\,{M}_{\odot }{(r/\mathrm{pc})}^{1.33}$) for high-mass star-forming clumps discovered by Kauffmann & Pillai (2010). In addition, the clump mass of "c8" is also comparable to the masses of high-mass starless clumps with similar radii cataloged by Yuan et al. (2017). All of these conditions indicate that "c8" has the potential to form high-mass stars.

The collapsing massive starless clump candidate "c8" may represent the very initial conditions for high-mass star formation and deserves more detailed studies at higher angular resolution. Indeed, searching for the existence or absence of high-mass prestellar cores, as has been done in other massive starless clump candidates (Beuther et al. 2013; Sanhueza et al. 2013, 2017; Tan et al. 2013; Liu et al. 2017; Contreras et al. 2018), is of great importance given that in "c8" the magnetic field and infall speed at large scales are now both known, unlike for the other studies.