Formation of the SDC13 Hub-filament System: A Cloud–Cloud Collision Imprinted on the Multiscale Magnetic Field

In order to identify the ridges of filamentary structures in the JCMT 850 μm Stokes I image, we adopt the DisPerSE algorithm (Sousbie 2011). We use a contrast threshold of 20 mJy beam⁻¹ (∼5σ) for the filament identification, and exclude the identified filaments with lengths shorter than 1 arcmin to ensure the significance. The identified filaments are plotted in Figure 4(a). The four longest filaments, the northeast (NE), northwest (NW), north (N), and south (S) filaments, are shown to converge to the central hub, forming a "Y"-shaped hub-filament system. This is consistent with the four filaments identified in Peretto et al. (2014) from IRAM 30 m MAMBO 1.2 mm dust continuum data and also with the NH₃ observations in Williams et al. (2018). Several shorter filaments are merging into the longer filaments via local convergent points. Two bridging filaments seem to connect the filament NE and the eastern hub.

Zoom In Zoom Out Reset image size

We use the python package FilChap (Suri et al. 2019) to estimate the widths of the identified filaments. We fit individual radial intensity profiles extracted at each pixel position along the filament ridges. A bootstrap method is applied to fit these radial profiles with a Gaussian function, and the uncertainties are estimated using a Monte Carlo approach to generate 100 simulated profiles based on the observed intensities and uncertainties. Those fits with a width smaller than three times the uncertainties are excluded from the further analysis. The fitted Gaussian widths (σ_F ) are converted to the FWHM width (ΔF) via ${\rm{\Delta }}F=\sqrt{8{ln}2}{\sigma }_{F}$ which are shown in Figure 4(b). We note that the intensity profiles along the filament N overlap with filament NE and NW leading to highly uncertain fitting results. They are, therefore, not shown in Figure 4(b). Generally, the filament widths grow with the local intensity, from ∼0.3 pc in the diffuse regions to ∼0.5 pc near the dense hub and clumps. We note that the smaller filament widths in the diffuse regions might be consistent with a universal 0.1 pc filament width discovered by the Herschel Gould Belt survey (e.g., Arzoumanian et al. 2011), although our resolution of 0.24 pc is insufficient to clearly resolve a 0.1 pc width.

In order to investigate the occurrence of star formation in the SDC13 hub-filament system, we overplot the starless and protostellar cores from Peretto et al. (2014) on the filaments in Figure 4. Without any exception, all the cores are located either on filament spines or in filament convergent points. The protostellar cores are mostly found near convergent points, while more starless cores are distributed along the filaments. This distribution is suggestive of star formation taking place not only in the central hub, but also along the converging filaments. Moreover, local convergent points seem to have a higher probability of triggering star formation. These features have also been seen in Williams et al. (2018) where the filaments are identified from the NH₃ data. This suggests that NH₃ and continuum data likely trace the same filamentary structures.

In order to investigate whether gravity influences the formation of the converging filaments, we estimate the projected gravitational vector field from the JCMT 850 μm continuum data. As the large-scale emission beyond about 3' is filtered out by the JCMT POL-2 observations, we focus here on the local gravitational field within the denser SDC13. The inter-cloud-scale gravitational field, traced by Herschel observations, is discussed in Section 5.1.

Following the development of the polarization-intensity gradient technique in Koch et al. (2012a, 2012b), the local projected gravitational force acting at a pixel position ( F _G,i ) can be expressed as the vector sum of all gravitational pulls generated from all pixel positions over the map (Wang et al. 2020) as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{G,i}={{kI}}_{i}\sum _{j=1}^{n}\displaystyle \frac{{I}_{j}}{{r}_{i,j}^{2}}\hat{r},\end{eqnarray} \tag{ 3 }$$

where k is a factor accounting for the gravitational constant and conversion from emission to total column density. I_i and I_j are the intensity at the pixel position i and j, and n is the total number of pixels within the area of relevant gravitational influence. r_i,j is the plane-of-sky-projected distance between the pixel i and j, and $\hat{r}$ is the corresponding unity vector. The above equation assumes that the intensity distribution is a fair approximation for the distribution of the total mass, and that the mass components in SDC13 are roughly at the same distances. A constant factor k = 1 is used because we will only utilize the directions of the local gravitational forces and not their absolute magnitudes. When calculating the local gravitational field, a threshold of 20 mJy beam⁻¹ (∼5σ) is introduced below which any gravitational influence is neglected. This is justified because any gravitational force originating from diffuse and extended structures tends to be rather symmetrical which means that any already small gravitational pulls will largely cancel out.

Figure 5 displays the local gravitational vector field in SDC13. Looking for distinct directions in this vector field, local gravity visually appears as a combination of two modes: (1) pulling to the filament ridge (local vectors prevailingly orthogonal to the filament's ridge) and (2) pulling to a converging center (local vectors rather along the filament's ridge and local vectors pointing azimuthally symmetrically to a mass center.

Zoom In Zoom Out Reset image size

The relative importance of these two modes likely determines the orientation offsets between local gravity and filaments. In the NE and S filament, local gravity seems to be more efficiently dragging material orthogonally to the filament ridges, except for the two local converging centers at the tips of the filaments. In the NW filament, the gravitational force pulling toward the central hub is likely comparable to the force toward the filament ridge, and hence an orientation offset of ∼30°–60° can be seen between local gravity and filament. In the hub center, the gravitational field is mainly pointing toward the most massive protostellar core within the central hub (also known as MM1), and the orientation offsets between local gravity and filaments are ∼40°–90°. We note that this analysis is limited by the observational resolution, and the gravitational field originated from possible structures at scales smaller than our resolution cannot be directly probed.

To analyze the velocity structure in SDC13, we derive local velocity gradients from the NH₃ (1, 1) centroid line-of-sight (LOS) velocity map in Williams et al. (2018). Probing the velocity structures at the same physical scale as the magnetic field traced by our polarization data requires smoothing the NH₃ data to a pixel size of 7'' from the original beam size of 40 × 28 (Figure 6). We estimate the centroid velocity by using the fitting schemes in the CLASS software. ¹² The local velocity gradient is calculated by fitting the 2D spatial distribution of the centroid velocities within 3 × 3 pixel (Goodman et al. 1993) with

$$\begin{eqnarray}&&{V}_{\mathrm{LSR}}={c}_{x}x+{c}_{y}y+{c}_{0},\end{eqnarray} \tag{ 4 }$$

where x and y are the positions of each pixel, and c₀, c_x , and c_y are free parameters representing the first-order expansion of the velocity field.

Zoom In Zoom Out Reset image size

The resulting local velocity gradient field is plotted in Figure 6. The velocity gradients are along an east–west direction near the hub center. They become more complex in the filament regions, seemingly either pointing toward the local clumps along the filaments, or being perpendicular to the filament ridges. A more thorough statistical analysis is necessary to understand whether these trends are significant or not.

In this section we perform systematic statistical analyses to investigate how the physical parameters act in the hub-filament system (HFS) and to reveal possible correlations among the orientations/directions of filaments (F), magnetic field (B), gravitational force (G), and gas velocity gradient (VG). The magnitudes of these parameters are excluded from our further analysis here, because our focus is on understanding first the possible correlations in orientations. Moreover, the complete three-dimensional information to possibly correct the projected magnitudes is observationally not accessible.

For the selected filaments within SDC13, we estimate a filament orientation at every pixel in the following way. For the ith pixel along the ridge of a filament, we fit the positions of the (i − 2)th to the (i + 2)th consecutive pixels along the filament with a straight line to estimate the local filament orientation. In this way, the median fitting error in the local orientation of a filament is ∼2°.

Figure 7 presents the histograms of the local orientations of filaments (F), magnetic field (B), gravity (G), and velocity gradients (VG) for the entire SDC13, and separated into high- and low-density regions. In all the regions, F, B, G, and VG show rather clear single or multiple peaks in their distributions. This is clearly different from random, i.e., uniform distribution in orientations. Noticeably, the magnetic field orientations display a more pronounced peak in high-density regions. Other parameters also appear to change with density, though in less definite ways.

Zoom In Zoom Out Reset image size

To further isolate possible correlations, we perform an all-pairwise comparison of the four parameters (F, B, G, VG), looking at their relative orientations at different densities. For each measurement of one parameter, we select a nearest measurement of another parameter within a radius of 14'' (one independent beam). These two measurements are then defined as one associated pair for the two parameters. A relative orientation (ΔPA) is calculated for each associated pair. We perform a one-sample Kolmogorov–Smirnov (KS) test to the resulting six ΔPA distributions in each region, to examine whether these distributions differ from a uniform distribution, as the null hypothesis. For each distribution, a probability (p-value) that the observed distribution can be drawn from the hypothetical uniform distribution is calculated, and a threshold of p = 0.05 (95% confidence interval) is used to reject the null hypothesis. We note that a p-value higher than 0.05 does not automatically indicate a random distribution, but could also result from an insufficient number of data points. Table 1 lists the features of those parameter pairs with nonrandom distributions. The complete histograms for all the parameter pairs, additionally divided into hub and filament regions, are given in Section B.

As the dominating physical process within a cloud might evolve with local density (Soler et al. 2013; Planck Collaboration et al. 2016; Wang et al. 2020), we have introduced a density threshold of 2 × 10²² cm⁻² (250 mJy beam⁻¹, assuming a dust temperature of 27.2 K; Williams et al. 2018), defining low- and high-density regions. This threshold is motivated as the differences in the resulting ΔPA distributions become most obvious. The histograms of the relative orientations for these two regimes are shown in Figure 8. The most significant difference between these two regimes is that filaments tend to be more aligned with the magnetic field in low-density regions, but become more perpendicular in high-density regions (Figure 8(a)). Similar trends are seen in numerous molecular clouds based on PLANCK data (Planck Collaboration et al. 2016), though with a substantially larger beam these data are rather probing the transition between the diffuse interstellar medium and molecular clouds.

Zoom In Zoom Out Reset image size

In addition to the magnetic field trend, both gravity and velocity gradients appear to fall into a similar range of relative orientations, ∼20°–75°, with respect to the filaments (Figures 8(b) and (d)). This offset angle is likely representing the combination of the two modes in the local gravitational force map, as pointed out in Section 4.1.2. Gravity can locally both pull the gas toward a filament ridge (leading to a more perpendicular gravity-filament configuration) and also radially toward the central hub (leading to a more parallel gravity-filament configuration). The offset angle is likely determined by the relative importance of these two modes. Thus, it can still vary over a range. Similarly, the offset angles between the velocity gradients and the filaments might be linked to these two modes. We, however, note that the G–F offset angle is not necessarily the same as the VG–F offset angle, because a filament's radial and longitudinal collapsing timescale can be different, depending on the density and the geometry of a filament (Vázquez-Semadeni et al. 2009, 2017).

An additional finding is that the velocity gradient in the low-density regions is more correlated with gravity (being more perpendicular), while in the high-density regions it is more correlated with the magnetic field (being more parallel; Figures 8(e) and (f)). This growing alignment between magnetic field and velocity gradient might point at a role of the magnetic field in guiding gas motions as density increases. This possibly results from the enhanced magnetic field in the central hub, where the magnetic field shows the largest field strength and a small angular dispersion (Figure 3; Table 2). In contrast to that, the magnetic field in the low-density regions, where overall the angular dispersion is larger, is more likely to be locally distorted by turbulence, core fragmentation, or external pressure. As a consequence, this more complex and less organized magnetic field morphology presents a largely random orientation with respect to the velocity gradient. Therefore, the magnetic field is also less capable of constraining gas motion.

Table 2. Magnetic Field Strengths and Mass-to-flux Ratios in SDC13

Regions	${n}_{{{\rm{H}}}_{2}}$	${N}_{{{\rm{H}}}_{2}}$	σv,NT	δϕ	Bpos(DCF)	λ (DCF)	Bpos(ST)	λ (ST)
	cm−3	cm−2	(km s−1)	(deg)	(μG)		(μG)
Hub	(4.2 ± 0.4) × 104	(3.2 ± 0.3) × 1022	0.41 ± 0.05	18.5 ± 1.0	94 ± 5	0.87 ± 0.05	75 ± 2	1.08 ± 0.03
Filament NE	(2.5 ± 0.2) × 104	(1.8 ± 0.1) × 1022	0.32 ± 0.05	35.5 ± 1.4	31 ± 1	1.50 ± 0.05	34 ± 1	1.34 ± 0.03
Filament NW	(1.1 ± 0.2) × 104	(1.3 ± 0.1) × 1022	0.22 ± 0.05	15.5 ± 2.3	34 ± 5	0.95 ± 0.14	25 ± 2	1.29 ± 0.10
Filament S	(2.9 ± 0.3) × 104	(1.5 ± 1.0) × 1022	0.32 ± 0.05	20.2 ± 1.2	58 ± 4	0.66 ± 0.03	49 ± 1	0.79 ± 0.02

Note. Magnetic field strengths and mass-to-flux ratios derived from the DCF and ST method. The uncertainties listed here are obtained from propagating the observational uncertainties through the corresponding equations. Possible additional systematic uncertainties due to, e.g., the unknown dust opacity, are not included.

Download table as:  ASCIITypeset image

Finally, we acknowledge that some of the possibly emerging correlations and trends discussed here can be affected by unknown projection and integration effects. A more complete statistical analysis with a larger sample is needed to fully establish these trends.

The DCF method assumes that kinetic and magnetic energy are in equipartition. Hence, the level of magnetic field perturbation, traced by the polarization angular dispersion δϕ, is the result of transverse incompressible turbulent Alfvén waves, traced by the LOS nonthermal velocity dispersion σ_v,NT. Under this assumption, the plane-of-sky magnetic field strength (B_pos) can be estimated as

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}=Q\,\sqrt{4\pi \rho }\displaystyle \frac{{\sigma }_{v,\mathrm{NT}}}{\delta \phi },\end{eqnarray} \tag{ 5 }$$

where ρ is the gas volume density, and Q is a factor accounting for complex magnetic field and inhomogeneous density structures. Ostriker et al. (2001) suggested that Q = 0.5 yields a good estimation of the magnetic field strength in the plane of sky if the magnetic field angular dispersion is less than 25°.

As the role of the magnetic field is likely different in the filament and the hub regions of SDC13, we estimate the magnetic field strength in each region separately. Polarization segments are selected for each region, and a polarization angular dispersion is calculated. In order to obtain a polarization dispersion originating from turbulent perturbation solely, i.e., without any confusion from larger-scale magnetic field features, we calculate differences between polarization position angles using only nearest pixel pairs. In this way, the polarization angular dispersion is derived as

$$\begin{eqnarray}&&\delta \phi =\sqrt{\displaystyle \frac{1}{N-1}\sum _{i=1}^{N}\delta {{PA}}_{i}^{2}},\end{eqnarray} \tag{ 6 }$$

where δPA is the absolute difference between two polarization position angles for every possible nearest-polarization pair. N is the total number of pairs, and the factor 1/(N − 1) is to debias the population standard deviation estimator. The calculated δϕ for each region is listed in Table 2.

To estimate the mean volume density in each region, we first construct a column density map using the POL-2 850 μm continuum data assuming a constant temperature of 12.7 K, adopted from the mean NH₃ rotational temperature (Williams et al. 2018), and a dust opacity κ of 0.012 cm² g⁻¹ at 850 μm (Hildebrand 1983). The total mass of a region is then obtained by integrating the column density over the selected region. The volume of the hub region is estimated assuming a sphere with a FWHM diameter of 1.0 ± 0.1 arcmin, obtained from a 2D Gaussian fit to the hub. The volume of a filament region is derived adopting a cylinder with the filament's radius of 0.5 pc, given in Section 4.1.1. The mean volume densities are then calculated from the total masses and their respective volumes.

In order to obtain a mean nonthermal velocity dispersion, we first average the observed NH₃ line widths in each region. Assuming a gas kinematic temperature (T_kin) of 12.7 K, the thermal velocity dispersion for NH₃ is $\sqrt{\tfrac{{k}_{B}{T}_{\mathrm{kin}}}{{m}_{{\mathrm{NH}}_{3}}}}=0.10$ km s⁻¹. The thermal velocity dispersion is then removed from the observed line width to obtain the nonthermal velocity dispersion (σ_v,NT) as

$$\begin{eqnarray}&&{\sigma }_{v,\mathrm{NT}}^{2}={\sigma }_{\mathrm{obs}}^{2}-\displaystyle \frac{{k}_{B}{T}_{\mathrm{kin}}}{{m}_{{\mathrm{NH}}_{3}}},\end{eqnarray} \tag{ 7 }$$

where σ_obs is the observed NH₃ velocity dispersion, and ${m}_{{\mathrm{NH}}_{3}}$ is the molecular weight. The resulting magnetic field strengths for the four regions in SDC13 are listed in Table 2.

The mass-to-flux criticality (λ_obs) is commonly used to evaluate the relative importance between magnetic field and gravity (Nakano & Nakamura 1978), and calculated as

$$\begin{eqnarray}&&{\lambda }_{\mathrm{obs}}=2\pi \sqrt{G}\displaystyle \frac{\mu {m}_{{\rm{H}}}{N}_{{{\rm{H}}}_{2}}}{{B}_{\mathrm{pos}}},\end{eqnarray} \tag{ 8 }$$

where μ = 2.33 is the mean molecular weight per H₂ molecule, G is the gravitational constant, and ${N}_{{{\rm{H}}}_{2}}$ is the molecular hydrogen column density. To correct for the unknown projection effect, Crutcher et al. (2004) suggest that a statistical average factor of 1/3 can be used to better estimate the mass-to-flux ratio of oblate spheroid cores, flattened perpendicular to the magnetic field. The corrected mass-to-flux ratio (λ) becomes

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{\lambda }_{\mathrm{obs}}}{3}.\end{eqnarray} \tag{ 9 }$$

The resulting mass-to-flux ratios are around 0.7–1.5 for all regions in SDC13. This suggests that both the hub and the filament regions in SDC13 are about transcritical. Different correction factors, rather than the factor 1/3 used in Equation (9), are suggested for different cloud and magnetic field geometries, e.g., π/4 for a spherical cloud (Crutcher et al. 2004), and 3/4 for a prolate spheroid elongated along the magnetic field (Planck Collaboration et al. 2016). Adopting these numerical factors, the mass-to-flux ratios become transcritical to supercritical with values ranging from 1 to 3. We emphasize that the field strengths and mass-to-flux ratios derived here are based on values averaged over the selected regions. Therefore, they have to be interpreted as average global properties. Local higher-density regions, within regions with a criticality of 0.5 to 3, are likely (highly) supercritical, consistent with the presence of star formation activity.

Skalidis & Tassis (2021) point out that observational results show that turbulence in the ISM is anisotropic and that nonAlfvénic (compressible) modes may be important. Hence, they propose a new method to estimate the magnetic field strength in the ISM considering these compressible modes, leading to

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}=\sqrt{2\pi \rho }\displaystyle \frac{{\sigma }_{v,\mathrm{NT}}}{\sqrt{\delta \phi }}.\end{eqnarray} \tag{ 10 }$$

The resulting magnetic field strengths and the corresponding mass-to-flux ratios using the ST method are listed in Table 2. Field strengths estimated from the ST method are similar to the ones derived from the DCF method within about 30%. The mass-to-flux ratios are mostly around unity.

In this section, we apply the virial theorem to evaluate the relative importance between gravitational, magnetic, and kinetic energy in the central hub. As the gravitational potential and pressure of a cylinder follow a formalism different from the one for a sphere, the virial analysis for the filament regions is presented in Section 4.2.4. In Lagrangian form, the virial theorem can be written as

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\ddot{I}=2({ \mathcal T }-{{ \mathcal T }}_{s})+{ \mathcal M }+{ \mathcal W }\end{eqnarray} \tag{ 11 }$$

(e.g., Mestel & Spitzer 1956; McKee & Ostriker 2007), where I is a quantity proportional to the trace of the inertia tensor of a cloud. The sign of $\ddot{I}$ determines the acceleration of the expansion or contraction of the spherical cloud. The term

$$\begin{eqnarray}&&{ \mathcal T }=\displaystyle \frac{3}{2}M{\sigma }_{\mathrm{obs}}^{2}\end{eqnarray} \tag{ 12 }$$

is the total kinetic energy, where M is the total mass, and σ_obs is the observed total velocity dispersion. We neglect the surface kinetic term ${{ \mathcal T }}_{s}$ because we aim to study the self-stability of each region. Nevertheless, we note that the presence of any external pressure could suppress the support and enhance the cloud's instability. The magnetic energy term, without any force from an external magnetic field, is

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{1}{2}{{Mv}}_{A}^{2},\end{eqnarray} \tag{ 13 }$$

where ${v}_{A}=B/\sqrt{4\pi \rho }$ is the Alfvén velocity, and ρ is the mean density. We note that the magnetic field morphology is not explicitly accounted for in this magnetic energy term. As the DCF and SF method only constrain the plane-of-sky magnetic field component, we use the statistical average to correct and estimate the total magnetic field strength as B = (4/π)B_pos (Crutcher et al. 2004). The term

$$\begin{eqnarray}&&{ \mathcal W }=-\displaystyle \frac{3}{5}\displaystyle \frac{{{GM}}^{2}}{R}\end{eqnarray} \tag{ 14 }$$

is the gravitational potential of a sphere with a uniform density ρ and a radius R.

The resulting energy ratios for the central hub in SDC13 are a kinetic-to-gravitational energy $| { \mathcal T }/{ \mathcal W }| =0.04\pm 0.02$ and a magnetic-to-gravitational energy $| { \mathcal M }/{ \mathcal W }| =0.05\pm 0.01$. This suggests that globally gravity is dominating over both kinetic and magnetic energy in the central hub. Moreover, the magnetic energy is comparably important as the kinetic energy, with a derived Alfvénic Mach number of 0.6 ± 0.1. The combined ratio $| \tfrac{2{ \mathcal T }+{ \mathcal M }}{{ \mathcal W }}|$ is 0.09 ± 0.02, smaller than unity, suggesting that even the combined kinetic and magnetic energy cannot support the system, and hence the central hub is contracting globally.

Williams et al. (2018) show that the filaments in SDC13 are on average thermally supercricital, suggesting that the thermal energy is insufficient to support the filaments against gravitational collapse and fragmentation. In this section, we further investigate how the filament criticality varies spatially, considering the support from thermal, nonthermal, and magnetic energy. Based on the Virial theorem, the stability of filaments is commonly evaluated using the critical line density, (or critical line mass), ${M}_{\mathrm{line},\mathrm{critical}}=2({c}_{s}^{2}+{\sigma }_{v,\mathrm{NT}}^{2})/G$ (Fiege & Pudritz 2000), which considers the balance between gravity and the support from both the thermal and turbulent energy. In order to reveal how the filament stability in SDC13 spatially changes, we use the local column density, filament width, and NH₃ velocity dispersion from the previous sections to estimate the filament criticality pixel by pixel. The local filament line density is calculated as M_line = ΣW_F , where ${\rm{\Sigma }}=\mu {m}_{{{\rm{H}}}_{2}}{N}_{{{\rm{H}}}_{2}}$ is the filament peak surface density assuming a mean molecular weight μ of 2.33, and W_F is the deconvolved filament FWHM width. The filament criticality is defined as M_line/M_{line,critical}. Hence, a criticality larger than unity is suggestive of a collapsing and fragmenting filament.

The calculated filament criticality is visualized in Figure 10. We generally find that most filaments in SDC13 are supercritical. This is consistent with the presence of the numerous starless and protostellar cores along the filaments (Figure 4). A few filaments show subcritical locations in the outer diffuse areas. These are typically at the tips of the filaments. Around the dense cores in the filament NE and S, the local filament criticality is clearly increasing from subcritical in the outer areas to supercritical in the dense cores. As the change in filament criticality is mainly due to the increasing local column density, this likely indicates that these filaments transition from subcritical to supercritical via accumulating mass, possibly through converging filaments, which then destabilizes the filaments. Subsequently, this will trigger core fragmentation and (future) star formation, as witnessed by the presence of protostellar and starless cores.

Zoom In Zoom Out Reset image size

Additionally to the support from gas kinetic energy, including thermal and nonthermal, magnetic fields can play a role in stabilizing filaments. However, the exact magnetic support depends on the morphology of the magnetic field within a filament, for which our current data do not have sufficient resolution. Fiege & Pudritz (2000) model magnetized filaments based on the virial theory and suggest that poloidal-dominated fields help supporting filaments against gravity while toroidal-dominated fields destabilize filaments. The critical linear density of a magnetized filament is ${M}_{\mathrm{line},\mathrm{vir}}^{\mathrm{mag}}={M}_{\mathrm{line},\mathrm{critical}}\times {(1\pm { \mathcal M }/{ \mathcal W })}^{-1}$, where ${ \mathcal M }$ is the magnetic energy per unit length (with a positive sign for poloidal fields and a negative sign for toroidal fields), and ${ \mathcal W }$ is the gravitational energy per unit length. Adopting a $| { \mathcal M }/{ \mathcal W }|$ of 0.05 ± 0.01, as derived in Section 4.2.3, the filament criticality can change within ±5%. Hence, even if the magnetic field is poloidal-dominated, the filaments around the dense cores likely remain supercritical, and only small stretches of the filament's outer diffuse regions might transition from supercritical to subcritical. Therefore, the overall finding of the filaments generally being supercritical in SDC13 remains valid even in the presence of magnetic fields, confirming that the conclusion in Williams et al. (2018) is still valid even considering the magnetic and turbulent energy.

Finally, we note that the underlying assumption to estimate the filament criticality is that the observed velocity dispersion trace the internal turbulent motions. This assumption might not be fully valid within the central hub, where multiple filaments and cores partially overlap, and resulting complex velocity structures might all contribute to the observed velocity dispersion. Williams et al. (2018) suggest that the increasing velocity dispersion toward the dense regions in SDC13 is possibly driven by the gravitational fragmentation and infall motion. Such motion is guided by the local gravity and does not support the filaments as efficiently as isotropic turbulence. This is possibly the reason why the growing filament criticality, as seen toward the dense cores in filament NE and S, is not that evident in the central hub.

In order to reveal the 10 pc scale magnetic field, Figure 11 shows the 353 GHz PLANCK magnetic field segments, with a beam size of 5', overlaid on our JCMT POL-2 polarization data. The large-scale magnetic field traced by PLANCK reveals a partially spiraling and converging morphology, possibly also resembling two incoming wings from the east and west, pointing toward the center of SDC13. This large-scale converging magnetic field pattern tends to locally align with the filament NE and to be perpendicular to filament S and NW. Similarly, the large-scale magnetic field tends to be either parallel or perpendicular to the nearby compact clouds detected in our POL-2 map (also see Section C).

Zoom In Zoom Out Reset image size

Magnetic fields either parallel or perpendicular to filamentary clouds have been found from statistics of PLANCK polarization (Planck Collaboration et al. 2016) and starlight polarization data (Li et al. 2013). A number of models have been proposed to explain the origin of this configuration, such as magnetic-field-channeled turbulence/shock compression (Nakamura & Li 2008; Inoue & Inutsuka 2009; Inoue et al. 2009; Chen et al. 2020) or filament–filament collision (Nakamura et al. 2014). In these models, magnetic fields are important in channeling and guiding the mass accretion and subsequent cloud collapse. Following these scenarios, the observed large-scale converging magnetic field morphology in Figure 11 is possibly guiding the large-scale gas flows converging toward SDC13. A similarly converging magnetic field, also guiding converging accretion flows, has been seen in another hub-filament system, G33.92+0.11 (Wang et al. 2020), although at a smaller parsec scale.

In order to reveal how the magnetic field is impacted by the possible converging events, we calculate the relative orientations between the small- and large-scale magnetic fields traced by PLANCK and POL-2 as displayed in Figure 11 with the colored points. We find a trend that the relative orientations vary from ∼0° to ∼90° from north to south within SDC13. Similar variations also occur in the nearby compact clouds where the small-scale magnetic field is similarly oriented as the large-scale magnetic field on one side, but then differs in its orientation on the other side of the cloud. As these large changes in relative orientations occur smoothly over the entire parsec-scale cloud, it is unlikely that they originate from the smaller-scale 0.1–0.001 pc core fragmentation or related star formation processes.

The above discussed variation from large- to small-scale magnetic field morphology can be explained by a converging flow scenario. The induced shocks at the flow colliding layers could compress the local magnetic field, and lead to a magnetic field parallel to the compressing layer or perpendicular to the converging flow (Peretto et al. 2012; Inoue et al. 2018). In SDC13, the large-relative-orientation areas in Figure 11, covering the filament S and the southern part of the central hub, are perpendicular to the large-scale magnetic field. This is consistent with the expectation if the southern part of SDC13 originates from a flow collision guided by the large-scale magnetic field. In the nearby compact clouds, the large-relative-orientation areas are commonly either parallel or perpendicular to the large-scale magnetic field (see Section C for a brief description on the individual compact clouds).

In order to recover the large-scale density structures filtered out by the JCMT POL-2 observations, we used the Herschel archival 5-band continuum data within an area of 2 × 2 deg² around SDC13 to construct a column density map. We smoothed the 70–500 μm continuum data to a resolution of 35'', and fit them with a graybody function using a dust opacity κ of 0.012 cm² g⁻¹ and a β of 2 (Hildebrand 1983). A zoom-in column density map around SDC13 is shown in Figure 12. This map reveals that SDC13 is embedded in a large-scale north–south filament with a width of ∼8 pc. The PLANCK magnetic field outside of this filament appears to be relatively uniform and mostly perpendicular to this filament, while the converging magnetic field pattern, as described in Section 5.1.1, is within this filament. This suggests that the converging pattern is likely linked to and driven by the kinematics of the outer larger-scale filament.

Zoom In Zoom Out Reset image size

Simulations predict that the formation of filaments threaded by perpendicular magnetic fields can twist the magnetic field and form a toroidal magnetic field morphology (e.g., Li & Klein 2019). Furthermore, predicted helical magnetic fields have been observed in a number of clouds (e.g., Tahani et al. 2018). The spiral-like converging magnetic field within the SDC13-encompassing large-scale filament might also be explained by an inclined toroidal or semi-helical magnetic field twisted by the large-scale filament. Adopting such a scenario, this encompassing large-scale filament threaded by the twisted magnetic field might be replenished by gas flows that are in turn guided by the magnetic field to a central converging region. This mechanism has also been suggested by joint gas kinematics and polarization studies toward other large-scale filamentary systems. Striations near filaments have been commonly found converging toward filament crests along plane-of-sky magnetic fields, e.g., Taurus B211/213 (Palmeirim et al. 2013), Musca (Bonne et al. 2020a, 2020b). Possible helical magnetic fields have also been detected surrounding large-scale filaments, e.g., Orion A, California, and Perseus (Tahani et al. 2018). Hence, the collision of converging flows might be the origin of the compact clouds embedded within the large-scale encompassing filament.

To test whether the large-scale gravitational field is consistent with the above outlined scenario, we calculate the projected gravitational field using the 2 × 2 deg² column density map, following our approach in Section 4.1.2 (Figure 12). The gravitational field outside of the encompassing large-scale filament appears to first converge toward the filament ridge along the large-scale magnetic field from northeast and west. Getting closer to the ridge of the large-scale filament, the gravitational field is turning toward south to become predominantly aligned with the orientation of the large-scale filament. The gravitational field also gradually becomes aligned with the magnetic field as it approaches the filament from the eastern and western side (between the first and second contour). This is consistent with the above scenario where the magnetic field is channeling the accretion flows. However, we also note that the gravitational field in the diffuse areas, far from the filament ridge, can show an orientation very different from the magnetic field. This might be explained by a magnetic field being more important in binding the gas than gravity in such low-density regions. This is possible because the typical mass-to-flux ratios of clouds with densities (N ${}_{{{\rm{H}}}_{2}}$) less than 10²¹ cm⁻² are likely still magnetically subcritical (Crutcher 2012).

In order to trace the gas possibly flowing along the large-scale magnetic field, we use IRAM 30 m C¹⁸O (1–0) data (G. M. Williams et al. 2022, in preparation) to identify the velocity components that are likely associated with SDC13. This data set (Project ID: 024-13, PI: Peretto) has a beam size of 246, a spectral resolution of 0.13 km s⁻¹, and an rms noise of 0.17 K. Figure 13(a) shows the identified three components, 5–20 km s⁻¹ (blue), 32–40 km s⁻¹ (green), and 42–58 km s⁻¹ (red). The green component clearly displays a "Y-shape" morphology, spatially overlapping with the "Y-shape" structure of SDC13 seen in the 850 μm continuum. The blue component reveals four filamentary arms converging to the dense core embedded in the filament S. One of these arms is aligned with the arc-like magnetic field pattern and also the eastern boundary of the central hub. The red component is connecting the eastern with the central hub through a bridge.

Zoom In Zoom Out Reset image size

The spatial arrangement of these C¹⁸O components supports our earlier outlined scenario where the observed changes from large- to small-scale magnetic field originate from cloud–cloud collisions. The "U-shape" morphology of the blue component, aligned with the "U-shape" magnetic field structure, is identical to a shock-compressed layer where the magnetic field is bent by the collision, as predicted in cloud–cloud simulations (Inoue & Fukui 2013; Inoue et al. 2018). The red component presents a possible incoming flow moving toward SDC13 from the northeast. Additionally, it probably also acts as a channel from the eastern hub.

To further understand the origin of the C¹⁸O velocity components, we searched for even larger-scale "reservoirs". We, therefore, looked for giant molecular clouds (GMCs) identified from Atacama Pathfinder EXperiment (APEX) ¹³CO (2-1) data in the SEDIGISM catalog (Duarte-Cabral et al. 2021; Schuller et al. 2021). GMCs were selected based on overlapping with SDC13 in position-position–velocity. To extract only the major cloud components, we excluded GMCs with an angular size smaller than 5'. The final selected GMCs are SDG013.178-0.0950 (blue), SDG012.840-0.2041 (green), and SDG013.222+0.0076 (red), which have a velocity range of 12–20 km s⁻¹, 35–38 km s⁻¹, and 47–59 km s⁻¹, respectively, consistent with the three velocity components found in the SDC13 IRAM 30 m C¹⁸O data. Figure 13(b) presents the integrated intensity maps of these three selected GMCs in different colors. We note that another GMC, SDG013.098-0.0821, also overlaps with SDC13. However, its corresponding C¹⁸O component, in the range 20–28 km s⁻¹, is associated with the compact cloud SDC13.121-0.091 (south of SDC13, see Figure 1). Hence, we do not further discuss it here.

SDG012.840-0.2041 (green) is consistent with both the large-scale north–south filament revealed by the Herschel column density map, and also with the main structure of SDC13 in the C¹⁸O data. Its velocity range is consistent with the NH₃ and C¹⁸O velocity component tracing the main structure of SDC13. SDG013.178-0.0950 (blue) and SDG013.222+0.0076 (red) appear with an extended morphology likely winding around the north–south filament along the magnetic field. The velocity difference between these two GMCs and SDG012.840-0.2041 is about 10–20 km s⁻¹, which is considered a typical range enabling cloud–cloud collisions and the formation of a massive cluster (∼10–40 km s⁻¹, Inoue et al. 2018; Cosentino et al. 2019; Dobbs & Wurster 2021; Fukui et al. 2021).

The positions of the compact clouds detected in our JCMT POL-2 data (Figure 1) are mostly falling onto the conjunction areas between two or three GMCs. SDC13 is embedded in SDG012.840-0.2041, and possibly compressed by the other two GMCs from the south-western and eastern side. The two compact clouds (SDC13.121-0.091 and SDC13.123-0.157) in the south of SDC13 are located at the boundaries of these GMCs, with their major axes along the boundaries. A series of small compact clouds north of SDC13 (SDC13.225-0.004 and SDC13.246-0.081) is located along the extending structures of SDG013.222+0.0076 overlapping with SDG012.840-0.2041, and connecting to the two compact clouds (SDC13.190-0.105, SDC13.198-0.135) east of SDC13.

In summary, this spatial consistency suggests that the interaction between these GMCs is highly correlated with the formation of the compact clouds. We, however, note that the sensitivity of the molecular line data is not yet sufficient to also detect the even diffuser gas that very likely is present within the velocity gaps between these GMCs. Such bridging structures connecting velocity components in the position–velocity diagrams are commonly considered as evidence of cloud–cloud collision (e.g., Habe & Ohta 1992). Future high-sensitivity observations are needed to understand and possibly model the interaction between these GMCs in more detail.