© The Authors 2025

1 Introduction

Star formation predominantly occurs in a clustered mode (Lada et al. 2003), with hub-filament systems (HFSs) (Myers 2009; Peretto et al. 2014) now widely recognized as the progenitors of stellar clusters (Schneider et al. 2012; Kumar et al. 2020). HFSs naturally encompass a broad range of filament densities (Hacar et al. 2023), as they gather gas from the surrounding diffuse cloud via low-density filaments (Kirk et al. 2013; Palmeirim et al. 2013; Schisano et al. 2014), while simultaneously building high-density filament networks within the central hub (Myers 2009; Kumar et al. 2022; Hacar et al. 2023). It has been established that most massive stars form within hubs (Kumar et al. 2020), as the high-density filament network within the hub (Kumar et al. 2022) creates the conditions necessary to achieve the highest critical densities required for massive stellar seeds to form. These seeds are further nourished to grow in mass (Kumar et al. 2020, 2022), as the hub serves as the central convergence point for gas drawn from the surrounding cloud (Peretto et al. 2013, 2014).

Myers (2009) defined the hub as the central region of the HFSs where the column density exceeds 10²² cm⁻², noting that in NGC 1333, the hub had the morphology of a segmented ellipse with a major axis of approximately 0.6 pc. In Mon R2, radial density profiles revealed a distinct turnover in slope at a radius of 0.8 pc, which was used to define the hub radius (Treviño-Morales et al. 2019; Kumar et al. 2022). Morphological analyses of other nearby (<1 kpc) HFSs prompted Kumar et al. (2020) to describe hubs as elliptical structures, with their foci marking the centres of star formation activity. However, the general properties of hubs within the broader context of HFSs remain in need of a more precise definition.

Several observational indicators, such as magnetic field (B-fields) alignment with filaments and flows, suppression of fragmentation (Chen et al. 2019), hour-glass morphologies (Beltrán et al. 2019), and magnetic braking and angular momentum transfer (Hull & Zhang 2019) indicate that B-fields play a significant role in facilitating massive star formation (see also Girart et al. 2009; Pattle et al. 2023). Despite the challenges in constraining their influence, significant progress has been made over the past decade in understanding the role of B-fields in the star formation process (cf. Hennebelle & Inutsuka 2019; Pattle et al. 2023, for reviews). Polarimetric measurements in the far-infrared (FIR) to (sub)-millimetre (sub-mm) wavelengths have become routine with many single-dish and interferometric telescopes, and have enabled detailed mapping of B-fields in star-forming regions. These observations span spatial scales ranging from clouds (~10 pc) (e.g. Planck Collaboration Int. XXXV 2016) to dense cores (~0.1 pc) (e.g. Kwon et al. 2018; Liu et al. 2019). Sub-mm polarization measurements of thermal dust emission from cold, dense molecular filaments reveal specific trends (Palmeirim et al. 2013; Cox et al. 2016; Soler et al. 2016; Planck Collaboration Int. XXXIII 2016; Pillai et al. 2020) in plane-of-the-sky (POS) B-fields:

• (a)

  Parallel to non-self-gravitating, low-density filaments,

• (b)

  Perpendicular to high-density, self-gravitating filaments,

• (c)

  Parallel within very high-density, self-gravitating filaments connected to massive hubs.

These results are based on observations conducted across diverse environments in both low- and high-mass star-forming regions. The relative orientation of the filament axis and B-fields at par-sec scales generally follows the patterns described above but can be altered at core scales (~0.1 pc) due to local phenomena such as outflows, ionizing sources, or the interplay of gravity, turbulence, and B-fields (Pattle et al. 2023).

Recent studies of B-fields in HFSs provide new insights into their structure and dynamics. For instance, Arzoumanian et al. (2021) traced three orders of magnitude in polarization intensity across the 10 pc-long NGC 6334. This work highlights the influence of young massive stars on local B-fields and explores the variation in energy balance – gravity, turbulence, and magnetic forces – throughout the HFSs. In other regions, such as Mon R2 (Hwang et al. 2022) and IRAS 18089-1732 (Sanhueza et al. 2021), striking toroidal magnetic fields have been observed, which reveal an intricate relationship between the morphology of the B-fields and the structure of the HFSs.

Therefore, it is evident that understanding the properties of hubs and magnetic fields in HFSs is key to unraveling the process of massive star formation. HFSs are expected to remain longer than the timescales of local individual star formation events.

This follows the scenario conjectured in the filaments-to-clusters (F2C) paradigm (Kumar et al. 2020), and is consistent with the theoretical expectation that star formation continues and accelerates over many millions of years (Inutsuka et al. 2015; Inutsuka 2017). Since the formation of massive stars generates feedback that can reshape the structure of both the HFSs and the associated magnetic fields, it is equally important to examine how HFSs evolve as star formation progresses. The goal of this study is to explore the global properties of HFSs, particularly the influence of B-fields resulting from the evolutionary stage of the HFSs and its star formation history. The evolution of HFSs is directly influenced by star formation within the hub, especially the development of the OB stellar content (Kumar et al. 2020). To achieve this, we identified three targets within a 2 kpc distance that allowed us to map the full extent of the HFSs while maintaining a spatial resolution of approximately 0.1 pc, necessary to trace the filaments and resolve the hubs. These targets – W3(OH), W3 Main, and S 106 – are well studied in the literature for their distinctive massive star formation activity and evolutionary signposts.

Targets of study: The W3 molecular complex is well known for its rich star formation activity (Rivera-Ingraham et al. 2011), with young massive stars located in W3-Main and the very young, OH and H₂O maser-rich region W3(OH). The entire W3 complex is systematically associated with a distance of approximately 2 kpc. Gaia Data Release 2 (DR2) analysis by Navarete et al. (2019) estimates a mean distance of 2.14 kpc for the complex, with W3-Main and W3(OH) at distances of 2.3 kpc and 2.14 kpc, respectively.

W3(OH) is an ultra-compact HII region (e.g. Wilner et al. 1995), famously associated with the embedded source to its east, an H₂O maser known as the Turner-Welch (TW) object (Turner & Welch 1984), or W3(H₂O). These sources are located within a bright FIR/sub-mm peak, which represents the hub of a network of filaments. Rivera-Ingraham et al. (2013) estimate a luminosity of 2.3×10⁴ L_⊙ and a mass of 1.6×10³ M_⊙ for this source. The estimated diameter of the hub is 0.42 pc, and it contains at least two forming O-type stars as represented by W3(OH) and W3(H₂O). These sources drive multiple outflows (e.g. Zapata et al. 2011) and have been found to drive a rich chemistry in the hot cores (e.g. Giese et al. 2024, and references therein), with various molecular lines tracing the dense gas. Therefore, W3(OH) represents the youngest hub in our study.

W3 Main hosts several hyper- and ultra-compact HII regions with diameters ranging from <1000 au to >20 000 au (Tieftrunk et al. 1997), which is indicative of massive stars in various evolutionary states. This dense region, approximately 1 pc in size, also hosts very luminous infrared sources, the brightest being IRS 5 (L~3×10⁵ L_⊙) (Campbell et al. 1995). IRS 5 is surrounded by a rich cluster of stars (Megeath et al. 2005), including many sub-stellar objects (Ojha et al. 2009). Several H₂O masers are found in correspondence with radio continuum sources around IRS 5, exhibiting radial velocities indicative of powerful molecular outflows (Claussen et al. 1994). Herschel observations (Rivera-Ingraham et al. 2013) reveal the HFSs nature of W3-Main, with a spiral structure (at least five arms). The luminous infrared sources and compact HII regions are located within the dense hub. Early sub-mm observations by Ladd et al. (1993) only traced the elongated dense hub, where they identified the ‘two centres” as SMS1 and SMS2 peaks.

The entire W3 complex hosts clusters of stars that represent a large age spread, as shown by Feigelson & Townsley (2008), who suggest that star formation has occurred in a prolonged fashion, with a time-dependent initial mass function (IMF) in this region.

Within W3-Main, Bik et al. (2014) show that the young low-mass star population is at least 2–3 Myr old, while the massive star population is younger, confirming the well-known sequence of low and high-mass star formation (Kumar et al. 2006; Ojha et al. 2010).

The S 106 molecular cloud is a nearby, 1.09–1.3 kpc (Xu et al. 2013; Zucker et al. 2020) star-forming region, known for its distinct bipolar radiation bubble cavity, wind driven by a late O-type star. This is evident from optical and infrared images from the Hubble Space Telescope (HST) (Bally et al. 1998, 2022). The young massive star, with dust luminosity inferred from sub-mm emission > 10⁴ L_⊙ (Adams et al. 2015), is located in a ‘cold dust bar’ (Mezger et al. 1987), which is clearly a dense hub of an HFS, as evidenced by recent observations of cold dust and molecular gas (Schneider et al. 2007; Adams et al. 2015). The nearly symmetric bipolar radiation bubble implies that the HFS is at a low inclination angle to the plane of the sky, with the hub viewed almost edge-on.

In summary, the three O-star-forming targets – W3(OH), W3 Main, and S 106 – represent an evolutionary sequence, from the youngest (W3(OH)) to the oldest (S 106), respectively. While W3(OH) and S 106 represent HFSs with relatively isolated O stars in early and evolved stages, W3 Main is more complex due to its total mass, multiple OB sources, and rich cluster. This study discusses new sub-mm dust polarization observations of these targets. Centred on the most luminous peak (hub) of the HFSs, the observations cover an area of ~10 pc×10 pc (W3(OH) and W3 Main, d=2.0 kpc) and ~5 pc×5 pc (S 106, d=1.3 kpc), with a spatial resolution of 0.14 pc (W3(OH) and W3 Main) and 0.09 pc (S 106). We examined archival data from the Very Large Array (VLA), and other near- and far-infrared facilities to examine the star formation history in all targets.

This paper is organized as follows: in Sect. 2, we describe the new JCMT observations and the archival data towards the three targets analysed in this work. In Sect. 3, we present the intensity and B-field structures of our three targets. Section 4 presents the observed properties towards the hub and filament regions of the HFSs. In particular, in Sect. 4.1, we propose a method to identify the hub, in Sect. 4.2, we identify different sections in the radial column density profiles of the HFSs centred on the hubs, and in Sect. 4.3 we identify an evolutionary pattern in the filament line mass histograms towards the three HFSs. We then present the observed polarization properties in Sect. 5. We discuss our results and their interpretations in Sect. 6, in particular the relation of hubs to clusters and the identification of the two peaks of star formation within the hubs (Sect. 6.2) and the observational signatures to trace the evolution of magnetized hub-filament systems (Sect. 6.3). Section 7 summarizes and concludes the paper.

2 Observations

2.1 JCMT SCUBA-2/POL-2 observations

We observed W3(OH), W3Main, and S 106 at 850 µm using SCUBA-2/POL-2 (Bastien et al. 2011; Holland et al. 2013; Friberg et al. 2016) installed on the JCMT. The observations were carried out between August and November 2021 under dry (band 2) weather conditions with the standard SCUBA-2/POL-2 CV DAISY polarization mapping mode with a constant scanning speed of 8″ s⁻¹ and a data sampling rate of 2 Hz. Single DAISY fields were observed for W3(OH) and W3 Main, while 3 DAISY fields were mosaicked to map S 106. The spatial distributions of the Stokes I, Q, and U parameters are derived from their time-series measurements using the pol2map data reduction pipeline (including skyloop) that is based on the Star-link routine makemap (Chapin et al. 2013; Currie et al. 2014) and optimized for the SCUBA-2/POL-2 data.

We adopted the flux calibration presented in Mairs et al. (2021). The data are projected onto grid maps with pixel sizes of 4″ at the 14″ HPBW spatial resolution of JCMT at 850 µm (Dempsey et al. 2013; Mairs et al. 2021). The mean uncertainties of the Stokes I maps at 850 µm are δI ~ 2.4 mJy/beam, 2 mJy/beam, and 3.8 mJy/beam, for W3(OH), W3 Main, and S 106, respectively, for the 14″-beam (see also Appendix A.1).

The polarized intensity (PI) and polarization fraction (PF) are calculated using the following relations: $PI=\sqrt{Q^2+U^2}$ and PF = PI/I. PI and PF are further debiased using the modified asymptotic estimator (Plaszczynski et al. 2014). The mean uncertainties of the polarized emission maps are δPI ~ 1.5 mJy/beam for W3(OH), 1.2 mJy/beam for W3 Main, and 1.7 mJy/beam for S 106, for the 14″-beam.

The polarization angle ψ = 0.5 arctan(U, Q) is calculated in the IAU convention, i.e. north to east in the equatorial coordinate system. The POS magnetic field (B_POS) orientation is obtained by adding 90° to the polarization angle ${\chi }_{B_{POS}}=\psi +90°$ (Andersson et al. 2015).

2.2 Column density maps

We estimated the column density ($N_{{\text{H}}_2}$) from the total intensity Stokes I values at 850 µm (I₈₅₀) with the relation $N_{{\text{H}}_2}=I_{850}/\left(B_{850}\left[T\right]{\kappa }_{850}{\mu }_{{\text{H}}_2}{m}_{\text{H}}\right)$, where B₈₅₀ is the Planck function, T = 20 K is the mean dust temperature, κ₈₅₀ = 0.0182 cm²/g is the dust opacity per unit mass of dust + gas at 850 µm (e.g. Ossenkopf & Henning 1994), ${\mu }_{{\text{H}}_2}=2.8$ is the mean molecular weight per hydrogen molecule (e.g. Kauffmann et al. 2008), and m_H is the mass of a hydrogen atom. We adopt T = 20 K, which is the mean temperature of the dense gas (traced by the 850 µm emission) derived from Herschel data towards the W3 (Rivera-Ingraham et al. 2013) and Cygnus (Duarte-Cabral et al. 2013; Schneider et al. 2021) complexes. The column density maps are shown in Fig. A.4.

2.3 Archival data

Near-to-far-infrared archival data of all the three targets were processed to produce overlays comparing stellar and dust continuum distribution. Radio free-free emission in the centimetre wavelengths obtained by the VLA are used to trace the ionizing sources and HII regions.

Ks-band (2.2 µm) images of W3(OH) and W3 Main were created by combining several exposures of the region obtained using WIRCAM on the 3.6 m Canada France Hawaii Telescope (CFHT). The archival data was retrieved by searching the Canadian Astronomy Data Center (CADC). The K-band (2.2 µm) image for the S 106 was obtained from the UKIDSS Galactic Plane Survey (Lucas et al. 2008) which used the 3.8 m United Kingdom Infrared Telescope (UKIRT) and the Wide-Field Camera (WFCAM). The mean full-width half maximum (FWHM) of the point sources are better than 1.0″ in these images.

Spitzer data in the 3.6 µm, 4.5 µm, 5.8 µm, and 8 µm bands were obtained from the Spitzer Heritage Archive. Basic calibrated data (BCD) were retrieved and mosaicked using the MOPEX software with standard mopex scripts for each band. The observations were made in the high-dynamic-range (HDR) mode, with the long exposures having 10.6 seconds.

Herschel data of all the three targets were obtained by searching on the ESA Herschel science archive and processed using the Herschel interactive processing environment (hipe) pipelines specific to SPIRE and PACS data. We combined data from multiple programs to produce a deeper image while effectively recovering saturated pixels in the centre of W3 Main region.

Very Large Array (VLA) data were retrieved from the NRAO data archive by searching for all past observations on the three targets and selecting high-sensitivity and high-spatial-resolutions observations. A large number of observations exist in the archive on these popular targets, observed in several centimetre bands and different interferometer configurations. To match the resolution of the infrared data above and maintain uniformity between targets we retrieved the 4.89 GHz continuum data observed in the B and C configurations. These data have a beam size of ~0.7″.

Fig. 1 Left: total intensity Stokes I map at 850 µm observed with the JCMT SCUBA-2/POL-2 towards W3(OH) in units of mJy beam−1. The contours are I = 30, 100, 500, 1000, and 5000 mJybeam−1. The lowest contour of I = 30 mJybeam−1 is equivalent to I/δI ~ 13. The red lines show the orientation of the POS B-field angle (${\chi }_{B_{\text{POS}}}$) for I/δI > 5, PI/δPI > 1, and PF < 25%. The lengths of the red lines are proportional to PF. A line showing PF = 10% is indicated on the plot. The data are at an angular resolution of 14″ or ~0.14 pc at the 2 kpc distance of the source. Right: contours and red lines are the same as in the left panel. The red lines are plotted with the same length to better show the B-field orientation. The background is a colour image composed using Spitzer 3.6 µm, 4.5 µm, and 8.0 µm images for the blue, green, and red, respectively.

3 Total intensity images and POS B-Field angles

The Stokes I images of W3(OH), W3 Main, and S 106 at 850 µm are displayed in the left panels of Figs. 1, 2, & 3, respectively. The POS B-field is plotted using red lines showing its morphology in each of the three targets. Here, we show the POS B-field lines for I/δI > 5, PI/δPI > 1, and PF < 25%. From Figure A.4 in Appendix A.2, one can see that the data points with PI/δPI > 1 are consistent and coherent with those with PI/δPI > 5. We thus have used PI/δPI > 1 for display purposes while a cut of PI/δPI > 5 is used for the analysis (see Sect. 5). In the right panels of these figures, the normalized POS B-field lines and Stokes I contours are overlaid on Spitzer three colour composites of each target field where the 8.0 µm, 4.5 µm, and 3.6 µm band images are coded red, green, and blue, respectively. These colour images roughly depict the general features of these star forming regions; a) W3(OH) appears as a deeply embedded red object, with a blue companion to its left (east), surrounded by infrared dark regions some of which trace the surrounding filaments. b) W3 Main appears as bluish-white with embedded red sources, the entire region mostly illuminated by active star formation, with infrared dark patches at the outskirts of the region. The exception is a an object to the south east which appears compact and dark, with an embedded red source. We name it foreground HFS because the filaments appear dark due to extinction of background diffuse emission from the W3 Main HII region. c) S 106 prominently reveals the bipolar shaped infrared nebula, which is largely surrounded by the 850 µm emission, the peak coinciding with the IRS 4 massive protostar.

3.1 W3(OH)

W3(OH) (Fig. 1) appears as a centrally condensed peak in Stokes I, with a slight elongation to the east leaning towards the infrared peak that appear bluish in colour (right panel of Fig. 1). The central peak or hub region is surrounded by three prominent filaments, two to the north and one to the south. To the north east, the Stokes I emission reveals a circular shell like structure (dotted circle in Fig. 1 right panel), and adjacent to this another circular shell like structure (marked by black circle) that appears as a bright nebular emission in the Spitzer colour image. Some scattered filaments (marked on the image) can also be seen to the west of the main HFS. The POS B-field lines appear as large (PF) and roughly orthogonal to the filament axis at the outer edges of the X-shaped HFS, smoothly streaming towards the central peak as they become small (PF). The northern shell like features are threaded by an organized POS B-field pattern that smoothly follows the curvature of the shells. Only half of the infrared bright nebular shell is traced by the sub-mm emission but the POS B-field traces the curvature of this shell remarkably well. The smaller filament fragments (marked on the image) display randomly scattered orientations for the B-field.

Fig. 2 Same as Fig. 1 for W3-Main. The data are at an angular resolution of 14″ or ~0.14 pc at the 2 kpc distance of the source.

Fig. 3 Same as Fig. 1 for S 106. The data are at an angular resolution of 14″ or ~0.09 pc at the 1.3 kpc distance of the source.

3.2 W3 Main

W3 Main (Fig. 2) reveals itself as a hub region with two peaks of Stokes I emission, and it is surrounded by lower intensity emission filaments that stretches to much larger area compared to W3(OH). The POS B-field lines appear in groups of relatively scattered large (PF) lines spread all over the field, however appear as well-organized streams of small (PF) lines in the hub region, that seem to focus on the two peaks. W3 Main is known to consist of young stars from multiple episodes of star formation (cf. Sect. 1), and the clumps of scattered large PF values spread all over the observed field might be one of the signatures of multiple generation star formation. A small well condensed peak in Stokes I to the south east traces the foreground HFS. A rough spiral or radial shaped B-field pattern is centred on this object, especially evident from the normalized B-field lines (Fig. 2 right panel). As this object is newly discovered here, we list some of its estimated properties, in Tables 1 & 2, however is omitted from the main analysis given its smaller angular extent and therefore less number of pixels to analyse. (See Sect. 4.2).

To better display the details hidden in the B-field morphology, we zoom-in on the hub region in Fig. 4 (left panel). Here the normalized POS B-field lines are overlaid on the Stokes I image. The thin grey contours represent the VLA 6 cm continuum emission tracing the free-free emission in the compact HII regions in the hub. The UCHII regions A, B, C, D are marked following the notation from Tieftrunk et al. (1997). The youngest most luminous source IRS 5 is also marked. In this figure displaying the central 2 pc region, the B-field structure appears relatively more organized compared to the larger (5 pc) scales shown in Fig. 2. The white arcs aid to visualize this organized pattern. The B-field lines to the top left and bottom right of the image, represent the filaments from the larger HFS, where the B-field appears to form arc shaped patterns pointing towards the hub. This is similar to the commonly observed pattern of the field dragged in by gravity (e.g. Sanhueza et al. 2021). A similar symmetric arc shaped pattern aligns very well along a bipolar massive outflow detected by Li et al. (2019). The B-field pattern in this direction is likely influenced by the powerful massive (4.92 M_⊙) outflow. The driving source of this outflow is not identified, but the coordinates for the outflow centre and its direction are represented by the white star symbol and arrow, respectively. Considering the low angular resolution of the CO outflow data, and absence of any sub-mm peak at this position, it is very likely that the actual driving source of the massive outflow is either IRS 5 or the UCHII region B. The estimated dynamical age of the out-flow is ~1× 10⁵ yr (Li et al. 2019), representative of either of the sources IRS 5 or UCHII region B. Besides the possible influence of the outflow, the B-field lines in the densest region appear to be parallel to the elongations of the total intensity emission around the UCHII regions A and D, which appear to be linked/connected through a horizontal magnetic stream.

The white box inset in the left panel of Fig. 4 represents a foreground region that is further zoomed in the right panel.

Here we use a Spitzer three colour composite background which uncovers a previously unknown mini-HFS with filaments represented by dark lanes converging to the central hub with a cluster of infrared point-sources, of which two are dominant sources. The greenish colour of this cluster implies an excess emission in the 4.5 µm band and therefore indicative of EGO’s (Cyganowski et al. 2008) and outflow activity. This HFS appears to extend up to 1 pc in size and the B-field lines display a rough spiral pattern centred on the cluster of point sources or the hub. In Fig. 2, this foreground HFS appears as a concentrated peak of circular contours to the south-east of the main-hub. In the Herschel study by Rivera-Ingraham et al. (2013), this mini-HFS is identified as the coldest clump W3 Main SE.

3.3 S 106

In S 106 the brightest peak of Stokes I corresponds to the IRS 4 source that is the junction of at least four filaments to the east (marked in Fig. 3). A second peak to the west is connected by three filaments (marked as western filaments). The southern long filament connected to this peak, also appears to trace the walls of the bipolar cavity. Indeed, the central region of this target is significantly shaped by the bipolar radiation bubble/HII region. There are two important features traced by the POS B-field morphology in this region. First, the B-field pattern is aligned along the dense edge-on torus-like structure within which the driving source of the radiation bubble/outflow, namely IRS 4, is located. However, at the outer regions of the filament, at lower Stokes I intensities, the B-field lines appear relatively more scattered, and close to perpendicular to the filament axis at the extreme north-eastern filament. Both on the eastern and western filament groups, the scattered B-fields with large PF values appear to smoothly stream towards the central regions, especially close to the IRS 4 peak as PF decreases. Second, the POS B-field appears to be aligned along the arc-shaped walls of the bipolar cavity, both to the north and the south. The southern dense filament appears to be mostly curved along the edge of the bipolar cavity and even in this filament the B-field lines align along the filament axis. Overall, it suggests that the B-field morphology is significantly re-arranged by the bipolar outflow/radiation bubble.

Fig. 4 Left: blown-up image of the 850 µm Stokes I emission towards the hub of W3 Main (in units of Jy/beam), overlaid by normalized POS B-field lines. Grey contours display VLA 6 cm emission at 0.1″ resolution; the symbols A, B, C, and D mark the known compact HII regions. The youngest most luminous source IRS 5 is also indicated. The white arcs aid visualization of the organized B-field pattern. The white star and arrow indicate the centre of a CO outflow and its direction, respectively (see text for details). Right: W3 foreground mini-HFS, also referred to as W3 Main SE, indicated with the square on the left panel. The Spitzer three colour image is overplotted with Stokes I contours and normalized POS B-field line.

4 Properties of the hub and filament regions towards HFSs

4.1 Defining and identifying the hub centre

The total intensity (Stokes I) maps discussed in the previous section were used to derive the column density maps (Sect. 2.2) and identify filamentary structures in each target. The method to identify filamentary structures is presented in Appendix A.2. Myers (2009) defined hubs as high column density low-aspect ratio objects compared to filaments that are low column density, high-aspect ratio structures. However, given the two central peaks of Stokes I emission, especially in W3 Main and S 106, it is not straightforward to define the boundaries of a hub. To aid the identification of a hub in the same way in all three targets, we have used further deliberations made by the F2C paradigm of Kumar et al. (2020). The F2C paradigm advocates the presence of two embedded peaks within the hub as a result of the evolution of the hub from stage II to stage III as defined by Kumar et al. (2020). Considering the existence of two peaks (for evolved sources) is necessary to first define the centre of the HFSs, and therefore it is the first constraining parameter of the hub.

For this purpose, we examined the literature for all known sign-posts of star formation in the central regions of the targets and display a zoom-in view of the central regions in Fig. 5. Even though W3(OH) and S 106 do not possess two major peaks in the column density maps unlike W3 Main, it turns out that both targets have well-catalogued compact embedded infrared clusters adjacent to the sub-mm peak. In W3(OH), this embedded cluster (Carpenter et al. 2000) reveals itself as the bluish group of stars with a bright blue-white star. W3(OH) itself is surrounded by young deeply embedded infrared stars, but the focus of several detailed studies, mostly using (sub)mm interferometry, has been on the young, luminous pair of objects W3(OH) and W3(H₂O), which are marked in Fig. 5. In W3 Main, the “two peaks” or the “two nodes” of the hub region are well-identified at both sub(mm) and infrared wavelengths. In S 106, the embedded cluster is identified as a loose group of infrared excess young stellar objects (Saito et al. 2009).

We define the central position of the hub as the projected midpoint of the two peaks. Thus the central positions of the hubs are the midpoints between the W3(OH) sub(mm) peak and the bluish-white star of the embedded cluster for W3(OH), the two sub(mm) peaks for W3 Main, and IRS 4 and the centre of the infrared cluster in S 106. The centre of the hubs are marked using a black box-circle symbol in Fig. 5 and listed as the central position in Table 1.

Free-free emission at centimetre wavelengths observed with the Very Large Array (VLA) is an important indicator of the evolutionary state of the ionized region, and used as a standard tool to study compact (ultra and hyper) HII regions. In Fig. 5, we over-plot the best available centimetre continuum data from the VLA data archive to visualize all the compact HII regions in each of the targets. As evident from the figure, they are the most compact in W3(OH), and the most extended in S 106, further justifying our choice of the evolutionary sequence from the youngest to the oldest. Owing to the well-known multiple episodes of star formation in W3 Main, other HII regions from OB stars surround the present episode. These are marked as W 3k and BIRS30 in Fig. 5.

4.2 Measuring the hub size from the HFSs column density radial profile

Having identified the centre of the hub, it is necessary to define a radius or boundary of the hub. To aid this purpose, we produced radial plots of azimuthally averaged column density for each target as shown in Fig. 6. The observed data was fit with between one and three power-laws while allowing fitting ranges to obtain least fitting errors. In all three cases, attempts to fit the data with a single power-law clearly failed. For W3(OH) and W3 Main, three power-laws could be fit, but the third one represented the very outer regions with poor signal-to-noise, which we therefore omitted. Also, in W3 Main and S 106 we deliberately truncated the inner-limit, as the dip within the 0.1–0.2 pc is created by the lower density mid-point region between the two nodes, which is chosen as the hub centre. As evident from the figure, the radial plots are best fit with two distinct power-law slopes, which separate the entire HFSs into two distinct regions. The turn-over point is a natural outcome of minimizing the errors on the two power-law fits. We define this turn-over point in the power-law slopes as the boundary of the hub. With the turnover point as reference, we call the inner and outer regions hub and filaments, respectively. Using the centre we fitted an ellipse to equal column density contours by setting the semi-major axis of the ellipse approximately equal to the scale defined by the turnover point in the radial profiles. Only the minor-axis and PA of the ellipse were fit. The projected major and minor axes and the position angles of the fitted ellipses are listed in Table 1 as hub properties, and are shown as black ellipses in Figs. 5 and 7. Apart from the notions inspired by the F2C paradigm, we also note that the 2D projected geometry in Stokes I emission is approximately elliptical which prompted the choice of an ellipse against other contours or geometries. For the foreground mini-HFS W3 Main SE however, we have not conducted a similar analysis owing to the lack of sufficient pixels encompassed within the JCMT 14″ beam. Only an equal density contour was used to fit an ellipse, the properties of which are listed in Table 1 and shown in Fig. 7. We notice an increasing size of the hub for the early to evolved HFSs evolutionary stage.

The column density radial profiles reveal several important features:

• (a) The central region, not fitted with power-laws, shows a flattening in W3(OH) (≤0.2 pc) and a dip with positive slopes in W3 Main (≤ 0.35 pc) and S 106 (≤ 0.2 pc). This dip is due to the low density valley between the “two peaks” of star formation that are represented by the pairs of sub(mm) peaks and infrared clusters.

• (b) The regions separated by the turnover points (between the two power-laws) are smoothly connected in the early evolutionary stage target W3(OH), and become increasingly distinct as a function of evolutionary time, in W3 Main and S 106. The difference in the slopes of power-laws on either side of the turnover point increases from the early stage W3(OH) target to the oldest region S 106, indicating how the two regions, the hub and filaments, evolve as star formation progresses.

• (c) If a power-law slope of −2 or lower were to be interpreted as suggestive of gravitating gas (Ward-Thompson et al. 1994), and slopes higher than −2 as non-gravitating gas, the hub in W3(OH) is undergoing gravitation collapse, while that in W3 Main and S 106 are not, indeed they may be even relaxing. On the contrary the filament regions in all three targets are gravitationally collapsing to the hub, albeit at different rates. These features are similar to the density profiles discussed for the MonR2 region (Kumar et al. 2022). Even though these observed variations of the power-law slopes contain important physics, we refrain from a detailed interpretation at this point.

Fig. 5 Hub regions of W3(OH), W3-Main, and S106 (clockwise from top left). The colour compositions are the same as the right panels of Figs. 1, 2, and 3. Here the cyan contours trace column densities of 1022, 1023 , and 5×1023 cm−2. The black ellipses represent the hubs as described in Table 1, and their centres are marked with a black box-circle symbol. The blue contours represent VLA 4.89 GHz continuum emission obtained in B/C configurations. This emission is shown with a beam size of 0.7″ for W3(OH), however, it is compact at the scale of 0.08″. In the ‘two node’ system, the OpCl* (Carpenter et al. 2000) represents the more evolved node compared to the younger node composed of W3(H2O) and W3(OH) (e.g. Qin et al. 2016; Wyrowski et al. 1999). In W3-Main, IRS 5 surrounded by compact HII regions (Tieftrunk et al. 1997) A & B represent the younger node, while the evolved node corresponds to the compact regions C & D. In S 106, the younger and older nodes correspond to IRS 4 and the infrared cluster, respectively (Saito et al. 2009).

Fig. 6 Left to right: circularly averaged radial profiles of column densities and associated power-law fits displayed for targets W3(OH), W3 Main, and S106. The shaded region indicates the standard deviation of the values in each annulus. The vertical dotted lines are marked at 0.6, 0.75, and 0.8 pc, respectively for the left to right panels, representing the approximate transition point between power laws. In W3-Main and S 106, the centre of the hub is located in between the ‘two nodes’ (see Sect. 4.1 for details and Fig. 5).

Fig. 7 Maps of filament line mass, Mline (in M⊙ /pc) for W3(OH), W3 Main, and S 106 from left to right. All three panels have the same physical scale of 5.8 pc × 5.8 pc and the same colour scale. The black ellipses indicate the hub sizes as given in Table 1.

4.3 Filament line mass maps of HFSs

The column density maps were used to produce line mass maps. The line mass ($M_{\text{line}}={\mu }_{{\text{H}}_2}{m}_{\text{H}}{N}_{{\text{H}}_2}{W}_{\text{fil}}$) is derived from the $N_{{\text{H}}_2}$ map (Sect. 2.2) and assuming a constant filament width of W_fil = 0.1 pc (Arzoumanian et al. 2011, 2019). In each region, the emission from extended structures has been masked to select only the filamentary structures using curvature maps as described in Appendix A.2. The filament line mass maps are displayed for all the three targets in Fig. 7.

As evidenced in Fig. 7, the depth of our observations are such that the faintest detected structures in all the three targets represent filaments ≳10 M_⊙/pc. Theoretically, filaments with line masses larger than the critical equilibrium value for isothermal cylinders (e.g. Inutsuka & Miyama 1997), $K_{\text{line,crit}}^T=2c_{\text{s}}^2/G~27{\text{M}}_{\odot }/\text{pc}$ (Ostriker 1964), where c_s ~ 0.24 km s⁻¹ is the sound speed at T ~ 20 K, are gravitationally unstable and fragment to form stars. In addition, recent observational studies proposed that the typical mass of the cores/stars formed along a given filament is proportional to the filament line mass (André et al. 2019; Shimajiri et al. 2023; Arzoumanian et al. 2023). Consequently, higher line mass filaments would form preferentially higher mass stars. It can be seen that the central hub regions (represented by the brightest emission and therefore high line mass filaments in Fig. 7) correspond to the centre of (high mass) star formation activity as evidenced by Fig. 5. These findings are in accordance with the formation of massive stars only in hubs postulated by the F2C paradigm (Kumar et al. 2020). Note that the only other region which includes line mass of 100 M_⊙/pc or more in the studied target fields is the W3 Main SE mini-HFS.

The hub region of W3(OH) displays an elliptical condensation enclosing a few hundred M_⊙/pc with a clear peak representing nearly a 1000 M_⊙/pc, which coincides with the well-studied OH and H₂O maser sources (Turner & Welch 1984; Giese et al. 2024), and multiple outflows (Zapata et al. 2011). In W3-Main, the “two nodes” of the hub are distinctly visible, with peaks representing M_line~1000 M_⊙/pc. The left node coincides with the massive young source IRS 5. The right node includes two distinct peaks aligned north-south, the northern peak coinciding with the hyper-compact HII region C (see Fig. 5), while the southern peak may be a younger core/massive stellar object, that can be resolved only with higher angular resolution studies. The hub region of S 106 is relatively more dispersed including mostly M_line~50–100 M_⊙/pc structures, but the peak region corresponding to IRS 4 has two roundish condensations of M_line~300 M_⊙/pc (red colour in Fig. 7) which also drives the bipolar radiation bubble (Fig. 3).

In order to quantify and compare the line mass properties of the hub and filament regions, we used, respectively, all the pixels inside and outside the black ellipses marked in Fig. 7. For each target, the number of pixels enclosed in these regions are listed in Table 2. For each target, Fig. 8 displays the pixel-wise histograms of line masses for the hub and filament regions compared to that of the total field of view. These histograms represent an equivalent filament line mass function (FLMF) (cf. André et al. 2019; Abe et al. 2021, for observational and simulation results, respectively) in these targets¹. Also, in light of the finding that hub is a network of short and very high-density filaments (Kumar et al. 2022), the peaks in Stokes I and therefore column-density actually represent filaments of high line mass that facilitate massive star formation. These filaments within the hub region have been identified using curvature maps as described in Appendix A.2 It can be seen that only in the youngest source W3(OH), the FLMF shows a relatively smooth power-law (at high line mass end) that includes both the filament and hub FLMFs, starting roughly from the 20 M_⊙/pc up to the limit of 2000 M_⊙/pc. However, there is a small departure from the continuous function around 100 M_⊙/pc corresponding to the peak of the hub FLMF. The power-law indices of the filament and hub regions are −1.36 ± 0.26 and −1.2 ± 0.23, respectively. In W3 Main, the histogram shows two distinct peaks, one at M_line~20 M_⊙/pc with a broad distribution, and another narrower at M_line~200 M_⊙/pc that correspond to the filament and hub regions, and are fitted with slopes of −1.04 ± 0.2 and −2.21 ± 0.5, respectively. Similarly in S 106, the histogram has two peaks each corresponding to the peak of filament and hub histograms with power-law indices −0.86 ± 0.17 and −1.41 ± 0.28, respectively. Comparing the FLMFs between these regions it appears that at the earliest evolutionary stage represented by W3(OH), a Salpeter-like FLMF is witnessed in the hub which gradually weakens as star formation progresses represented by the −0.86 power law for S 106. In contrast, the FLMF in the filament region appears to be similar in the earliest evolutionary stage and the most evolved stage as represented by W3(OH) and S 106. However, the FLMF has a steeper power-law in W3 Main when the star formation activity is highest, as evidenced by the numerous compact HII regions and the young massive stars in the hub.

In summary, similarly to the hub size and the radial column density profile, we find distinct signatures of the filament line mass distribution for each of the HFSs reflecting their different evolutionary stages. Next, we examine how polarization properties change as a function of the evolution of the HFSs and their star formation activity.

Fig. 8 Line mass histograms for the full image (black dashed line), filament (blue), and hub (red) regions for W3(OH), W3-Main, and S 106 from left to right. The fitted power laws are shown by blue (filament) and red (hub) lines, and the fitted slopes are indicated in each panel.

5 Polarization properties

In this work, because we have assembled targets at different evolutionary stages with a variety of star formation sign-posts from stars to outflows, we can compare the polarization properties as a function of evolutionary stage to better understand the impact of the star formation activity and feedback on the polarization properties of the dense gas.

5.1 Distributions of PF and PI in hubs and filaments

Having defined and identified the hub and filament regions of the HFSs in each target, we now examine the behavior of the observed polarization properties in these two regions and compare with that of the region as a whole for each HFS and among the targets. In Figs. 9 and 10 respectively, we plot the histograms of polarized intensity (PI) and polarization fraction (PF). In producing these histograms we used maps with 4″ size pixels, and have applied the combined criteria that the plotted data points have PI/δPI > 5 and error in PF, δPF < 3%. We have, however, examined the histograms with a more relaxed constraint of error in PF with δPF < 3%, which only adds a few pixels at high PF values arising at the edges of filaments. On the contrary the plotted data points are robust and represent the crest regions of filaments and within the hub.

It can be seen from Fig. 9, that the highest and lowest values of polarized intensity are respectively found in the youngest and oldest targets namely W3(OH) and S 106. Even though W3 Main is the most massive HFS of all the three targets, the hub region of W3(OH) has larger number of pixels at PI >80 mJy/14″ beam. The PI histograms of S 106, due to its late stage of evolution, stand out from the other two targets, owing to a spatially resolved bipolar radiation bubble.

For W3(OH) and W3-Main, the PF histograms show a broader distribution for the filaments with larger values (up to 14%) reached for W3(OH) compared to the narrower distributions for the hub region peaking at values around 1% (Fig. 10). The PF histograms of S 106 in the hub region shows a broader distribution peaking at 3% with values reaching up to 8–10%.

Fig. 9 Polarized intensity (PI) histograms for the total (black), hub (red), and filament (blue) regions for W3(OH), W3-Main, and S 106 from left to right. The vertical axes indicate the number of pixels per bin.

Fig. 10 Polarization fraction (PF) histograms for the total (black), hub (red), and filament (blue) regions for W3(OH), W3-Main, and S 106 from left to right. The vertical axes indicate the number of pixels per bin.

5.2 Spatial variation of PF and PI

The spatial variation of PF and PI in all three targets were examined and compared with the location of dominant massive young stellar sources and the sub-mm peaks in each target. Because PF diminishes inside the highest density regions, we computed 1/PF and plotted it, along with PI, on the Herschel SPIRE 250 µm images as shown in Fig. 11. The most luminous (massive) young stellar objects and known massive outflows are marked.

In W3(OH) the sources OH and W3(H₂O) are embedded in a disc-like structure traced by the 1/PF contours (black) while PI (white contours) reveals a distinct bipolar shaped structure centred on this disc. Using SMA observations, Qin et al. (2016) uncover two outflows originating from W3(H₂O) which is roughly aligned north-south but tilted more towards north-east compared to the cyan vector shown in Fig. 11. Understandably, the disc-shaped structure traced by the black contours hosts several massive young stellar sources at very early evolutionary stage, including the OH source and the three cores associated with the H₂Omaser Taylor–elch object. The bipolar pattern from our 14″ resolution maps indicate the net result of the multiple outflows in the region, all of which are effectively in the north-south direction (see Fig. 4 of Qin et al. 2016) The disc-like structure is well centred on the peak of 250 µm emission. The OH maser source is associated with a hyper-compact HII region as uncovered by the VLA observations (Dzib et al. 2014), it is a relatively evolved source with no outflow detection from Qin et al. (2016).

The PI and I/PF contours are more complex in W3 Main owing to a number of massive ionizing and young stellar sources crowded in the two nodes of the hub. In the left node, the I/PF contours reveal an elongated structure (marked with a blue vector) on either side of IRS 5 (coinciding with the 250 µm peak) perpendicular to which lies a massive outflow mapped in CO by Li et al. (2019) (see also Sect. 3.2). The PI contours are not striking in tracing this bipolar outflow (cyan vector), however, two lobes can be visualized that roughly match with the axis of the outflow. In the right node of the hub, the PI and I/PF contours reveal complementary structures; PI contours trace a clear bipolar structure, at the centre of which 1/PF forms a condensed peak, coinciding with the VLA source C. The 1/PF contours can also be seen to trace the walls of the northern bipolar lobe. The 250 µm emission here is split between the condensation traced by 1/PF and also with the northern lobe of the bipolar pattern traced by PI. These two peaks of 250 µm may correspond to the two small peaks observed in the Stokes I and consequently the line mass maps (see Sect. 4.2). SMA and IRAM 30m observations of the full hub-region of W3 Main in both continuum and spectral lines (Wang et al. 2012) reveal that the region contain several massive sources and multiple outflows that render this region complex. Our 14″ maps here are tracing the major out-flow cavity walls and associated cores where the driving sources are located.

In S 106, the peak traced by 1/PF coincides well with the source IRS 4, and while both 1/PF and PI trace similar structure, the peak traced by PI is shifted to the left of IRS 4. The region to the north and south, covered by the well-studied S 106 bipolar radiation bubble/outflow lacks significant 850 µm emission thus no comparison can be made.

Fig. 11 Clockwise from top left: spatial correlation between total intensity, PF, and PI in W3(OH), W3 Main, and S 106 displayed over Herschel SPIRE 250 µm image. Black contours represent 1/PF with levels 0.25, 0.5, 1, 1.5, and 2 (for W3(OH) and W3 Main, top) and 0.1, 0.3, 0.6, and 0.9 (for S 106, bottom). White contours display PI with levels from 20 to 90 mJy/beam in intervals of 10 mJy/beam for W3(OH) and W3 Main (top), and levels 5, 10, 20, 30, 40, 50, and 60 mJy/beam for S 106 (bottom). The cyan vectors in the top two panels mark bipolar patterns in PI that correspond to outflows (see Sect. 5.2). The blue vector in W3 Main identifies the low PF elongated structure centred on IRS 5. This is similar to the elongated disc structure in W3(OH) traced by 1/PF contours.

5.3 POS B-field angle (${\chi }_{B_{\text{POS}}}$) versus filament orientation (θ_fil)

In order to examine how ${\chi }_{B_{\text{POS}}}$ behaves with respect to the different filamentary structures over the entire region, we computed the difference between ${\chi }_{B_{\text{POS}}}$ and the filament orientation θ_fil derived from our filament map (cf. Appendix A.2). These maps are shown in Fig. 12 for all the three targets. It can be seen that the densest filaments, especially those close to the central regions of the HFSs where the hub is located, appear in dark blue in Fig. 12 representing a small difference between ${\chi }_{B_{\text{POS}}}$ and θ_fil. This suggests that the ${\chi }_{B_{\text{POS}}}$ is more aligned with the axis of the filaments towards the hubs. In contrast, the large difference between θ_fil and ${\chi }_{B_{\text{POS}}}$ (represented by pink and redder colours) is found in the lower density areas of the filaments. In W3(OH) the dark bluer colours are concentrated towards the hub with little or no redder colours. In W3 Main the small angle differences are prominent both in the central hub and also in the foreground mini-HFS SE source. However, there is also a significant mixture of larger angle differences in the W3 Main hub, unlike W3(OH). This is likely due to the presence of multiple compact HII regions and young bipolar outflows in the central hub of W3 Main that is in the process of altering the density structure of the hub region. This mixed features of both large and small angle differences in the central hub region is most prominent in S 106 where a well-resolved bipolar outflow/radiation bubble driven from the IRS 4 has significantly reshaped the surrounding material.

To further explore these maps, we also plot the histograms of the difference angle, ${\theta }_{\text{diff}}=\left|{\chi }_{B_{\text{POS}}}-{\theta }_{\text{fil}}\right|$, for pixels located in the hub, filament, and the full field as shown in Fig. 13. In producing these histograms only pixels with errors on ${\chi }_{B_{\text{POS}}}$ smaller than 20° are selected. The choice of ${\chi }_{B_{\text{POS}}}$ error was made by comparing the distribution of PI/δPI > 1 and > 5 with the corresponding ${\chi }_{B_{\text{POS}}}$ error which showed that the major concentration of data points was within an error of <20° with no particular correlation to the PI/δPI cuts. In W3(OH), the histogram of θ_diff for the entire region shows a significant peak around 10° indicating that the POS B-field is mostly parallel to the filaments, both inside and outside the hub, albeit with a larger number of pixels with perpendicular or random orientation along the filaments outside the hub. The histogram of the full region of W3-Main also peaks at ~ 10–15° however with a broad power-law like distribution up to ~80–90° indicating larger mis-alignments and perpendicularity between the POS B-field and the filaments. This trend (peak towards parallel orientations and power-law distribution of large θ_diff values) is observed both for the hub and the filament regions. Finally, in the most evolved S 106 HFS, the results are different with the presence of two peaks, the first largest peak at angles <45^° and the second smaller peak for >45^°. This distribution is mostly dominated by the distribution in the hub, while the angle difference for the filaments shows a random distribution. These histograms show noticeable differences between our three targets possibly tracing the evolution stages of the HFSs, as discussed further below.

Fig. 12 Maps of the difference between the position angle of the filaments θfil and the ${\chi }_{B_{\text{POS}}}$ angle for W3(OH), W3-Main, and S 106 from left to right. The black ellipses indicate the hub sizes as given in Table 1 and shown in Fig. 7.

Fig. 13 Histograms of the difference between the filament orientation and ${\chi }_{B_{\text{POS}}}$ for the total (black), hub (red), and filament (blue) regions, for W3(OH), W3-Main, and S 106 from left to right. The vertical axes indicate the pixel number. Here, pixels with the error on ${\chi }_{B_{\text{POS}}}$ smaller than 20 degrees are chosen.

6 Discussion

In this work, we examined the properties of three star cluster-forming HFSs in different evolutionary stages in an attempt to understand their evolution. We have combined observations of star-forming dense gas and young stellar population to obtain our results. The total intensity images at 850 µm trace the dense gas in the HFSs, identifying the filament and hub regions in each target, while allowing us to quantify column densities and line masses. The changes in polarized intensity and polarization fraction with the progress of star formation and evolution of the HFSs have been examined. The HII region emission from VLA and the information on outflows from the literature have served to assess the influence of feedback during the evolution.

6.1 Polarization and magnetic field properties of hub-filament systems

We analysed JCMT POL-2 data of polarized dust emission at 850 µm for three selected HFSs in consecutive evolutionary stages, from early to evolved, to investigate the global evolution of the magnetic field structure. The histograms of relative orientation, ${\theta }_{\text{diff}}=\left|{\chi }_{B_{\text{POS}}}-{\theta }_{\text{fil}}\right|$ between the B-field angle (${\chi }_{B_{\text{POS}}}$) and the filament orientation (θ_fil) shows a variation; 1) between the hub and the outer filamentary parts of the HFSs, and also 2) as a function of the evolutionary stage of the HFSs (Fig. 13). At early stages (for W3(OH)), ${\chi }_{B_{\text{POS}}}$ and θ_fil are mostly aligned with a small dispersion of the θ_diff histogram. At the intermediate stage (for W3-Main), the θ_diff distribution broadens with a peak at small angles of parallel orientation and a power-law distribution up to θ_diff~90º for perpendicular relative orientation. At the most evolved stage of S 106, the θ_diff histogram towards the hub shows two peaks, with a stronger peak for the parallel configuration and a weaker peak for the perpendicular configuration.

The current theoretical models propose the formation of high line mass filaments through accumulation of matter through flows along the B-field lines resulting in filaments perpendicular to the B-field (e.g. Inoue et al. 2018; Pineda et al. 2023; Abe et al. 2025). This theoretical prediction is compatible with the observations of high line mass filaments not connected to hubs (e.g. Palmeirim et al. 2013; Cox et al. 2016; Ching et al. 2022; Pattle et al. 2023, to refer to a few examples). In HFSs however, observations show that the relative orientation between the B-field and the filaments changes from mostly perpendicular to mostly parallel as the density increases towards the hub region (e.g. Wang et al. 2020; Pillai et al. 2020; Arzoumanian et al. 2021; Wang et al. 2024). The reorientation of the B-field from perpendicular to parallel along the high-line mass filaments connected to massive hubs has been interpreted as resulting from dynamical flows along the filaments carrying the matter and the frozen-in B-field towards the hub (as suggested by MHD simulations, e.g. Gómez et al. 2018; Suin et al. 2025). Consequently, we interpret the parallel configuration of filaments and B-field in our three targets as the influence of the gravity of the hub onto the surrounding filaments and B-field. In addition, the parallel relative orientation between the filaments and the B-field could result from projection effects as discussed for example in Doi et al. (2020).

Previous observational studies of star forming clouds have shown that the POS B-fields are reorganized tangentially along the swept-up material in the shells formed by the expanding HII regions (Tang et al. 2009; Planck Collaboration Int. XXXIV 2016; Arzoumanian et al. 2021; Cortés et al. 2021; Tahani et al. 2023). We observe a similar effect in the northeastern region of W3(OH) (Fig. 1), though it is not directly associated with a known HII region in this case. In contrast, W3 Main exhibits a B-field geometry modulation that can be partially attributed to the compact HII regions present within its hub (Fig. 2). Moreover, an hourglass-shaped pattern embedded within the hub (Fig. 4 left) aligns with the most massive molecular outflow identified in the region (Li et al. 2019).

We also find that the POS B-fields can also be reshaped by massive outflows such as in Lyo et al. (2021) and their associated radiation bubbles, as evident in S 106 and the central region of W3 Main (Fig. 4 left). In S 106 (Fig. 3), the B-field lines trace the cavity walls of the well-documented massive outflow and radiation bubble driven by IRS 4. This reordering of the B-field along the outflow cavity walls may likely explain the high PF values observed in the hub of S 106 (cf. Fig. 10).The increase of PF in this context may result from a combination of the increase of grain alignment efficiency from the anisotropic radiation from the massive central star (according to the RAT theory Lazarian & Hoang 2007; Hoang et al. 2021) and the reduced depolarization due to the reordering of the B-field structure. Outflows reshaping the B-fields have been previously seen in the study of numerous other targets (Hull et al. 2017; Fernández-López et al. 2021; Pattle et al. 2022; Karoly et al. 2023; Encalada et al. 2024).

Across all observed HFSs targets, on average lower PF values are consistently found towards the hubs as opposed to the lower column density emission outside the hubs. Such a decrease of PF with increase in column density is compatible with previous observations (cf. Pattle et al. 2023, for a comprehensive review) interpreted as resulting from the combined effects of loss of grain alignment, variation of grain properties, change of the 3D B-field configuration, and depolarization due to averaging along the line-of-sight and within the beam.

A particularly interesting feature in our data is the correlation of low PF regions, tracing the dense cores and/or discs/toroids, with the geometric centres of bipolar outflow-shaped patterns traced by PI, most prominently seen in W3(OH) and W3 Main near ionizing sources (Fig. 11). The edge-on disc-like morphology of the 1/PF contours in W3(OH) likely represents magnetized, self-gravitating toroidal structures, analogous to those observed in the G31.41 HFS (Beltrán et al. 2019, 2024). The G31.41 HFS remains one of the most thoroughly studied example of massive star-forming HFSs, with detailed observations across multiple scales, angular resolutions, wavelengths, and tracers. Beltrán et al. (2024) report that PF values are approximately 2% around four spatially resolved individual massive sources and increase up to 10% in the broader hub. The observations of G31.41 are modeled as a magnetized toroidal structure hosting massive sources, with an hourglass-shaped B-field pattern emerging along the polar direction. The PF and PI patterns observed in Fig. 11 align with the scenarios depicted by Beltrán et al. (2024). These findings underscore the complex interplay between gravity, magnetic fields, massive outflows, and the radiation environment within HFSs, offering new perspectives on the role of B-fields in massive star formation and cluster evolution.

6.2 Star formation in hubs

6.2.1 Hub definition and its relation to cluster formation

In this study, for the first time, we define the ‘hub’ based on column density radial profiles of multiple targets. We do so by demonstrating that the column density distributions in our three targets can be approximated as elliptical structures, with the radial density profiles showing distinct transitions at the boundaries of these ellipses (cf. Fig. 6). Similar definition has been previously used for Mon R2 by Treviño-Morales et al. (2019) and Kumar et al. (2022) resulting in a radius of 0.8–1.2 pc (see Table 1). The systematic turnover of column-density slopes at the hub boundaries in different targets suggest a deeper implication that may shed light into the physics of star formation.

One intriguing feature of young stellar clusters (YSCs) in the Milky Way is that the cluster radii are highly constricted (median ${\text{R}}_c={0.58}_{-0.20}^{+0.32}\text{pc}$ for d < 2 kpc, Kharchenko et al. (2013)) which is larger for massive (> 10⁴ M_⊙) YSCs (${\text{R}}_c={1.5}_{-0.42}^{+2.35}\text{\hspace{0.17em}pc}$) (Portegies Zwart et al. 2010). In addition, the cluster radii have a weak dependence (Krumholz et al. 2019) on mass (R_c ∝ M^1/3, up to 10⁵ M_⊙). We note that the observed cluster sizes are in close agreement with the hub sizes (~1 pc), derived both in this work and other regions (example Mon R2 Treviño-Morales et al. 2019; Kumar et al. 2022). This correlation is intriguing, given that star formation takes place in the filaments of the HFSs which are several times larger, ranging from 3–10 pc (Kumar et al. 2020).

6.2.2 Double-nodes or two peaks of star formation in hubs

In Kumar et al. (2020), we first highlighted that nearby star-forming regions such as Orion and NGC 2264 exhibit two distinct centres of star formation within their central hub regions, each demonstrating a clear evolutionary difference: one centre being younger than the other. These nodes within the hub represent the densest and most massive gas concentrations, displaying characteristics indicative of ongoing star cluster formation. In this study, we extend the initial identification of the two star formation peaks by Kumar et al. (2020) through a clear detection of the corresponding nodes in all three targets analysed. We have shown that the two nodes are consistently associated with discernible concentrations of young stars at distinctly different evolutionary stages (Sect. 4.1). Remarkably, these features remain evident across all three targets studied here, regardless of their evolutionary state.

In Table 3, we compile a list of well-studied star-forming regions that we propose as examples of these “two star formation nodes.” The separation between these nodes is generally less than 1 pc, with exceptions such as NGC 2264 (Peretto et al. 2006) and the Aquila Rift HFSs (Bontemps et al. 2010). While there is some debate and discrepancy in distance estimates for regions like W40 and Serpens South (see Bontemps et al. 2010, for discussion), it is noteworthy that both sources are centrally located in the large-scale Herschel maps presented by these authors. The smallest node separation is observed between Orion’s BN/KL and the Trapezium cluster, at 0.14 pc.

The double node of star formation has not yet been systematically explored, and its prevalence as a common feature of HFSs remains unclear. However, in the well-studied regions listed in Table 3, the presence of both young and older nodes appears systematic. These clumps are distinguishable not only through clear indicators of youth and evolved states but also by notable temperature differences, even when located in close proximity. For instance: NGC 6334 I(N) is colder than NGC 6334 I (Arzoumanian et al. 2022, see Fig. 1). In NGC 1333, the SSV 13 cluster is embedded within a colder clump compared to the SVS3/BD+30°549 region (Hacar et al. 2017, see Fig. 2). Despite these observations, the physical origin of such well-defined dual centres of star formation activity remains largely unexplored and poorly understood.

We hypothesize that the two observed star-forming nodes may represent different stages in a sequential star formation process within the theoretical expectation that star formation continue (and accelerate) over many million years (Inutsuka et al. 2015; Inutsuka 2017). Given the timescale over which stellar feedback disperses gas, it is plausible that only two such nodes – one relatively younger than the other – are visible within a single dense hub at any given time.

We further suggest that earlier star formation episodes, now traced by young star clusters, are located near – but no longer physically associated with – the currently active star-forming hub, having lost their connection to dense gas. A notable example is the Orion Trapezium cluster, located near BN/KL, which is no longer associated with dense gas, likely due to its dispersal by stellar feedback (see Fig. 1 of Hacar et al. 2018). Yet, it represents a major star formation node from the previous episode. Conversely, the Orion South dense clump, situated adjacent to BN/KL, could represent the next or emerging star formation node (also found in Fig. 1 of Hacar et al. 2018).

Extending this idea, we propose that dense hubs in very early evolutionary stages may host only a single, active star-forming node alongside a second, pre-stellar one. A potential example is the Mon R2 hub: it contains a cluster of five IRS sources forming one node, and a cold, dense pre-stellar clump or filament located ~0.3 pc southeast of the IRS cluster, observed in sub-millimetre emission (the elongated highest density (white colour) feature in Fig. F.1 of Didelon et al. 2015), forming the second.

Based on these observations, we propose that the ‘double-node’ configuration of star formation represents a continuous evolutionary sequence. The spatial configuration of dense gas (the hub) and associated star clusters changes over time. At a given epoch, one node appears young while the other is older. In the previous generation, the younger node may have been pre-stellar; over time, as stellar feedback clears out dense gas or accretion reorganizes the hub’s structure (e.g. through shock compression or filament collisions), the older node evolves into an exposed star cluster located outside the current dense hub.

Future theoretical modelling and more systematic observational studies of star-forming nodes - both within and around dense hubs – will be crucial to better understand the nature and implications of this proposed ‘double-node’ configuration in HFSs.

6.3 Evolution of magnetized hub-filament systems

Hub-filament systems are now well established as progenitors of star cluster formation. The initial mass function (IMF) reliably characterizes the mass spectrum of stars within these clusters, where low-mass stars dominate not only in number but also in mass, because the IMF power-law exponent is smaller than −2. The relatively few massive stars are known to form inside hubs, possibly ubiquitously (Kumar et al. 2020). The goal of this study has been to identify indicators of the evolution of the HFSs and in particular hubs, as star formation proceeds.

6.3.1 Line mass evolution

As elaborated in Sect. 4.3, the central intensity peaks that constitute the hubs naturally enclose the highest line mass filaments and the overall slopes of the FLMF (represented by line-mass histograms in Fig. 8) display discernible variations as the HFSs evolves. In particular a smooth Salpeter-like slope enclosing both the hub and non-hub regions in W3(OH) breaks out into two different slopes in the intermediate and evolved targets W3 Main and S 106 respectively. This likely indicates that star formation proceeds significantly in the hub-region turning dense gas into stars. As the highest line mass filaments in the hub-region result in massive star formation and fragmentation producing other low- intermediate-mass stars, one can expect the FLMF slope turning steeper as evidenced in W3 Main and S 106 histograms. In contrast, the external filament regions do not display significant alterations in their line-mass distributions. At very evolved stages nearing the end of star formation and destruction of the HFSs, even the filament regions may change their mass distributions.

6.3.2 Hub sizes

Next, as evident from Table. 1 (see Sect. 4.2), there is a noticeable increase in the size of the hub from the early to evolved HFSs evolutionary stage, in-spite of the significant differences in their total mass and luminosities. Even though this effect cannot be generalized given the limited sample size, the trend is important because the “hub size” effectively reflects the final agglomeration of the dense gas that is responsible for the formation of the young stellar cluster. As the HFSs evolve, the process of accumulating dense gas from the outermost parts of the HFSs to the central hub is nearly complete, effectively representing the maximum size of the hub. At the same time, the feedback effects from massive star formation in the hub would have its impact on the lowest density matter, especially at the outer regions of the HFSs. This is because feedback propagates most easily in the low density inter-filamentary cavities, eroding the low density material. The highest density structures in the hub are the last structures to wane at the end of star formation.

One interesting feature from Fig. 13 (see Sect. 6.1) is the distinct peaks in B-field orientation angles with respect to the filament in the most evolved stage. In the hub, the orientation is mostly parallel. Considering again the limited sample, and projection effects, we cautiously suggest that this observation may be suggesting that the dragging of B-fields along the highest density filaments is mostly complete. While this process of gravity dragging the B-fields begins at the earliest stages of the HFSs, finds it maximum effect at the most evolved stage.

6.3.3 Evolution of PF and PI

The distributions of PF for the hub region is found to broaden with the evolutionary state of the target (see Sect. 5.1 and Fig. 10). In comparison the outer filament regions maintains similar widths for the PF distribution in the case of W3(OH) and W3 Main but display a large scatter in S 106. Given that low PF values at high column density is a product of several possible effects such as loss of grain alignment, variation of grain properties, and depolarization due to line-of-the-sight averaging, we can only speculate possible reasons. In the early stages such as W3(OH), the presence of nearly edge-on toroidal structures (Fig. 11) can result in depolarization effect from line-of-sight averaging. As star formation proceeds, and massive young stars emerge, shedding much uv radiation on the surroundings, increased PF in S 106 can be the result of re-alignment of dust grains due to radiation. Also, the highest and lowest PI values in the hub are found in the youngest and oldest targets (Sect. 5.1 and Fig. 9). In essence, at early evolutionary stages, the high density hub structure is more ordered resulting in less depolarization in W3(OH) compared to the more evolved W3 Main. At later stages, the emergence of young stars and associated radiation and waning dense structures can result in higher PF values with broader distributions and decreasing PI values, as evident from S 106.

6.3.4 B-fields in HFSs and their evolution

Massive stars preferentially form in hubs Kumar et al. (2020). While magnetic fields initially aid the formation of massive stars, eventually they get reshaped from their stellar feedback. There is mounting observational evidence (Beltrán et al. 2019) and theoretical considerations (see, for example, Krumholz 2015) pointing to the critical role of magnetic fields in the assembly of massive stars. Inside the hubs, magnetic fields become especially significant by creating the ideal conditions for regulated mass accumulation and providing support against gravitational collapse at the highest densities. An increasing number of multi-scale observations of magnetic fields (Liu et al. 2023; Beltrán et al. 2024; Koch et al. 2022) have revealed their structure from HFSs to core scales. The few-parsec scale observations encompassing the three HFSs presented in this study, combined with detailed B-field analyses of G31.41 HFS (Beltrán et al. 2024; Law et al. 2025) and NGC 6334 (Liu et al. 2023) at sub-parsec to circumstellar scales, offer some of the most comprehensive insights to date on the role of magnetic fields in massive star formation.

From the roughly orthogonal orientation in the outer filament regions, the B-fields reorient to become parallel at the highest densities inside the hub. This shift in orientation underscores the importance of magnetic-pressure support – even more than turbulent support – in massive star-forming hubs. An increasing number of observational studies corroborate the idea that magnetic support plays a more dominant role than turbulence in these environments (Beltrán et al. 2019, 2024; Rawat et al. 2024). The unique configuration of HFSs is critical in this context, as they naturally assemble the highest density network of filaments within their hubs, a prerequisite for high-mass star formation (see e.g. Kumar et al. 2022). As higher densities are assembled within the hubs of HFSs – via mechanisms such as filament coalescence - B-fields undergo modulations in both strength and morphology. These changes appear to provide the magnetic support necessary (Girart et al. 2009; Krumholz 2015) to restrict fragmentation (Chen et al. 2019) and sustain continued mass accretion onto massive stellar seeds. Observations of magnetically supported toroidal structures in G31.41 HFS (Beltrán et al. 2019) and disc-like magnetic structures around W3(OH) offer compelling evidence supporting this proposition.

Once HFSs begins forming young stellar clusters, especially with massive stars, feedback mechanisms - such as outflows, radiation bubbles, and HII regions – reshape the surrounding B-fields (Tahani et al. 2023; Khan et al. 2024). For instance, Liu et al. (2023) find that in NGC 6334, at the highest densities, B-fields exhibit no preferred direction or show slightly perpendicular orientations relative to the dense gas. This is largely attributed to feedback-driven outflows. Our observations further suggest that this magnetic field reshaping predominantly occurs in the lower-density media surrounding the dense filaments. Within the densest structures themselves, B-fields remain largely unaltered, continuing to provide critical support for mass accumulation and regulating the star formation process. These results reinforce the notion that magnetic fields are integral to both the initiation and evolution of massive star formation within hub-filament systems, offering new insights into the multi-scale dynamics of these complex environments.

7 Summary and conclusions

In this study, we investigated the characteristics of three star-cluster-forming hub-filament systems – W3(OH), W3 Main, and S 106 – that represent early, intermediate, and advanced evolutionary stages, respectively. We presented new observations of dust-polarized emission with JCMT/SCUBA-2/POL-2 and reanalysed archival near-infrared and radio data from the VLA, WIRCAM, Spitzer, and Herschel. Our main findings are summarized below:

• Column density and line mass maps derived from the 850 µm dust emission reveal the hub-filament nature of these regions, with central hubs showing the highest column densities and line masses ranging from several hundred to a thousand M_⊙ pc⁻¹. The plane-of-sky magnetic field angles (${\chi }_{B_{\text{POS}}}$) were compared to Stokes I maps, showing that B-flelds are predominantly perpendicular or randomly orientated with respect to lower-density filaments in outer regions, which suggests the possibility that B-fields are mostly perpendicular to the filaments in 3D (Doi et al. 2020). However, they become aligned with higher-density structures within the hubs;

• Archival VLA and near-infrared data, together with column density maps, reveal two primary star formation nodes within the hubs. These nodes exhibit clear evolutionary differences, such as young stellar clusters visible in near-infrared adjacent to deeply embedded massive YSOs or compact HII regions;

• Ellipses were fitted enclosing the two nodes of star formation. Azimuthally averaged radial column density profiles, centred at the midpoint between these two nodes, were used to define the hub. These profiles exhibit a two-power-law structure, transitioning at r~0.6–0.8 pc, with power-law indices > − 2 inside and < − 2 outside the hub boundary (Fig. 6). As the HFSs evolve, we notice an increase in the hub size, and a steeper FLMF slope inside the hub;

• Histograms of polarization fraction (PF) and polarized intensity (PI) for both hub and filament regions illustrate an evolutionary trend from the youngest system, W3(OH), to the most evolved, S 106, indicating the influence of star formation feedback. Contour maps of 1/PF and PI reveal disc-and bipolar outflow-shaped patterns centred on the most luminous sources in each region (Fig. 11). We observed that massive bipolar outflows and radiation bubbles significantly reshape the magnetic fields, which align along the outflow cavity walls. This feature is consistent across all three targets. The sources driving these outflows correspond to regions with the lowest PF values, likely due to large grain destruction or increased grain alignment from enhanced anisotropic radiation of young stars;

• Histograms comparing the difference between (${\chi }_{B_{\text{POS}}}$) and filament orientation angles were generated for both hub and filament regions. The results align with observed morphologies, indicating that B-fields are generally parallel to dense hub structures and perpendicular in lower-density filament areas. In the evolved system S 106, these angle differences display both parallel and perpendicular alignments within the hub, which suggest modifications driven by outflow activity;

• In W3 Main, we identified a previously unreported spiral-shaped mini-HFS in the foreground, which had been recognized as the W3 Main SE clump in prior Herschel studies;

• Alongside our current targets, we compiled a list of nearby star-forming regions where a ‘double-node’ pattern of star formation is evident, consistently showing evolutionary differences between adjacent nodes. This pattern may be a common feature in many HFSs and could reflect the physical conditions governing their formation and evolution.

This is the first study to examine HFSs properties as a function of evolutionary stage using only three targets. To better understand the change in the properties of HFSs as they evolve in time and as star formation proceeds, a statistical analysis of a large sample is required. Such a study should examine the distribution of hub sizes, line-mass functions, and magnetic field orientations as a function of evolution. High angular resolution data are needed to resolve the structure of the hubs (e.g. ALMA, JWST) and molecular line observations are needed to describe the kinematics. Quantitative comparisons of results from MHD numerical simulations at different epochs of the simulation will also be important to understand the physical processes responsible in the formation and evolution of star-cluster-forming hub-filament systems.

Acknowledgements

We thank the anonymous referee for a very thoughtful, encouraging and hawk-eyed review that greatly enhanced the clarity, bibliographic completeness and quality of this manuscript. This research is partially supported by Grants-in-Aid for Scientific Researches from the Japan Society for Promotion of Science (KAKENHI 18H05437, 19H01938 and 21H00045). MSNK was supported as a visiting scientist at the Theoretical Astrophysics Laboratory, Nagoya University through the aid of KAKENHI 18H05437. RSF was supported by the Visiting Scholars Program provided by the NAOJ Research Coordination Committee, NINS (NAOJ-RCC-23DS-050). NKB acknowledges the support of the China Postdoctoral Science Foundation through grant No. 2025M773187. Data analysis was in part carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the National Astronomical Research Institute of Thailand; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations in the United Kingdom, Canada, and Ireland. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community.

Appendix A Stokes I maps and selection of the filamentary emission

Appendix A.1 Stokes I maps

The left panels of Figs. A.1, A.2, and A.3 show the full Stokes I map at 850 µm observed with the JCMT SCUBA-2/POL-2 towards the three regions studied in this paper. The mean δI uncertainties of the Stokes I maps at 850 µm are 2.4 mJy/14″-beam for W3(OH), 2mJy/l4″-beam for W3-Main, and 3.8 mJy/14″-beam for S106. The mean δPI uncertainties of the polarized emission maps are 1.5 mJy/14″-beam for W3(OH), 1.2 mJy/14″-beam for W3-Main, and 1.7 mJy/14″-beam for S 106. From the Stokes I maps at 850 µm and assuming a uniform temperature of T = 20 K, we calculate the column density $N_{{\text{H}}_2}$ maps (see details in 2.2). Fig. A.4 show the column density maps for the three targets with the POS B-field angles.

Appendix A.2 Selection of the filamentary emission

To remove the extended emission of the clouds and select the filamentary structures, we compute the Hessian matrix H of the second-order partial derivatives for all pixels of the Stokes I maps with a pixel size of 4″. The solution of the characteristic equation of the Hessian matrix provides two eigenvalues. The smaller of the two eigenvalues, λ⁻, corresponds to the minimum curvature of the intensity map (Arzoumanian 2012; Schisano et al. 2014; Planck Collaboration Int. XXXII 2016). The minimum curvature maps derived from the Stokes I are shown in the middle panels of Figs. A.1, A.2, and A.3, for W3(OH), W3-Main, and S 106, respectively. Filaments are elongated structure which corresponds to a relatively small curvature along their main axis. Consequently, filamentary structures are enhanced in the minimum curvature maps with respect to the surrounding extended emission. To derive “filament maps”, we further mask the pixels of the λ⁻ maps with I/δI < 5 (cf. the right hand side panels of Figs. A.1, A.2, and A.3).

Fig. A.1 Left: Stokes I map at 850 µm observed with the JCMT POL-2 towards W3(OH) in unit of mJy beam−1. The HPBW resolution of the maps shown on this figure is 14″ and the pixel size is 4″. The red contours correspond to I/δI > 5. The cyan and blue contours correspond to PI/δPI > 1 and PI/δPI > 5, respectively. The black dashed circle have a diameter of 12′. Middle: The minimum curvature map (in unit of mJy beam−1) of W3(OH) derived from the Stokes I map. Right: W3(OH) minimum curvature map (in unit of mJy beam−1) of W3(OH) where pixels with I/δI < 5 are masked. The filamentary structures are enhanced by the minimum curvature map. The black rectangle is the field shown in Fig. 1. The dashed circle is the same as on the left panel.

Fig. A.2 Same as Fig. A.1 for W3-Main. The black rectangle on the right hand side panel corresponds to the area shown in Fig. 2.

Fig. A.3 Same as Fig. A.1 for S106. The black rectangle on the right hand side panel corresponds to the area shown in Fig. 3.

Fig. A.4 Column density map, in units of 1021 cm−2, of W3(OH), W3 Main and S 106 from left to right. The contours are 1.4, 4.6, 23, 46, and 231 ×1021 cm−2 equivalent to the Stokes I contours shown in Figs. 1, 2, and 3. The blue and red lines show the ${\chi }_{B_{\text{POS}}}$ angle for 1 < PI/δPI < 5 and PI/δPI ≥ 5, respectively, and I/δI > 5. The lengths the lines are normalized to better represent the POS B-field orientation. The data are at an angular resolution of 14″.

Using the Hessian matrix H, we also compute θ_fil the local POS angle of the filamentary structures for all the selected pixels of the filament maps as follow: $${\theta }_{\text{fil}}=\frac{1}{2}\arctan \left(\frac{2H_{xy}}{H_{x^2}-H_{y^2}}\right)+90°,$$(A.1)

where the x and y axes are the horizontal and vertical axes of the Stokes I map in Cartesian coordinates, and $H_{x^2}={\partial }^{2}I/\partial x^2,H_{y^2}={\partial }^{2}I/\partial y^2$ and H_xy = H_yx = ∂² I/∂x∂y (see also Arzoumanian et al. 2019).

References

All Tables

All Figures

Fig. 1 Left: total intensity Stokes I map at 850 µm observed with the JCMT SCUBA-2/POL-2 towards W3(OH) in units of mJy beam−1. The contours are I = 30, 100, 500, 1000, and 5000 mJybeam−1. The lowest contour of I = 30 mJybeam−1 is equivalent to I/δI ~ 13. The red lines show the orientation of the POS B-field angle (${\chi }_{B_{\text{POS}}}$) for I/δI > 5, PI/δPI > 1, and PF < 25%. The lengths of the red lines are proportional to PF. A line showing PF = 10% is indicated on the plot. The data are at an angular resolution of 14″ or ~0.14 pc at the 2 kpc distance of the source. Right: contours and red lines are the same as in the left panel. The red lines are plotted with the same length to better show the B-field orientation. The background is a colour image composed using Spitzer 3.6 µm, 4.5 µm, and 8.0 µm images for the blue, green, and red, respectively.

In the text

	Fig. 2 Same as Fig. 1 for W3-Main. The data are at an angular resolution of 14″ or ~0.14 pc at the 2 kpc distance of the source.
In the text

	Fig. 3 Same as Fig. 1 for S 106. The data are at an angular resolution of 14″ or ~0.09 pc at the 1.3 kpc distance of the source.
In the text

Fig. 4 Left: blown-up image of the 850 µm Stokes I emission towards the hub of W3 Main (in units of Jy/beam), overlaid by normalized POS B-field lines. Grey contours display VLA 6 cm emission at 0.1″ resolution; the symbols A, B, C, and D mark the known compact HII regions. The youngest most luminous source IRS 5 is also indicated. The white arcs aid visualization of the organized B-field pattern. The white star and arrow indicate the centre of a CO outflow and its direction, respectively (see text for details). Right: W3 foreground mini-HFS, also referred to as W3 Main SE, indicated with the square on the left panel. The Spitzer three colour image is overplotted with Stokes I contours and normalized POS B-field line.

In the text

Fig. 5 Hub regions of W3(OH), W3-Main, and S106 (clockwise from top left). The colour compositions are the same as the right panels of Figs. 1, 2, and 3. Here the cyan contours trace column densities of 1022, 1023 , and 5×1023 cm−2. The black ellipses represent the hubs as described in Table 1, and their centres are marked with a black box-circle symbol. The blue contours represent VLA 4.89 GHz continuum emission obtained in B/C configurations. This emission is shown with a beam size of 0.7″ for W3(OH), however, it is compact at the scale of 0.08″. In the ‘two node’ system, the OpCl* (Carpenter et al. 2000) represents the more evolved node compared to the younger node composed of W3(H2O) and W3(OH) (e.g. Qin et al. 2016; Wyrowski et al. 1999). In W3-Main, IRS 5 surrounded by compact HII regions (Tieftrunk et al. 1997) A & B represent the younger node, while the evolved node corresponds to the compact regions C & D. In S 106, the younger and older nodes correspond to IRS 4 and the infrared cluster, respectively (Saito et al. 2009).

In the text

Fig. 6 Left to right: circularly averaged radial profiles of column densities and associated power-law fits displayed for targets W3(OH), W3 Main, and S106. The shaded region indicates the standard deviation of the values in each annulus. The vertical dotted lines are marked at 0.6, 0.75, and 0.8 pc, respectively for the left to right panels, representing the approximate transition point between power laws. In W3-Main and S 106, the centre of the hub is located in between the ‘two nodes’ (see Sect. 4.1 for details and Fig. 5).

In the text

	Fig. 7 Maps of filament line mass, Mline (in M⊙ /pc) for W3(OH), W3 Main, and S 106 from left to right. All three panels have the same physical scale of 5.8 pc × 5.8 pc and the same colour scale. The black ellipses indicate the hub sizes as given in Table 1.
In the text

	Fig. 8 Line mass histograms for the full image (black dashed line), filament (blue), and hub (red) regions for W3(OH), W3-Main, and S 106 from left to right. The fitted power laws are shown by blue (filament) and red (hub) lines, and the fitted slopes are indicated in each panel.
In the text

	Fig. 9 Polarized intensity (PI) histograms for the total (black), hub (red), and filament (blue) regions for W3(OH), W3-Main, and S 106 from left to right. The vertical axes indicate the number of pixels per bin.
In the text

	Fig. 10 Polarization fraction (PF) histograms for the total (black), hub (red), and filament (blue) regions for W3(OH), W3-Main, and S 106 from left to right. The vertical axes indicate the number of pixels per bin.
In the text

Fig. 11 Clockwise from top left: spatial correlation between total intensity, PF, and PI in W3(OH), W3 Main, and S 106 displayed over Herschel SPIRE 250 µm image. Black contours represent 1/PF with levels 0.25, 0.5, 1, 1.5, and 2 (for W3(OH) and W3 Main, top) and 0.1, 0.3, 0.6, and 0.9 (for S 106, bottom). White contours display PI with levels from 20 to 90 mJy/beam in intervals of 10 mJy/beam for W3(OH) and W3 Main (top), and levels 5, 10, 20, 30, 40, 50, and 60 mJy/beam for S 106 (bottom). The cyan vectors in the top two panels mark bipolar patterns in PI that correspond to outflows (see Sect. 5.2). The blue vector in W3 Main identifies the low PF elongated structure centred on IRS 5. This is similar to the elongated disc structure in W3(OH) traced by 1/PF contours.

In the text

	Fig. 12 Maps of the difference between the position angle of the filaments θfil and the ${\chi }_{B_{\text{POS}}}$ angle for W3(OH), W3-Main, and S 106 from left to right. The black ellipses indicate the hub sizes as given in Table 1 and shown in Fig. 7.
In the text

	Fig. 13 Histograms of the difference between the filament orientation and ${\chi }_{B_{\text{POS}}}$ for the total (black), hub (red), and filament (blue) regions, for W3(OH), W3-Main, and S 106 from left to right. The vertical axes indicate the pixel number. Here, pixels with the error on ${\chi }_{B_{\text{POS}}}$ smaller than 20 degrees are chosen.
In the text

Fig. A.1 Left: Stokes I map at 850 µm observed with the JCMT POL-2 towards W3(OH) in unit of mJy beam−1. The HPBW resolution of the maps shown on this figure is 14″ and the pixel size is 4″. The red contours correspond to I/δI > 5. The cyan and blue contours correspond to PI/δPI > 1 and PI/δPI > 5, respectively. The black dashed circle have a diameter of 12′. Middle: The minimum curvature map (in unit of mJy beam−1) of W3(OH) derived from the Stokes I map. Right: W3(OH) minimum curvature map (in unit of mJy beam−1) of W3(OH) where pixels with I/δI < 5 are masked. The filamentary structures are enhanced by the minimum curvature map. The black rectangle is the field shown in Fig. 1. The dashed circle is the same as on the left panel.

In the text

	Fig. A.2 Same as Fig. A.1 for W3-Main. The black rectangle on the right hand side panel corresponds to the area shown in Fig. 2.
In the text

	Fig. A.3 Same as Fig. A.1 for S106. The black rectangle on the right hand side panel corresponds to the area shown in Fig. 3.
In the text

Fig. A.4 Column density map, in units of 1021 cm−2, of W3(OH), W3 Main and S 106 from left to right. The contours are 1.4, 4.6, 23, 46, and 231 ×1021 cm−2 equivalent to the Stokes I contours shown in Figs. 1, 2, and 3. The blue and red lines show the ${\chi }_{B_{\text{POS}}}$ angle for 1 < PI/δPI < 5 and PI/δPI ≥ 5, respectively, and I/δI > 5. The lengths the lines are normalized to better represent the POS B-field orientation. The data are at an angular resolution of 14″.

In the text