© The Authors 2026

1 Introduction

Protoclusters, or embedded star clusters, are gas-dominated regions in which stars are actively forming. Unlike clusters, where gravity is dominated by the stars themselves, protoclusters are defined by the gravitational effect of the dense gas from which the stars emerge (Stutz 2018). It is crucial to study these regions for understanding the early stages of star cluster formation, particularly the formation of high-mass stars, which significantly affect their surroundings through outflows, radiation, and supernovae, and which predominantly form in clusters (Motte et al. 2018). These massive stars play a vital role in shaping galaxy evolution and estimating star formation rates, highlighting the importance of studying protoclusters and their internal assembly processes to understand key astrophysical phenomena.

Molecular gas is often organized in filaments: elongated and dense gas structures that are the birth sites of prestellar and protostellar cores. These cores are embedded within the filaments and represent the sites of star formation. They provide a natural link for investigating the kinematic relation between dense gas and forming protostars. Filaments can be detected at the scale of molecular clouds, but also at small protocluster scales (<1 pc); they constitute the immediate gas reservoir from which cores accrete (e.g., Kirk et al. 2013; Stutz & Gould 2016; Zhou et al. 2022; Hacar et al. 2023; Louvet et al. 2024). In cluster-forming environments, filaments are not expected to be static structures. High-mass star formation (HMSF) environments are characterized by deep gravitational potentials, high gas surface densities, stellar feedback, and strong gas motions (e.g., inflow, turbulence, and outflow), all of which can significantly modify the kinematic state of the filamentary gas. Observational signatures encoded in molecular spectra such as longitudinal and transverse velocity gradients (e.g., Álvarez-Gutiérrez et al. 2021, 2024; Sandoval-Garrido et al. 2025), velocity-coherent substructures (e.g., González Lobos & Stutz 2019; Cunningham et al. 2023), or multiple velocity components (e.g., Csengeri et al. 2011a,b) can trace gas-accretion flows inside filaments (e.g., Kong et al. 2019; Liu et al. 2023), gravitational collapse, rotational motions, or the effect of ionizing feedback and gas ejection (e.g., Galván-Madrid et al. 2010, 2024; Towner et al. 2024; Armante et al. 2024). It is particularly important to distinguish between these scenarios in massive protoclusters, where different theoretical frameworks predict different roles for filamentary gas. For example, models of global hierarchical collapse suggest that filaments funnel material toward cluster centers through gravitation-driven inflows (e.g., Vázquez-Semadeni et al. 2019), while alternative scenarios emphasize converging turbulent flows (e.g., Padoan & Nordlund 2002) or rotationally supported structures as mechanisms shaping the filament dynamics (e.g., André et al. 2014). Therefore, the study of the kinematics of dense gas in filaments at protocluster scales provides a direct way to test how mass is assembled and redistributed during the formation of stellar clusters.

To study these processes in a massive cluster-forming environment, we focused on the G012.80 protocluster (hereafter G012), which is one of the 15 nearby massive protoclusters observed by the Atacama Large Millimeter/submillimeter Array (ALMA) Initial Mass Function (IMF) Large Program¹ (ALMAIMF LP, Motte et al. 2022). Specifically, G012 is an active and massive (M_{870 μ}m = 1.7 × 10³ M_⊙, see Motte et al. 2022) star-forming region centered on the W33 main clump, with a trigonometric parallax distance of 2.4 kpc (Immer et al. 2013) and a local standard of rest velocity of V_LSR = 37 km s⁻¹ estimated from maser velocities (Immer et al. 2014). A central star cluster with spectral types from O7.5 to B1.5 has been reported (Immer et al. 2013), along with evidence of dissociation of complex molecules on small scales (e.g., Immer et al. 2014). Detections of HC₃N, CO outflows, class I methanol masers, H``II``` regions, OB star clusters, and a high abundance of cores indicate that G012 is an optimal testbed for studying HMSF in an evolved protocluster, particularly the mechanisms regulating the distribution and kinematics of dense gas in a feedback-rich environment (e.g., Haschick & Ho 1983; Yu et al. 2019; Xie et al. 2023; Armante et al. 2024). To probe the dense gas kinematics in G012, we focused on diazenylium (N₂H⁺; e.g., Bergin & Langer 1997; Caselli et al. 1995; Bergin & Tafalla 2007; Lippok et al. 2013; Tatematsu et al. 2008; Busquet et al. 2011; Gómez et al. 2022). This nitrogen-bearing molecule is characterized by lower depletion levels onto dust grains and was first detected in the interstellar medium by Thaddeus & Turner (1975). For the N₂H⁺(1−0) transition, the critical density ranges from 2.0 × 10⁴ cm⁻³ at 100 K to 6.1 × 10⁴ cm⁻³ at 10 K (Shirley 2015), indicating the conditions under which this line becomes an effective tracer. Variations in the observed N₂H⁺ emission can therefore be linked to changes in the physical conditions of the protocluster (Tobin et al. 2013; Tafalla et al. 2021; Yu et al. 2022).

When the kinematics and dense gas emission traced by molecules like N₂H⁺ are scrutinized in individual protoclusters, multiple evolutionary stages (or generations of stars) emerge within a single region (Cunningham et al. 2023; Pouteau et al. 2023; Armante et al. 2024). This motivates detailed studies of evolved protoclusters with internal structures at potentially different evolutionary stages. The search for correlations between the kinematic properties of these dense gas structures (e.g., Stutz & Gould 2016; Álvarez-Gutiérrez et al. 2021; Xu et al. 2023; Álvarez-Gutiérrez et al. 2024; Reyes-Reyes et al. 2024; Sandoval-Garrido et al. 2025) and their line-mass profiles (e.g., Stutz 2018) sets the stage for developing observationally driven physical models of protocluster evolution as stellar mass is assembled. In this context, we focused on the two dominant internal cold dense gas structures in G012 identified with N₂H⁺ (1−0). These structures have the form of coherent filaments that show evidence of either rotation or collapse. We characterize the dense gas kinematics and key physical parameters in the protocluster, with particular emphasis on these two main filaments.

This paper is organized as follows. In Sect. 2 we detail the ALMA-IMF dataset and core catalogs for the G012 protocluster. In Sect. 3 we analyze the N₂H⁺ moment maps, PV diagrams, column density, and N₂H⁺ core velocities in the region. We perform a detailed analysis of the two dominant cold dense gas filamentary structures in Sect. 4, focusing on average velocity gradients and line-mass profiles. In Sect. 5, we apply a rotation toy model to the observed dense gas kinematics in the R1 filament. In Sect. 6, we discuss potential evolutionary scenarios by comparing the kinematics and other observables (e.g., core incidence). Finally, we summarize our main conclusions in Sect. 7.

2 Data

We used observations from the ALMA telescope provided by the ALMA-IMF LP. The line parameters for N₂H⁺ and additional tracers are summarized in Table 1.

2.1 ALMA-IMF datacubes

We employed N₂H⁺(J=1−0) observations at a frequency of 93.1734 GHz. We followed a similar data reduction procedure as Álvarez-Gutiérrez et al. (2024) and Sandoval-Garrido et al. (2025). N₂H⁺ observations were cleaned using the ALMA-IMF imaging pipeline² (Cunningham et al. 2023) and the version 5.6.0 of The Common Astronomy Software Applications package (CASA, CASA Team 2022). We used the imcontsub task to subtract the continuum emission from the N₂H⁺ line emission. To estimate the continuum emission, we only considered emission-free channels in the range of 12−20 km s⁻¹, where we applied a linear fit using fitorder = 0. To recover the cloud emission at all available scales, we combined observations of the 7m and 12m arrays with total power data (TP). To do this, we used the task feather. This combination also allowed us to recover the missing flux, visible as negative bowls in the inter-ferometric spectral data produced by the lack of zero-spacing. We obtained a fully integrated multi-scale emission dataset. The resulting N₂H⁺ datacube contained a total of 217 velocity channels over a range from 12 km s⁻¹ to 61 km s⁻¹, with a spectral resolution of 0.23 km s⁻¹ and a final beam size of ~2.3" (see Table 1), corresponding to ~5.5 kau at the source distance. In Fig. 1 we display Spitzer observations at 8 μm (red), 4.5 μm (green) and 3.6 μm (blue) along with an N₂H⁺ integrated-intensity contours at 25 and 100 K km s⁻¹ (left panel) and N₂H⁺ integrated-intensity map at a signal-to-noise ratio (S/N) > 12 (right panel; see Appendix A for the S/N selection criteria). The N₂H⁺ integrated-intensity map reveals a filamentary and clumpy morphology that traces the dense gas in the region, with a clear absence of emission toward the protocluster center, consistent with the potential chemical destruction of N₂H⁺ in the vicinity of the young stellar cluster (see Sect. 6.1).

In addition to the N₂H⁺ data, we used complementary spectral lines previously analyzed in ALMA-IMF studies. These include the H41α recombination line (Galván-Madrid et al. 2024), which traces HII regions associated with ionizing sources in the protocluster center; DCN(J=2−1; Cunningham et al. 2023), which traces more compact and warmer emission than N₂H⁺; SiO (J=5−4; Towner et al. 2024), used to analyze potential outflow features; and C¹⁸O (J=2−1; Koley et al. 2025), which traces compact emission in the central region of the protocluster and more extended and fainter emission in the surroundings.

Fig. 1 Multiwavelength view and molecular gas distribution of the G012 protocluster. Left panel: Spitzer RGB composite figure of the G012 protocluster at 8 μm (red), 4.5 μm (green), and 3.6 μm (blue). The two cyan contours trace the N2H+ integrated-intensity emission at 25 and 100 K km s−1, respectively. Right panel: N2H+ integrated-intensity map with contours at 25 and 100 K kms−1, corresponding to an S/N > 12, shown in black (same levels as in the left panel). The blue boxes highlight the two main filamentary structures, R1 and R2, with estimated lengths of 0.56 pc for both. We show the average spectra of the two regions in the insets. The black curve represents the data within the black boxes. The red markers represent ionizing regions (crosses), radio sources (plusses), and infrared sources (dots) detected in the region by Haschick & Ho (1983). The green contours show continuum emission from the Atacama Large Millimeter/submillimeter Array Band 6 observations (Ginsburg et al. 2022b). The black ellipse in the bottom right corner represents the beam size of the N2H+ data.

2.2 ALMA-IMF dense core catalogs

Louvet et al. (2024) constructed the first core catalog of the full ALMA-IMF sample based on the 1.3 mm and 3 mm continuum emission. They used the toolkit Getsf (Men'shchikov 2021) to identify compact continuum sources. In G012, they detected 57 dense cores with estimated masses of 0.6–6.2 M_⊙ and temperatures of 25–100K. As a follow-up, Cunningham et al. (2023) used the DCN molecular line emission to study the core population from Louvet et al. (2024). For G012, they were able to estimate spectral parameters of 38 cores, including their line-of-sight velocities. Deeper, Armante et al. (2024) constructed a catalog of G012 cores complementing dust continuum emission (at 1.3 mm and 3 mm) with molecular emission from the ¹²CO(2−1), SiO(5−4), CH₃OCHO, and CH₃CN spectral lines. They classified the cores as prestellar or protostellar, finding a total of 94 detections. Additionally, Motte et al. (2025) provided new mass, temperature, and luminosity estimates for cores in the remaining ALMA-IMF protoclusters using radiative transfer modeling. For G012, they found protostellar and prestellar cores distributed over mass ranges of 0.3–9.5 M_⊙.

To complement these data and provide a more complete characterization of G012, we used all these core catalogs in our analyses. In particular, we estimate in the Appendix F new velocities for 48 cores from these previous catalogs.

3 Analysis of the N₂H⁺ data cube

By modeling the hyperfine structure of the N₂H⁺ (as described in Appendix A), we derived key kinematic parameters and simultaneously identified multiple velocity components in the N₂H⁺ spectra. We adopted a similar fitting procedure as Sandoval-Garrido et al. (2025) to model the N₂H⁺ emission. We used the hyperfine line structure model provided by version V.1.0.1 of the PySpecKit spectroscopic analysis toolkit (Ginsburg et al. 2022a). The vast majority of spectra that we inspected are well characterized by either one or two velocity components (see Appendix A). Hence, we applied a fitting model with up to two components, using parameters such as excitation temperature (Tex), optical depth (τ), velocity centroid (Vc), and velocity dispersion (σ). The model-fitting results showed that G012 is characterized by two main velocity structures separated by ~2.5 km s⁻¹, which we defined as the first velocity component (FVC) and the second velocity component (SVC). We adopted 35.6 km s⁻¹ as the reference velocity (see Appendix A) to separate them: FVC refers to spectra with velocities below this value, and SVC refers to those above it. Using the FVC and SVC structures, we further explored the kinematics of G012. First, we reviewed moment maps, which provide a detailed look at the velocity centroid and dispersion patterns in the region and of the particular filamentary structures R1 and R2 (see Fig. 1). We then constructed PV diagrams of the entire region to capture the velocity gradients and dynamical behavior of the protocluster. We also estimated column densities and masses in G012 in order to set constraints on the gas mass and density profiles in the region.

3.1 Moment maps

In Fig. 2, we display the N₂H⁺ integrated-intensity (moment zero map, upper panels), velocity center (moment one map, middle panels), and velocity dispersion (moment two map, bottom panels) maps of the FVC and SVC. In the N₂H⁺ integrated-intensity map, G012 is characterized by a filamentary distribution in the plane of the sky (POS). Specifically, we highlight the two prominent filaments with blue boxes (see also Fig. 1), which are the focus of this study, hereafter called R1 and R2. These structures were identified based on the visual inspection of prominent and kinematically well-separated filamentary structures, where the kinematic patterns described below are more clearly distinguished. However, to complement this more visual inspection, we also applied a more automated method, as described in Appendix B.

The R1 filament is largely dominated by the FVC (see Fig. 2, left panel), which displays one pronounced velocity gradient across the filament. The R1 mean centroid velocity and line width are 〈V〉 = 34.4 km s⁻¹ and 〈σ〉 = 0.81 km s⁻¹, respectively. In contrast, the R2 filament is dominated by the SVC (see Fig. 2, right panel), which exhibits more homogeneous velocity structures. The mean centroid velocity and line width for R2 are 〈V〉 = 36.8 km s⁻¹ and 〈σ〉 = 0.75 km s⁻¹, respectively (see Table 2 for global parameters of the main filaments). Additionally, at the center of G012, we observe an absence of N₂H⁺ emission that may be related to the destruction of the molecule in the region (see Sect. 6). To analyze this absence, we reviewed the integrated-intensity maps of the fitted complementary tracers in Fig. A.3. The C¹⁸O emission (upper left panel) traces more extended gas than N₂H⁺ and is primarily distributed around the H41α bubbles (bottom left panel). In addition, we identified high C¹⁸O integrated-intensity structures in the surroundings of the R1 and R2 filaments. The SiO integrated intensity reveals elongated features in the R2 filament (upper right panel). Specifically, one perpendicular outflow feature is observed at the top of R2, related to a previously detected hot core (Armante et al. 2024). Furthermore, DCN, which is associated with increased star formation, traces some of the densest parts of the filamentary structures in R1 and R2 (bottom right panel).

3.2 N₂H⁺ PV diagrams

We used the technique developed in González Lobos & Stutz (2019) to construct the N₂H⁺ intensity-weighted position velocity diagrams (e.g., Álvarez-Gutiérrez et al. 2021, 2024; Sandoval-Garrido et al. 2025). To obtain better-defined structures in PV space, we spatially rotated the N₂H⁺ data cubes (45^° with respect to the protocluster center) by vertically aligning the most predominant filaments. In Fig. 3, we show PV diagrams of the first velocity component (blue shades) and the second velocity component (green shades) structures. The upper left panel shows the integrated-intensity map after spatial alignment. In the upper right and bottom left panels, we display the PV diagrams perpendicular and parallel to the region, respectively. Additionally, we include the positions of previous cores detected in G012 (Cunningham et al. 2023; Louvet et al. 2024; Armante et al. 2024; Motte et al. 2025). For the PV diagrams, we show a representative timescale of 0.2 Myr³ based on the PV structures that are more strongly spread in velocity (ΔV > 3 km s⁻¹). In addition to the main features mentioned above, we also highlight an elongated structure with high integrated intensity in the bottom left diagram, located at a velocity range (ΔV) of approximately −2 to 0 km s⁻¹. This structure is spatially related to the R1 filament (see Sect. 4). The key features in the PV diagrams are listed below.

1. R2 reveals a set of characteristic twisting features (Fig. 3, bottom right panel) that are confined to a narrow velocity range (< 2 km s⁻¹), in contrast to the broader velocity spread observed in R1.

2. The PV features are more spread out in velocity (ΔV > 3 km s⁻¹ ; see the upper right panel in Fig. 3) in regions near H41α and SiO emission (see Fig. A.3), indicating short associated timescales. This potentially highlights the effect of stellar feedback close to the center of G012.

3. The R1 and R2 filaments reveal the most predominant and intricate PV structures.

4. The R1 filament exhibits a double-helix feature (Fig. 3, bottom right panel), spanning approximately ~3 km s⁻¹ in velocity. This type of feature has been associated with filament rotation in other regions (e.g., Álvarez-Gutiérrez et al. 2021). We address this point in Sect. 5.

Fig. 2 N2H+(1−0) integrated intensity (upper panels), velocity centroid (middle panels), and velocity dispersion (bottom panels) of the FVC (left panels) and SVC (right panels). The black contours trace the N2H+ integrated-intensity emission at 25 and 100 K km s−1. The black boxes in all panels display the spatial location of the R1 and R2 filaments. The black ellipse in the bottom right corner represents the beam size of the N2H+ data.

3.3 N₂H⁺ column density

The column density and mass in protoclusters provide insights into the amount of material available for star formation. In order to estimate the N₂H⁺ column density, we used the T_ex, τ, and σ PySpecKit output parameters (see Appendix A). We applied the column density approximation outlined in Caselli et al. (2002a) and Redaelli et al. (2019), which is described as $$\text{N}\left({\text{N}}_2{\text{H}}^{+}\right)=\frac{4{\pi }^{3/2}{v}^{3}Q\sigma \tau \text{\hspace{0.17em}exp}\left(\frac{E_u}{k_B{T}_{\text{ex}}}\right)}{\sqrt{\text{ln\hspace{0.17em}}2}{c}^{3}Ag\left(\text{exp}\left(\frac{hv}{k_B{T}_{\text{ex}}}\right)-1\right)}.$$(1)

Here, ν represents frequency of the N₂H⁺(1−0) emission, g = 3 is the statistical weight (degeneracy) of the upper energy level of the transition level, c is the light speed, A = 3.6e⁻⁵ s⁻¹ is the Einstein coefficient, E_u = 4.47 K corresponds to the upper energy level of the transition, Q is the rotational partition function of N₂H⁺, k_B is the Boltzmann constant, and h is the Planckconstant (Shirley 2015; Redaelli et al. 2019).

To mitigate potential biases introduced by poorly defined fitting parameters (see Appendix A) and to prevent anomalies in the column density and relative abundance estimates, we excluded all pixels affected by systematic biases in τ and T_ex. Specifically, we removed pixels where T_ex = 150 K (the fitting upper limit) and τ < 0.2, which resulted in discarding 30% of the affected pixels (see Appendix D). We consider the final column density map as the sum of the FVC and the SVC column densities. Overall, the R1 column density distribution appear more interconnected, compact, and filamentary. In contrast, the R2 column density values are distributed more homogeneously, taking high values of about 1 × 10¹⁴ cm⁻² in most of the region. We converted the N₂H⁺ column density into mass, obtaining a respective N₂H⁺ mass of 1.72 × 10⁻⁵ M_⊙ in the entire protocluster. The two main filaments, R1 and R2, contribute approximately 12% and 17% of the N₂H⁺ protocluster mass, respectively. Moreover, when considering the extended emission associated with these filaments (see Appendix B), their total contribution accounts for about 60 % of the total N₂H⁺ mass, further emphasizing their importance in the kinematics and mass distribution of the region.

Fig. 3 Position-position and position-velocity diagrams of the N2H+ observations. Upper left panel: G012 integrated intensity of FVC (blue shades) and SVC (green shades). We display the DCN cores from Cunningham et al. (2023) catalog with black pluses. The Armante et al. (2024) core catalog is divided into four categories: cores only detected with N2H+ (red crosses), cores only detected with DCN data (red circles), cores detected with N2H+ and DCN (red triangles), and cores detected neither in N2H+ nor DCN (black triangles). The symbol sizes are proportional to the estimated core mass. The black ellipse in the bottom right corner represents the beam size of N2H+ data. Upper right panel: PV diagram along the y-axis; the velocity axis is subtracted from the N2H+ systemic VLSR of the protocluster (35.5 km s−1). The black arrow illustrates a slope of 5 pc (km s−1 )−1, which corresponds to a timescale of ~0.2 Myr, that is, the timescales that approximately correspond to some of the extended structures in this PV diagram. At the top of R1, lies a wrapped (or double-helix) type velocity field with spreads of more than ~3 km s−1. On the other hand, R2 appears very compact in velocity along its extent, with small-scale spatial wiggles. Bottom left panel: PV diagram in the perpendicular direction compared to the upper right panel. This shows the emergence, albeit somewhat hidden in the overall velocity field, of an approximately uniform and extended gradient in R1 that is apparent near ΔX ~ 0.60 to 0.75 pc, ΔV ~ −2 to 0 km s−1. In contrast, R2 appears as the compact blob of cores (triangles and x-symbols) on the right-hand-side of the panel, characterized by the absence of an obvious gradient in position and velocity. Bottom right panel: zoomed PV diagrams of the main filaments R1 (left panel) and R2 (right panel), enclosed in the black boxes of the upper left panel. R1 presents a wrapping double-helix signature that is most obvious toward the top of the diagram and which is dominated by the FVC velocities. R2 exhibits comparatively compact velocity variations (ΔVmax~1.5 km s−1) along the filament and contains a high number of massive cores (1–3 M⊙).

Fig. 4 Left panel: relative abundance map in G012. The black contour highlights the mask of $\frac{\tau }{e\left(\tau \right)}>2$ applied to the N2H+ column density (see Sect. 3.4); most of the pixels removed by this mask do not affect the main filaments significantly. The areas lacking data (top of R1 and bottom of R2) are due to the H2 column density map lack of coverage in the filaments. The representative relative abundance value is 0.93 × 10−10. The black circle in the bottom left corner represents the beam size of the H2 map. Right panel: relative abundance histogram of the values inside the black contour in the left panel. The dashed red line represents the mode of the distribution. In the upper-right corner, the values of the mode and standard deviation are shown, which we consider as the representative value and error of the sample.

3.4 Relative abundance

In G012, previous estimates of the H₂ column density map (Dell'Ova et al. 2024) did not provide the full spatial coverage achieved by the N₂H⁺ map. Nevertheless, despite these limitations, particularly the incomplete coverage of the main filaments in the region, the map from Dell'Ova et al. (2024) remains the most suitable tool currently available for estimating the relative abundance in the region. We used the N₂H⁺ column density to estimate the total mass (M_tot) in the protocluster as a whole and in areas where the H₂ map lacks complete spatial coverage. First, we calculated the relative abundance in the N₂H⁺ and H₂ column density maps (see Fq. (2)) where we had coverage in both. To estimate this ratio, we reprojected the N₂H⁺ image to match the pixel scale of the H₂ data, which is slightly larger (with a pixel scale of 0.83"). Subsequently, we calculated the relative abundance map by taking the pixel-to-pixel ratio of the two maps, $$X\left(N_2{H}^{+}\right)=\frac{N\left(N_2{H}^{+}\right)}{N\left(H_2\right)}.$$(2)

Here, N(N₂H⁺) is the N₂H⁺ column density map described above, and N(H₂) is the H₂ column density map from Dell'Ova et al. (2024). For the N(H₂) map, we considered an S/N > 3 based on the ratio of the column density and the associated errors. In Fig. 4 (left panel) we present the resulting relative abundance map. We identify a clear trend in which the relative abundances values decrease toward regions associated with ionization (see Fig. 8 left panel, for a comparison between the N₂H⁺ emission and the H41α distribution). This behavior likely reflects the effect of stellar feedback on the molecular gas distribution. We discuss this point further in Sect. 6.

We determined representative N₂H⁺ relative abundance values using the mode of its distribution, a choice motivated by its stability under different binning and masking conditions (see Appendix D). We obtained a relative abundance of X(N₂H⁺) = (0.93 ± 0.10) × 10⁻¹⁰ (see Figure 4, right panel), consistent with previous determinations (Caselli et al. 2002c; Sandoval-Garrido et al. 2025). This representative value was then applied to derive total H₂ mass. Since these estimates will underpin the line-mass profiles presented in Sect. 4.2, we first validated our choice of mass map by calculating the line-mass profile in a section of the R2 filament covered by the N₂H⁺ and Dell'Ova et al. (2024) column density maps (see Appendix F). We found variations of less than 28% between the two methods, indicating that the use of N₂H⁺ to re-estimate the total H2 mass, and consequently, the line-mass profiles discussed in Sect. 4.2, is robust.

Fig. 5 Average velocity gradients perpendicular to R1 (left panel) and R2 (right panel) filaments. In both diagrams, the ΔV and Δr axes have the same range, allowing for a direct comparison of the slopes between the panels. In R1 the red line represents the linear fit weighted by the integrated-intensity points of the velocity gradient (VG). The associated VG timescale (τ) estimate is displayed in the upper left corner. In contrast, R2 lacks a clear velocity gradient structure as in R1 and is instead comparatively compact in its velocity distribution. In addition, we identified two different velocity distributions that spatially correspond to different regions in the R2 filament: 1) the most compact and central structure (compact feature) is related to the densest region of the filament; and 2) the most elongated feature is spatially related to the eastern edge of the R2 filament, characterized by elongated SiO emission. In the diagram, the lower intensity and scattered points represent regions surrounding the filament.

4 R1 and R2 characterization

The R1 and R2 filaments are the predominant N₂H⁺ integrated-intensity structures in the G012 protocluster. These regions present different velocity structures in PV diagrams (see Sect. 3), suggesting the presence of two different star-forming environments taking place in the same protocluster. In this section, we focus on the characterization of these two filaments with the aim to set constraints on their star-forming stage.

4.1 Average filament velocity gradients

As shown in Fig. 3 (bottom right panel), the R1 filament is associated with a potential signature of rotation, while R2 presents a relatively compact velocity distribution. We analyzed the velocity gradients perpendicular to the two filaments through average-velocity gradients estimations. We adapted the method described by Álvarez-Gutiérrez et al. (2021) for the filament rotation analysis and divided the procedure into the following steps:

1. We determined the total integrated intensity in each region (i.e., the moment 0 of the FVC and SVC, see Sect. 3) and calculated peak integrated-intensity values along the long axis of the two filaments in order to determine a representative filament ridge line (see Appendix C).

2. For each slice along the y-axis, we straightened the two filaments by centering the peak integrated intensity from the ridge line to a fixed arbitrary position. Then, each point of the integrated intensity, velocity centroid, and velocity dispersion maps had a projected radius to the center of the filament.

3. We generated ridge line-averaged PV diagrams as follows. For each point, at each radius from the ridge line, we plot the radial velocities as shown in Fig. 5. Here, the color of each point indicates their integrated intensity. These velocity versus radius diagrams therefore capture PV structures integrated over the length of the filament. If a given filament has a prominent velocity gradient approximately perpendicular to the ridge line, these diagrams will reveal these structures. Since the projected radius is measured perpendicular to the filament ridge line, a coherent velocity gradient in these velocity versus projected radius diagrams is naturally interpreted as rotation around the filament axis. This is the case for R1, while for R2, no such equivalently prominent structure is observed.

4. We applied a linear fit to the velocity gradient in the R1 filament using the integrated intensity as a statistical weight.

In Fig. 5, we show the R1 (left panel) and R2 (right panel) ridge line-averaged PV diagram described above. In this representation, a linear and elongated structure reflects a systematic change in velocity with projected radius, as expected from rotation around the filament axis. The global R1 velocity gradient (VG) has a magnitude of VG ~ 10.4 km s⁻¹ pc⁻¹ (with an associated statistical fitting error of ± 0.29 km s⁻¹ pc⁻¹). In R2, in sharp contrast to R1, we do not observe an equivalent global average velocity gradient. In contrast, R2 is dominated by a compact velocity structure, clumping at velocities around ~1 km s⁻¹ in the center of the diagram. This compact feature traces the central and densest part of the filament. We also observe a more spread-out structure in projected radius and velocity extending up and to the left in the diagram and ending near - 0.05 pc, 0.75 km s⁻¹ (elongated feature). This structure is approximately associated with the eastern edge of the R2 filament and has associated elongated SiO emission that likely traces protostellar outflows and shocks (Towner et al. 2024).

4.2 Line-mass profiles and associated 3D model quantities

The characterization of the mass distribution, and hence, of the gravitational potential, is a requirement for interpreting the kinematics in these systems since the gravitational field of the gas serves as governor of the dynamics of a system (e.g., Stutz 2018; González Lobos & Stutz 2019; Álvarez-Gutiérrez et al. 2021; Reyes-Reyes et al. 2024). We followed the formalism outlined in Stutz & Gould (2016) to estimate these line-mass profiles in the main filaments of the G012 protocluster. We first aligned the filament mass maps with respect to their ridge lines (see Appendix C). We used a box region to capture the densest parts of the two filaments (see Fig. 1) and constructed N₂H⁺ cumulative mass distributions for both regions. As expected for elongated and dense structures, the cumulative mass distribution profiles in R1 and R2 are nearly linear (see Fig. F.1), indicating that the mass distribution along the filament is approximately uniform. This allowed us to derive an averaged line-mass profile as a function of projected radius from the ridge line over the entire filament. While the two distributions share this linear behavior, they differ in normalization: R2 is systematically shifted toward higher masses than R1. Following the formalism in Stutz & Gould (2016), we found that the line-mass profiles are well-described by power laws of the form $$\lambda \left(\omega \right)=\zeta {\left(\frac{\omega }{\text{pc}}\right)}^{\gamma },$$(3)

where ω is the projected radius in the POS (or impact paramenter from the ridge line), γ corresponds to the index in the enclosed mass over length (M/L) versus projected radius diagram (see Fig. 6), and ζ is the M/L normalization constant. For R1 and R2, we obtained $$\begin{array}{c}{\lambda }_{\text{R}1}\left(\omega \right)=5660\frac{{\text{M}}_{\odot }}{\text{pc}}{\left(\omega /\text{pc}\right)}^{0.30};\\ {\lambda }_{\text{R}2}\left(\omega \right)=6943\frac{{\text{M}}_{\odot }}{\text{pc}}{\left(\omega /\text{pc}\right)}^{0.20}.\end{array}$$(4)

In Fig. 6, we display the M/L profiles of the R1 and R2 filaments (red and blue lines, respectively), along with other star-forming regions (dashed black lines) described below and in Table 3.

We followed Stutz & Gould (2016); Stutz (2018), and Álvarez-Gutiérrez et al. (2021) in estimating the density, gravitational potential, and acceleration assuming cylindrical 3D geometry. All quantities were estimated in the POS since we lack access to inclination information. We refer to these as "apparent" profiles following Álvarez-Gutiérrez et al. (2021). We estimated the apparent volume density as $$\rho \left(\text{r}\right)=\beta {\left(\frac{\text{r}}{\text{pc}}\right)}^{\gamma -2},$$(5)

where for the R1 and R2 filaments, we obtained β_R1 = 211.3 M_⊙ pc⁻³ and β_R2 = 181.3 M_⊙ pc⁻³. We estimated the gravitational potential as $$\varphi \left(\text{r}\right)=\psi {\left(\frac{\text{r}}{\text{pc}}\right)}^{\gamma },$$(6)

where, considering the R1 and R2 parameters described above, we estimated ψ_R1 = 123.1 km² s⁻² and ψ_R2 = 269.3 km² s⁻².

Finally, we estimated the gravitational acceleration as follows: $$g\left(\text{r}\right)=-\xi {\left(\frac{\text{r}}{\text{pc}}\right)}^{\gamma -1}.$$(7)

For R1 and R2, we obtained ξ_R1 =37.4 km² s⁻² pc⁻¹ and ξ_R2 = 51.3 km² s⁻² pc⁻¹, respectively.

Relative to other star-forming regions (see Fig. 6), the filaments in G012 exhibit significantly higher line-mass profiles that are up to three times more dense than extended regions such as the G351.77 protocluster (Reyes-Reyes et al. 2024). This trend remains even in more extended regions along the filaments, as detailed in Appendix B. The corresponding line-mass profile parameters, including the filament extensions, are listed in Table 3.

Fig. 6 Line-mass profile of the R1 and R2 filaments (red lines). The red shaded area represents the uncertainty in the two estimates, derived from the relative abundance uncertainty. This percentage therefore provides a lower limit to the error of 11% in both distributions. We include profiles of the Orion ISF (Stutz & Gould 2016), the ONC (Stutz 2018), the California L1482 region (Álvarez-Gutiérrez et al. 2021), and the G351.77 protocluster (Reyes-Reyes et al. 2024) with black lines. The main filaments of G012 reveal denser line-mass distributions than those found in other star-forming regions.

4.3 Star formation rate and efficiency

The star formation rate (SFR) and efficiency (SFE) are critical parameters for understanding the current star formation activity, offering insight into the efficiency with which gas is converted into stars. To estimate the SFR (Eq. (8)), we assumed lifetimes based on observational studies of pre- and protostellar cores. Specifically, the lifetimes of prestellar cores vary significantly with environment and initial conditions, ranging from 0.5 to 2 Myr in low-mass regions (e.g., Jessop & Ward-Thompson 2000; Merín et al. 2008; Evans et al. 2009; Könyves et al. 2015) to as short as 0.05–0.24 Myr in high-mass cores, where collapse is driven by fast-converging flows (Csengeri et al. 2011a,b; Motte et al. 2018; Valeille-Manet et al. 2025). As cores evolve into protostellar phases, accretion becomes dominant, and the lifetimes decrease with increasing mass from ~0.5Myr for low-mass protostars (Megeath et al. 2022) to 0.1–0.3 Myr for high-mass cores (Duarte-Cabral et al. 2013; Valeille-Manet et al. 2025). Gravitational contraction, temperature increase, and accretion rates are strongly affected by the initial core mass and local conditions (Bontemps et al. 2010; Dunham et al. 2014).

For the SFR and SFE estimates, we used the formalism outlined in Megeath et al. (2022) and considered the G012 core masses provided in Motte et al. (2025). In R1 and R2, Armante et al. (2024) and Motte et al. (2025) identified a total of 16 cores. Fourteen of these are located in R2, and only two are in R1. These cores are classified as either prestellar or protostellar, with masses ranging from 0.3 to 3 M_⊙. For our analysis, we adopted a prestellar timescale of 1.2 Myr, which is the expected value for prestellar cores in the mass range of 0.003–10 M_⊙ (Könyves et al. 2015). We also used a protostellar lifetime of 0.5 Myr, assuming a low-mass protostellar regime around 0.5 M_⊙ (Megeath et al. 2022). The two cores in R1 are classified as prestellar, while the R2 filament exhibits five cores in this phase. In addition, the R2 filament contains nine protostellar cores. Considering the pre- and protostellar lifetimes described above, we estimated the SFR as follows: $$\text{SFR}\left({\text{M}}_{\odot }{\text{Myr}}^1\right)=\frac{{\text{M}}_{\text{cores}}}{{\text{t}}_{\text{cores}}}.$$(8)

Here, M_cores is the mass of the core population (prestellar and/or protostellar) and t_cores is the approximate lifetime of the cores in million years. For R1 and R2, we obtained

• SFR_{R1–prestellar} = (4.24 ± 0.29) M_⊙Myr⁻¹,

• SFR_{R2–prestellar} = (12.1 ± 1.90) M_⊙Myr⁻¹,

• SFR_{R2–protostellar} = (43.2 ± 6.67) M_⊙Myr⁻¹,

• SFR_R2–total = SFR_prestellar + SFR_protostellar = (55.3 ± 6.93) M_⊙Myr⁻¹.

Errors were estimated solely from the mass uncertainties reported in the Motte et al. (2025) catalog, without possible uncertainties on the adopted timescales. For instance, a lifetime of 0.3 Myr for the intermediate- to high-mass protostellar regime would increase the SFR to (71.9 ± 11.2) M_⊙ Myr⁻¹.

In addition, we estimated the SFE as $$\text{SFE=}\frac{{\text{M}}_{\text{cores}}}{{\text{M}}_{\text{cloud}}}+{\text{M}}_{\text{cores}}.$$(9)

In this case, M_cores represents the same core mass as we used for the SFR estimate, and M_cloud corresponds to the R1 or R2 filament mass (see Table 2). In R1 and R2, we obtained

• SFE_{R1–prestellar} = (1.6 ± 0.1) × 10⁻³,

• SFE_{R2–prestellar} = (3.3 ± 0.4) × 10⁻³,

• SFE_{R2–prestellar} = (4.8 ± 0.4) × 10⁻³,

• SFE_R2–total = SFE_prestellar + SFE_protostellar = (8.1 ± 0.5) × 10⁻³.

In this case, we adopted the same assumptions for the core mass uncertainties as for the SFR. For the cloud (filament) mass, we assumed an 11% error, corresponding to the uncertainty associated with the relative abundance values in the region. Our estimates of SFR and SFE should be regarded as lower limits, since they are based exclusively on the population of pre- and protostellar cores traced by some specific tracers. Unlike nearby star-forming regions, more distant regions such as G012 require a more complete core census for accurate determinations. Overall, the SFR and SFE in R1 are lower by roughly an order of magnitude than in R2. This trend even persists in more extended filamentary regions (see Appendix B).

5 A filament rotation toy model applied to R1

In the N₂H⁺ PV diagram, the R1 filament shows a helical pattern. This feature is also observed in the California L1482-S cloud and was attributed to filamentary rotation (Álvarez-Gutiérrez et al. 2021; Hsieh et al. 2021). In Fig. 3 (bottom right panel) we highlight this feature at the top of R1, which is also associated with a clear average-velocity gradient of ~ 10.4 km s⁻¹ pc⁻¹ (see Fig. 5, left panel).

To determine whether rotation or gravity is the dominant factor in the R1 filament, we compared the velocity gradient to the gravitational acceleration implied by the filament mass distribution. We applied the formalism from Álvarez-Gutiérrez et al. (2021) to estimate the centrifugal force associated with rotation of the R1 filament. We then calculated the ratio of the centrifugal and gravitational forces as follows: $$\begin{array}{cc}\frac{{\text{F}}_{\text{c}}}{{\text{F}}_{\text{g}}}=\frac{{\text{VG}}^2{\text{pc}}^{-1}}{\xi {\text{\hspace{0.17em}cos}}^3\left(\theta \right)}{\left(\frac{\text{r}}{\text{pc}}\right)}^{-\gamma };& \text{VG\hspace{0.17em}}=\text{\hspace{0.17em}10}.4\frac{\text{km}}{\text{s\hspace{0.17em}pc}},\xi =37.4\frac{{\text{km}}^2}{{\text{s}}^2\text{pc}}.\end{array}$$(10)

Here, F_c and F_g represent the centrifugal and gravitational forces, respectively. VG is the average-velocity gradient associated with the R1 filament rotation, ξ is the constant associated with the gravitational acceleration distribution (see Table 3), and θ represents the (unknown) inclination angle of the filament relative to the POS.

In Fig. 7, we show the resulting profile for the ratio of forces in the R1 filament (red line) compared to the L1482-S profile (black line) in the California cloud (Álvarez-Gutiérrez et al. 2021). For this comparison, we assumed cos(θ)=1, that is, the filament is not inclined with respect to the POS. For reference, an inclination of 45^° moves the R1 curve upward in this diagram by a factor of ~0.45, closer to the L1482-S curve. In other words, a nonzero POS inclination would increase the role of rotation in the filament. For cos(θ) = 1, the R1 profile indicates that the filament is predominantly affected by gravity rather than rotation, similar to the internal rotation pattern observed in the California L1482-S cloud. The high line-mass profile seen in the R1 filament, combined with the gravitational dominance in this profile, suggests that the filament might primarily be dominated by gravity, even in the presence of rotation. Thus, its most likely fate is to collapse into a state similar to that of R2.

Interestingly, despite the differences in global properties such as total mass and environment for R1 and the California L1482-S filament, both structures exhibit remarkably similar kinematic signatures, including helical patterns and comparable ratios of the centrifugal and gravitational forces. This suggests that such rotation-like features arise in structures filamentary covering a very wide range in line-mass values. Hence, these signatures might reflect more general kinematic processes inherent to filament evolution, potentially making them a common feature in low- and high-mass star-forming environments.

However, the features of potential rotation do not generically always appear in the same line tracer, although we only know of two examples to date (and one is presented here, in N₂H⁺). In contrast to California L1482-S, no rotation is observed in the C¹⁸O gas in the R1 filament. N₂H⁺ and C¹⁸O trace different physical conditions within molecular clouds, reflecting variations in density, temperature, and chemistry that evolve during star formation (e.g., Tafalla et al. 2004; Tanaka et al. 2013). C¹⁸O, a stable CO isotopolog, remains in the gas phase in relatively warm, low-density environments, and it thus traces material prior to significant CO depletion (e.g., Friesen et al. 2010; Punanova et al. 2016; Sabatini et al. 2022). As clouds contract and cool, CO freezes onto dust grains, creating conditions that favor the survival of N₂H⁺ (Caselli et al. 2002b). This chemical progression marks the transition from quiescent to dense, star-forming gas. Hence, rotation detected in different tracers provides clues about the evolutionary stage of a region. While L1482-S (Álvarez-Gutiérrez et al. 2021) and R1 share similar kinematics and low core formation rates, their contrasting linemass profiles suggest that R1 is in a slightly more advanced evolutionary phase.

Fig. 7 Ratio of the centrifugal (Fc) and gravitational (Fg) forces in the R1 filament (red line), where θ represents the inclination angle of the filament relative to the POS. We assumed cos(θ) = 1, corresponding to a filament aligned with the POS. The red shaded area represents the lower limit error of 11 %, derived from the relative abundance uncertainty. For N2H+, the gravitational force dominates, with a distribution similar to the internal rotation observed in the California L1482-S filament (Álvarez-Gutiérrez et al. 2021).

6 Discussion

Despite the morphological similarities between the two main G012 filaments, our analysis revealed differences in their gas kinematics, mass distributions, and core population. In the following discussion, we explore the different star-forming scenarios that may be associated with the R1 and R2 filaments and emphasize their potentially different evolutionary paths.

6.1 N₂H⁺ destruction

In Sect. 3, we compared N₂H⁺ emission with complementary tracers related to ionization, outflows, and dense gas in the G012 protocluster. Specifically, we highlighted a spatial anticorrelation between N₂H⁺, C¹⁸O emission, and the H41α recombination line. Previously, Tobin et al. (2013) estimated the abundances of N₂H⁺ in protostellar systems. In one of these systems (L1157), C¹⁸O traces regions without N₂H⁺. Specifically, the authors found that as the C¹⁸O peaks decrease, N₂H⁺ increases, indicating a clear anticorrelation between the two tracers. This anticorrelation was also detected by Tafalla et al. (2021), who also showed that N₂H⁺ is detected at high column densities (N(H₂) ≥ 10²² cm⁻²; see their Fig. 7). In addition, Yu et al. (2022) analyzed a H ``II``` bubble (S156) associated with the G305 star-forming complex and reported that N<sub>2</sub>H<sup>+</sup> increased far to the edges of S156. This implies that N<sub>2</sub>H<sup>+</sup> might be destroyed by the hot and expanding gas in the bubble. Together with increased CO near H ``II``` regions, the free electrons traced by recombination lines (e.g., H41α), ionized metal emission lines, and continuum emission provide additional methods for mapping these zones in which N₂H⁺ is destroyed (e.g., Yu et al. 2018; Ginsburg et al. 2022b; Armante et al. 2024).

Galván-Madrid et al. (2024) analyzed the 3 mm continuum emission and H41α recombination line data for the full ALMAIMF sample and reported that an increase in emission is mostly seen in evolved regions with OB associations. All of these studies, along with the effect of temperature on the critical density of the N₂H⁺ molecule (ranging from 6.1 × 10⁴ cm⁻³ at 10 K to 2.0 × 10⁴ cm⁻³ at 100 K), suggested a close relation between the decrease in N₂H⁺ and the evolutionary stages of the protocluster. The critical density indicates the conditions under which N₂H⁺ becomes an effective tracer, and they thus link its observed decrease to changes in the protocluster condition. In G012, the 1.3 mm continuum emission is concentrated at the center of the region and north of the R2 filament (see Fig. 1, right panel). Moreover, the 3 mm continuum exhibits high flux intensity distributions that appear as bubbles in the protocluster center (Ginsburg et al. 2022b), similar to the H41α distribution (bottom left panel of Fig. A.3). In addition, the relative N₂H⁺ abundance distribution decreases in regions near H41α bubbles and ionizing emission from hot cores, highlighting the effect of feedback mechanisms in the distribution and abundance of the dense gas traced by N₂H⁺. Because of the high concentrations of H41α and C¹⁸ O in the center of G012, it is likely that the central H ``II``` regions heat the gas in the protocluster and destroy the N₂H⁺ molecule through interactions with free electrons (Vigren et al. 2012) and with the CO molecule (Bergin & Langer 1997; Caselli et al. 2002b; Caselli & Ceccarelli 2012).

Fig. 8 DCN PV diagrams along the R2 filament and toward the central H ``II``` region. The left panel shows a zoom-in of the N2H+ moment 0 map, overlaid with DCN (green) and H41α (yellow) contours. The middle and right panels present DCN PV diagrams extracted from the red and brown boxes shown in the left panel, respectively. The red box highlights the filament extending toward the central H ``II``` bubble in the protocluster, and the brown box traces the R2 filament, including the northern region without N2H+ emission. The red stars indicate the cores detected in DCN (Cunningham et al. 2023). In the middle panel, the dashed black lines mark the center of the H ``II``` bubble (horizontal line) and the systemic velocity estimated from masers (Immer et al. 2014), corresponding to the approximate velocity of the H ``II``` bubble (37 km s−1 ; vertical line).

6.2 DCN accretion evidence in the R2 filament and toward the central OB star cluster

As N₂H⁺ is significantly affected by stellar feedback, particularly by the increase in temperature associated with the central massive OB star cluster, we further investigated the kinematics of the gas in these regions using DCN emission, which is more resilient to temperature enhancements. Although DCN does not trace filamentary structures as extended as those revealed by N₂H⁺, it highlights hotter and denser gas components and reveals filamentary features within the R2 filament and in structures located close to the central OB cluster and the respective H ``II``` bubble.

First, we assessed whether the R2 filament showed evidence of ongoing core accretion processes. Recently, Dell'Ova et al. (2024) derived a dust temperature map of the G012 protocluster, finding values of ~35 K at the top of R2. They also reported a spatial shift between the H2 and N₂H⁺ column density peaks, which they interpreted as the effect of a detected hot core (Armante et al. 2024) that heats the gas and destroys N₂H⁺ at the top of the filament. The DCN emission in this region reveals a clear V-shaped structure. This feature was observed in other filaments undergoing accretion (e.g., Liu et al. 2019; Bonne et al. 2020, 2023; Álvarez-Gutiérrez et al. 2024; Sandoval-Garrido et al. 2025) and is coincident with the hot core location (see Fig. 8, right panel), along with an SiO-elongated feature perpendicular to the filament that is related to an outflow component (see Fig. 3 in Armante et al. 2024, and Fig. A.3). Additional kinematic evidence is found along the filament. Fig. 5 (right panel) shows that the average-velocity diagram reveals a clumping of velocities around ~1 km s⁻¹ in the central and densest region of R2 that is spatially coincident with the integrated-intensity peak. Such clumps in PV space are commonly interpreted as signatures of accretion-driven core formation (e.g., Hacar et al. 2013; Kainulainen et al. 2016; Ladjelate et al. 2020), marking regions where gravity overcomes thermal support (e.g., André et al. 2010; Palau et al. 2013). Thus, the observed velocity clumping in addition to the small velocity spreads in the PV diagrams features along the R2 filament and the similar values in velocity between cores and dense gas further indicates an active role of cores in shaping the filament kinematics.

To place R2 in context, we compared its star-forming properties with those of the ISF OMC-2 filament. From previous catalogs (Megeath et al. 2012, Megeath et al. 2016; Furlan et al. 2016; Kainulainen et al. 2017; Stutz 2018; González Lobos & Stutz 2019), we derived an SFR of 95 M_⊙ Myr⁻¹ and an SFE of 0.008 for protostellar cores in OMC-2. In addition, previous studies estimated a line-mass profile of λ(ω) = 385 M_⊙pc⁻¹ (ω/pc)^0.38 (Stutz 2018). Although the R2 line-mass profile is even steeper, its SFR is lower than that in OMC-2 (see Fig. 6), while both regions exhibit similar efficiencies. This comparison suggests that R2 is already in an evolved phase of filamentary accretion, possibly slightly earlier than OMC-2, but progressing toward a comparable state. This makes R2 an excellent candidate for testing filamentary accretion scenarios with upcoming James Webb Space Telescope (JWST) observations of YSOs in these filaments.

In addition to the accretion evidence identified in the R2 filament, we show in Fig. 8 a set of coherent velocity structures traced by DCN emission extending across the southern portion of R1 (blue box in the left panel) and toward the central region of the protocluster where the OB cluster is located (red box in the left panel). These structures appear as elongated features in DCN (green contours in the left panel) and display organized velocity patters in parallel PV diagrams toward the central region (middle panel). In particular, the parallel PV diagrams reveal multiple velocity gradients that form characteristic V-shaped patterns. While such V-shaped structures are often interpreted as signatures of gas inflow, we note that similar PV morphologies can arise from different kinematic configurations (e.g., rotation or outflows), and therefore do not uniquely trace accretion. One of these structures is detected toward the southern extended portion of the R1 filament (blue box). Here, DCN traces a relatively shallow V-shape that coincides spatially with a similar but weaker velocity pattern seen in the N₂H⁺ emission. Unlike the more prominent structures discussed below, this feature does not appear to be directly associated with a compact core. This may indicate that the V-shaped structure traces gas motions along the filament that are not associated with a compact and localized gravitational potential (e.g., a prestellar or protostellar core, as observed in the R2 filament), and it might be consistent with large-scale gas flows along the southern portion of R1 that might include inflow.

More prominent V-shaped structures are observed closer to the central region of the cloud (red boxes). Specifically, we found two well-defined V-shaped patterns in DCN (see Fig. 8, middle panel), one of which appears to converge toward the position of a massive core detected in Cunningham et al. (2023). The coincidence of the vertex of the V-shaped and the location of the core suggests that this velocity gradient might trace gas flows feeding this object along the filamentary structure (Álvarez-Gutiérrez et al. 2024; Sandoval-Garrido et al. 2025). A second V-shaped velocity structure is detected nearby, whose vertex lies closer to the central region associated with the OB star cluster (represented by a horizontal dashed black line in the figure). The projected position of this feature is consistent with the location of one of the H ``II``` ionized bubbles in G012. While the spatial proximity to the ionized region raises the possibility that this observed gradient might be affected by stellar feedback, several observational characteristics suggest that it is more consistent with inflow motions than with outflowing gas. In particular, the relatively low velocity variations (around 4 km s<sup>−1</sup>) are significantly smaller than those typically observed in outflows, and the emission is traced by DCN, which probes cold and dense gas not usually associated with high-velocity ejected material. Together, these properties favor an interpretation in terms of gas motions directed toward the central H ``II``` region, possibly related to inflow, although contributions from feedback-induced perturbations or other complex velocity configurations cannot be entirely ruled out. We highlight that these types of patterns that connect R1 with the central region are not found for the R2 filament in N₂H⁺ or DCN.

6.3 Potential star-forming scenarios in G012

The formation of stars occurs over diverse timescales, reflecting the intricate processes governing their evolution (Bergin & Tafalla 2007; Evans et al. 2009; Dunham & Vorobyov 2012; Dunham et al. 2014; Motte et al. 2018). These variations have direct implications for the G012 protocluster. In the R1 filament, which contains only two prestellar cores (Armante et al. 2024), the low SFR of ~4.24 M_⊙ Myr⁻¹ and an SFF of ~0.002 suggest that star formation is proceeding slowly despite the high gas density. In contrast, the R2 filament contains prestellar and protostellar cores and shows a significantly higher SFR of ~55.3 M_⊙ Myr⁻¹ and an SFF of ~0.008 (see Sect. 4.3). Although these estimates were derived from subregions (R1 and R2) selected to isolate the kinematically coherent portions of the filaments and not from their full spatial extent, the same trend is also observed for the larger scale structures (see Appendix B): the extended R1 filament still shows lower SFR and SFF values than the extended R2 filament.

The presence of multiple protostellar cores indicates that R2 has so far been more active in forming cores and is a more efficient star-forming region. The lack of protostellar cores and very limited number of prestellar cores in R1, together with the kinematic differences previously reviewed, show the marked contrasts between the two filaments. The differences suggest that R1 is at an earlier evolutionary stage, prior to the filament concentration and active star formation already evident in R2. This implies that, while potentially less evolved than R2, the R1 filament can eventually reach a higher level of core formation as its dense gas reservoir becomes more active. However, the physical processes driving this activity remain uncertain. Similar contrasts in core populations have also been reported in young filaments, such as those in the G351 protocluster (Sandoval-Garrido et al. 2025). However, the extended structure of R1 suggests that its internal evolutionary state might not be uniform. In particular, the lower portion of the filament, located closer to the central region hosting the OB star cluster and its associated H ``II``` bubbles, shows a higher concentration of dense cores and more complex kinematic signatures than the upper part of the R1 filament. The DCN PV diagrams reveal several V-shaped velocity patterns in the extended area at the bottom of R1, some of which converge toward massive cores, while others appear to be spatially connected to the central ionized region. These features might trace localized gas inflows along filamentary structures toward either individual cores or the cluster environment. One possible interpretation is that the expansion of the nearby H ``II``` bubble compresses the surrounding molecular gas, potentially enhancing core formation in the lower portion of R1. Such feedback-driven triggering has been suggested in other star-forming regions where expanding ionized bubbles interact with dense filaments (e.g., Arzoumanian et al. 2022; Neupane et al. 2024). In this context, the higher abundance of dense cores in this bottom and extend part of R1 might reflect a locally accelerated phase of star formation (compared to the upper part, which exhibits potential rotation signatures) superposed on an otherwise less evolved filament.

Additional evidence supporting this evolutionary contrast comes from chemical differentiation in other regions where molecular tracers such as DCN and N₂D⁺ highlight different physical conditions (e.g., Sakai et al. 2022; Cunningham et al. 2023), suggesting that similar mechanisms might be at play in G012. Moreover, the estimated lifetimes of high-mass prestellar cores in Valeille-Manet et al. (2025) are comparable to the velocity gradient timescale in R1 (~0.1 Myr), suggesting that the observed kinematic conditions might set a characteristic timescale for core formation. This supports a causal view where gas kinematics drive core assembly and evolution and do not simply result from it. However, this interpretation likely applies primarily to the filament regions that are not strongly affected by external feedback. In the lower part of the extended R1 structure (the region traced by DCN), which lies closer to the central H ``II``` region, the gas kinematics might also be affected by the expansion of the ionized bubble, which might modify the local conditions for core formation (see above). The contrasting core populations and timescales in R1 and R2 thus highlight how local dynamical conditions regulate the star formation efficiency and rate. Comparisons with other regions reinforce this evolutionary interpretation. As shown in Fig. 9, the PV-diagram features and line-mass profiles of R1 and R2 resemble those of filaments at different stages: R1 shares characteristics with the early-stage quiescent L1482-S region (Álvarez-Gutiérrez et al. 2021), while R2 resembles the more evolved Orion ISF (González Lobos & Stutz 2019). For this comparison to be meaningful, it is important to verify that the datasets probe similar physical and kinematics scales. Although the angular resolutions differ among the above mentioned datasets, the physical resolutions are comparable. The N₂H⁺ observations of González Lobos & Stutz (2019) and Álvarez-Gutiérrez et al. (2021) have spatial resolutions of 17–22" at distances of 420–500 pc, corresponding to 7–11 kau, respectively. Our observations have a beam of 2.3" at 2.4 kpc, corresponding to 6 kau. Therefore, all datasets probe similar physical scales along the filaments. The spectral resolution in our data (0.23 km s⁻¹) is moderately higher than in those studies (0.12-0.16 km s⁻¹), which might smooth very fine velocity substructure but does not affect the identification of global velocity gradients or velocity spread features discussed here.

The coexistence of filaments in distinct evolutionary phases within the same protocluster has been observed in other regions (e.g., Tafalla et al. 2004; André et al. 2010; Schneider et al. 2012; Hacar et al. 2013; Peretto et al. 2014; Zhou et al. 2022; Sandoval-Garrido et al. 2025). For instance, in the Taurus cloud, Hacar et al. (2013) identified coherent filaments exhibiting quiescent and active regions. Similar findings were reported in SDC13 (Peretto et al. 2014) and the Rosette cloud (Schneider et al. 2012), where filaments with different densities, temperatures, and fragmentation levels coexist. Several factors can explain this diversity, including initial density variations, differential gas accretion, and local feedback from protostellar outflows and H ``II```. Denser filaments tend to collapse earlier, while diffuse structures can remain stable or be affected by nearby activity. These mechanisms can lead to asynchronous star formation, where young stellar objects in the same cluster span a wide range of ages (Peretto et al. 2014; Nony et al. 2021). Recognizing the presence of filaments at different stages within a single region like G012 is key to understanding the temporal and spatial complexity of star formation in protoclusters. An intriguing open question is whether the apparent rotation of R1 reflects a genuine large-scale rotational motion with potential inherited angular momentum from the parent molecular cloud or arises from other kinematic configurations. Several alternative scenarios can produce apparent transverse velocity gradients. We list the scenarios below.

• A coherent rotation of the filament would naturally produce a smooth and approximately linear velocity gradient perpendicular to the filament ridge line. This behavior is indeed observed in the R1 filament, where the ridge line-averaged velocity versus radius diagram reveals a well-defined and continuous gradient (Fig. 5, left panel);

• The observed gradient might result from unresolved dense substructures, such as fibers, forming cores, or accretion streamers. Simulations and observations suggest that young filaments can retain rotational signatures from their fragmentation process, with cores often exhibiting rotation axes preferentially aligned relative to the filament orientation (Kong et al. 2019; Hsieh et al. 2021; Lee et al. 2025). Alternatively, sub-filamentary networks such as fibers and streamers that are common in turbulent clouds can serve as conduits for organized flows and potentially impart angular momentum to host filaments (Hacar et al. 2017; Maud et al. 2017; Hacar et al. 2018; Pineda et al. 2023; Olguin et al. 2025). In this case, we would expect a more fragmented velocity structure, with discontinuities or localized deviations along the filament. While such substructures may be present below the current angular resolution, the smoothness of the averaged gradient in the top of R1 suggests that they do not dominate the large-scale kinematics;

• Organized inflows of gas along sub-filamentary networks can also imprint velocity gradients. These flows typically produce asymmetric or localized velocity patterns such as V-shaped structures in parallel PV diagrams. In R1, we analyzed parallel gradients, but found no V-shaped structures in the region that exhibits the proposed rotation feature. We only found a weaker V-shaped pattern that is detected in the bottom extended portion of R1 traced by DCN and N₂H⁺, and near the central H ``II``` region. This might suggest localized inflows toward the central H ``II``` bubble environment rather than along the entire filament (Fig. 8, middle panel);

• External feedback, such as the expansion of nearby H ``II``` regions, might strongly perturb filament kinematics by compressing gas and enhancing velocity dispersions (Zhou et al. 2022; Arzoumanian et al. 2022; Neupane et al. 2024; Zhao et al. 2025). Although these effects can contribute locally, the absence of a comparable transverse gradient in the nearby R2 filament, which lies within the same large-scale environment, suggests that feedback alone cannot explain the kinematic behavior observed in R1.

Of the possible explanations, a coherent rotation of the filament provides the most straightforward interpretation of the observed transverse velocity gradient. While localized inflows or feedback-related perturbations can contribute to the kinematics in specific regions, particularly in the bottom and extended portion of R1 near to the H ``II``` region, the overall smoothness and continuity of the transverse gradient along the filament are most consistent with large-scale rotational motion. Separating these other alternative scenarios requires higherresolution observations capable of resolving sub-parsec structures such as fibers, forming cores, and accretion streamers, and detailed comparisons with simulations of feedback-affected filaments.

Fig. 9 Comparative PV diagrams of California L1482-south (left panel, Álvarez-Gutiérrez et al. 2021), OMC-2 (right panel, González Lobos & Stutz 2019), and the R1 and R2 filaments (center panels). Similar to California, R1 presents a double-helix feature, low presence of cores, and low temperatures (<25 K, see Dell'Ova et al. 2024). Despite this, line-mass profile of both regions (estimated within the range of the red lines on the central panels) shows that R1 is more massive and dense than California L1482-south. Furthermore, R2 shows PV diagram features comparable to the ISF OMC-2 region. Specifically, R2 presents smooth undulations in velocities characterized with uniform velocities. R2 is also characterized by a high presence of cores.

7 Conclusions

G012 is one of the most evolved and massive regions in the ALMA-IMF protocluster sample. Its dense gas emission, traced by the N₂H⁺ molecule, reveals a morphology mainly composed of two large filaments, and additional irregular but similarly dense structures. The center of G012 presents an absence of N₂H⁺, which agrees with the existence of large-scale H ``II``` regions, which destroys the dense gas that forms a disrupted structure in the protocluster center. We focused on the kinematical analysis of the N₂H⁺ emission of the G012 protocluster, including a line-mass distribution analysis, and the kinematic characterization of its two main filamentary structures, which we named R1 and R2. Our main conclusions are summarized below:

1. We found multiple velocity components in G012. About 55% of the N₂H⁺ spectra are well fitted by two velocity components, and the remainder are well described by a single velocity component. We defined the two main velocity structures in the protocluster as the first and second velocity components, FVC and SVC, respectively;

2. The integrated-intensity map revealed two dominant filamentary structures: the R1 and R2 filaments. Despite their similarity in morphology, we found obvious differences in their kinematics;

3. R1 is characterized by a velocity gradient across its spine in the moment 1 map. We review the PV diagram features, which showed a helical structure at the top of the filament that might be associated with filament rotation. We found an average velocity gradient of ~ 10.4 km s⁻¹ pc⁻¹ (or ~0.1 Myr timescale) in this feature;

4. R2 is characterized by smooth, compact velocity undulations (<2 km s⁻¹) and the absence of a large-scale velocity gradient;

5. The mean line-mass profiles of the two filaments are well described by λ(ω) = 5660 M_⊙ pc⁻¹ (ω/pc)^0.30 and λ(ω) = 6943 M_⊙ pc⁻¹ (ω/pc)^0.20 for R1 and R2, respectively. When we compared these distributions with those from other well-studied regions, we found that R1 and R2 stand out as having significantly higher line-mass profiles, which is consistent with higher densities than the structures mentioned above. Based on these line-mass profiles, we calculated the gravitational field, force, and potential for R1 and R2, assuming cylindrical symmetry;

6. Using simple toy models, we estimated in R1 the ratio of the centrifugal to gravitational force (the latter based on the line-mass profile). We found that rotation in the R1 filament is subdominant compared to gravity, but as the radius increases, the role of rotation increases. This is very similar to the behavior reported for the California L1482 filament (Álvarez-Gutiérrez et al. 2021);

7. While the kinematic features appear to be similar between L1482 and R1, the line-mass profile normalization at 1 pc of R1 is almost 30 times higher than that of L1482. Hence, the similarities in the kinematics and the force-ratio profiles might indicate a common angular momentum evolution of filaments across wide ranges in mass and line-mass scales;

8. In terms of the estimated star formation activity, R1 lacks protostellar cores, and the few cores that are present are in the prestellar phase, implying an SFR ~ 4.24 M_⊙ Myr⁻¹. In contrast, R2 hosts a larger population of massive cores than R1 and contains an SiO outflow oriented perpendicular to the filament. Specifically, R2 forms protostellar cores at a rate of SFR ~ 55.3 M_⊙ Myr⁻¹ with a similar efficiency as the Orion ISF;

9. We identified filamentary structures traced by DCN exhibiting signatures of potential accretion toward the center of the protocluster, where the OB stellar cluster is located and N₂H⁺ is largely destroyed, suggesting that these structures might represent the remnants of a former hub-filament system;

10. Using N₂H⁺ (1−0) and DCN (3−2) spectral lines, we estimated updated velocities for 48 previously identified cores. In general, the cores appear to be kinematically coupled to the dense gas traced by N₂H⁺ and DCN from which they are forming;

11. Considering the differences of R1 versus R2 in gas kinematics, core incidence, and estimated SFRs, we proposed that R1 is still rotating and younger than the R2 filament, which has collapsed to a more efficient mode of star formation, likely by shedding its angular momentum.

The results we presented emphasize the role of dense filamentary structures in high-mass star formation, especially when found at different stages of evolution within the same protocluster. The detection of filaments in different evolutionary phases within a single protocluster allowed us to characterize multiple stages of star formation, from the earliest accretion processes to later stages of core collapse, but after feedback effects from H ``II``` regions. By uncovering these structures, we gained valuable insights into the diverse mechanisms governing the gas kinematics, including gas accretion and filament rotation evolution. This ability to observe a snapshot of the different stages in one region presents a unique window into the ongoing processes shaping Milky Way protoclusters.

Acknowledgements

ADS/JAO.ALMA#2017.1.01355.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The project leading to this publication has received support from ORP, which is funded by the European Union's Horizon 2020 research and innovation program under grant agreement No. 101004719 [ORP]. A.S., J.S., N.S.G., R.A.G., N.C.T., and G.B.M. gratefully acknowledge support by the Fondecyt Regular (project code 1220610), and ANID BASAL project FB210003. R.G.M. and J.S. acknowledge support from UNAM-DGAPA-PAPIIT project IN105225. R.A.G. gratefully acknowledges support from ANID Beca Doctorado Nacional 21200897. F.L. and F.M. acknowledge funding from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130) and from the French Agence Nationale de la Recherche (ANR) through the project COSMHIC (ANR-20-CE31-0009). A.K. gratefully acknowledges support from Fondecyt postdoctoral grant. L.B. gratefully acknowledges support by the ANID BASAL project FB210003. N.S.G. gratefully acknowledges support from ANID Beca Doctorado Nacional 21250244. G.B. acknowledges support from the PID2023-146675NB-I00 (MCI-AEI-FEDER, UE) program.

References

Appendix A N₂H⁺ (1−0) line fitting: one and two velocity components

In this appendix, we describe the data preparation steps prior to line fitting, including S/N analysis and input parameter selection. We also characterize the output modeling parameters and associated uncertainties, identifying the velocity structures used in the above analysis.

Appendix A.1 Data preparation

Before fitting, we evaluate the robustness of the N₂H⁺ data cube by inspecting negative bowls and noise levels. The most negative values are concentrated toward the compact central region of G012. To identify reliable spectra, we construct an S/N map, where the noise is estimated pixel-by-pixel from the RMS in emission-free channels between 12-20 km s⁻¹ and 49–61 km s⁻¹. The S/N ratio is defined as the peak intensity per pixel divided by its RMS noise. We adopt an S/N threshold of 12, preserving most of the N₂H⁺ structure while excluding noisy or unresolved spectra unsuitable for fitting. Preliminary moment 0 and 1 maps are constructed using the full line for the integrated-intensity and the isolated component for the mean velocity (see Álvarez-Gutiérrez et al. 2024, for details on isolated component identification). The integrated-intensity map (right panel in Fig. 1) reveals two main N₂H⁺ filaments and a disrupted central morphology in G012. The moment 1 map shows a velocity gradient perpendicular to the R1 filament ridge line (Sect. 4.1) and broad velocity dispersion in the central region, likely associated with multiple velocity components.

Appendix A.2 Choosing input parameters

To model the N₂H⁺ data, we used the specfit fitting tool based on the built-in n2hp_vtau fitter. The specfit task requires the following input values: fittype specifies the model type and requires the number of velocity components and model parameters; guesses are initial model parameter values (only applied to the first pixel fit); limits are lower and upper parameter bounds used across all pixels via the Levenberg-Marquardt algorithm optimizing the χ² function; limited applies or disregards the limits defined above; errmap is the error map, calculated as the standard deviation in the noise channels (see above); signal_cut defines the signal lower limit for the modeling, where S/N > signal_cut; and start_from_point is the initial pixel for fitting.

We test models with one and two velocity components, using four main parameters: the excitation temperature (T_ex), optical depth (τ), velocity centroid (V_c), and velocity dispersion (σ). For the initial guesses and limits, we select values based on the preliminary moment maps described above. The start_from_point is chosen as a pixel with high S/N and a well-defined spectrum. Testing various signal_cut values confirms that pixels below an S/N of 12 yield poor fits or large uncertainties in critical parameters like radial velocities errors. The input parameters for one- and two-component fits are listed in Table A.1.

Appendix A.3 Output parameters and dependencies

PySpecKit produces cubes containing the best-fit parameters (T_ex, τ, V_c, and σ) and their associated uncertainties for each spectrum and velocity component. Using these parameters, we construct a modeled N₂H⁺ cube and separate the first (FVC) and second (SVC) velocity components into individual cubes. To evaluate the dependence of the fitting on the input parameters, we analyze the two most prominent N₂H⁺ filaments (R1 and R2; Fig. 1). We test the start_pixel parameter while keeping all other inputs fixed, performing 716 fits for R1 and 576 for R2, each using a different initial pixel. We find that the derived parameters and uncertainties are largely independent of the initial pixel. However, when fitting multiple velocity components, the starting pixel must correspond to a spectrum where the components are clearly resolved. We also test different parameter limits and find that our adopted ranges are sufficiently broad to avoid biasing the results.

Fig. A.1 N2H+ two velocity component fitting example: The panels, from top to bottom, display the raw data (grey curve), the first velocity component (blue curve) overlaid on the raw data, the second velocity component (green curve) overlaid on the raw data, the total model (black curve) fitted to the raw data, and the residuals of the model (grey dashed curve). The PySpecKit parameters for the first and second velocity components are shown in each respective panel.

Appendix A.4 Output cleaning

With this modeling approach, some spectra may yield poor fits due to noise or an incorrect number of velocity components. Therefore, it is essential to identify the appropriate number of components and assess the quality of the fitted parameters and uncertainties. In particular, we analyze the error distributions of the fitted parameters. Errors in centroid velocity and velocity dispersion are mostly below 0.3 km s⁻¹, with only 5 and 26 pixels exceeding this value, respectively. Since this threshold is conservative relative to the N₂H⁺ velocity resolution, we adopt it as the maximum allowed uncertainty for V_c and σ. We also find 184 fits with σ uncertainties equal to 0 km s⁻¹, indicating that the best-fit values lie at the edge of the allowed parameter range. These spectra likely contain unresolved velocity structure and may require multiple velocity components.

The parameter τ is more difficult to constrain because it can become arbitrarily large. We therefore apply a mask requiring τ/e(τ) > 1, where e(τ) is the uncertainty in τ. Large τ values are typically associated with large uncertainties and flat-top spectral profiles, leading us to exclude ~860 pixels from the model. Most retained τ values are below 25. For T_ex, optically thin pixels (τ < 1) often converge to the upper input limit with uncertainties of 0 K. However, because T_ex has little impact on the spectral profile in this regime, we retain the original fitted values until the column density analysis (Sect. 3.3).

Inspection of the rejected spectra shows that many are associated with multiple velocity components, including cases where the components are difficult to resolve. To address this, we refit the entire cube using two velocity components and evaluate whether the fits improve inside and outside the previous selection criteria. For spectra containing a single velocity component, the two-component model often produces poor fits for either the bluest or reddest component. In these cases, we retain the single-component fit when its uncertainties are low and the fit quality is satisfactory. In addition, poor two-component fits are generally associated with velocity uncertainties larger than 0.3 km s⁻¹ (~15% of pixels), large relative τ uncertainties (~30%), zero-valued parameters (11% of pixels), and high τ values. We additionally impose an upper limit of τ = 100 for each velocity component, removing 52 pixels above this threshold.

In total, we obtain ~15 000 spectra with a well defined model for one (~45%) or two (~55%) velocity components. Two velocity components are spatially distributed in the whole region, without a preferred location, but with a slight concentration in the R2 filament and in regions with large σ values.

Appendix A.5 Model cube merging

We obtain the final modeled cube by merging the fits with single and double velocity component. To carry out this step we follow a similar method than the one performed by Sandoval-Garrido et al. (2025).

To separate the two velocity components, we define a boundary using the midpoint between the upper limit of the first component (mean plus one standard deviation) and the lower limit of the second component (mean minus one standard deviation). Components below 35.6 km s⁻¹ (~3535 pixels) are assigned to the first velocity component (FVC), while those above (~2984 pixels) are assigned to the second velocity component (SVC). This results in two velocity structures containing 7675 pixels in the FVC and 7124 pixels in the SVC. $$\text{limit\hspace{0.17em}}=\frac{\langle {\text{V}}_{\text{c},1}\rangle +\text{std}\left({\text{V}}_{\text{c,1}}\right)+\langle {\text{V}}_{\text{c,2}}\rangle -\text{std}\left({\text{V}}_{\text{c,2}}\right)}{2}\to 35.6{\text{kms}}^{-1}.$$(A.1)

Fig. A.2 Histograms of the velocity distributions for single-(grey) and two-component (blue and green) velocity fits to the spectra (see text). The black dashed line shows the velocity boundary, at 35.6 km s−1, used to distribute one velocity component fits in the merged model cube.

Here V_c,1 and V_c,2 correspond to the velocity center parameter of the first (bluest) and second (reddest) velocity component from the double component model, respectively. This approach ensures the limit reflects a balanced division between the two velocity distributions observed in the histogram in Figure A.2, resulting in a velocity boundary of 35.6 km s⁻¹.

Appendix A.6 Complementary tracers line fitting

For tracers not previously modeled (DCN, H41α, and SiO), we apply a single-Gaussian fit using PySpecKit. The model includes three free parameters: amplitude (Amp), centroid velocity (V_c), and velocity dispersion (σ), with initial guesses derived from the moment maps. Given the simplicity and low computational cost of the model, we do not constrain the parameter ranges. We additionally apply tracer-dependent S/N thresholds through the signal_cut parameter. The S/N is calculated from the peak intensity of each spectrum divided by its RMS noise, estimated from emission-free channels (Table 1).

Fig. A.3 Integrated intensity of the complementary tracers: C18O (J=2−1, upper-left panel; Koley et al. 2025), SiO (J=5−4, upper-right panel; Towner et al. 2024), H41α (bottom left panel; Galván-Madrid et al. 2024), and DCN (J=3−2, bottom right panel; Cunningham et al. 2023). The black contour shows N2H+ integrated intensity at 25 and 100 K km s−1. The scale bar indicates 0.3 pc at the distance of the G012 protocluster. The H41α recombination line and C18O emission are mostly concentrated at the center of G012, where N2H+ is mostly absent and an OB type stars cluster is located (see Fig. 1). In the areas surrounding the main N2H+ filaments, we observe a distribution of lower C18O integrated-intensity emission (~15 K km s−1 ). DCN traces some regions of the main N2H+ structures, and its emission is associated with cores and filaments. Specifically, DCN dominates the top of the R2 filament, central regions of the protocluster, at the edges of R1, and surroundings of the H41α bubbles.

Appendix B Alternative filament detection using FilFinder

Our current filament selection is based on the visual identification of prominent and kinematically well-separated filamentary structures. More specifically, we are interested in the filamentary structures where the most clear and significatively kinematic patterns are. However, the R1 and R2 selected regions could be part of more extended structures or be conected to other protocluster regions. To addres this, and in order to validate the robustness of our filament selection, we use the FilFinder algortihm to independently identify filamentary structures, which we then compare to our original selection.

Fig. B.1 Results from FilFinder filament detection over the FVC (left panel) and SVC (right panel) integrated-intensity maps. In both panels, the yellow contours represent the final masking area where we apply the algorithm (estimated by SNR and spatial extension considerations), green rigdelines represent the skeletons and branches identifying by FilFinder, black rigdelines shows the filament spine identification in Sect. 4.1 based on the integrated-intensity emission in the R1 and R2 regions. We add a zoom in the R1 and R2 regions to compare the differences between FilFinder (green lines) and our previous (black lines) rigdeline identification.

Appendix B.1 Method

We use the FilFinder python package (Koch & Rosolowsky 2016) to identify and characterize filamentary structures in 2D emission maps. The analysis is performed on the integrated-intensity maps of each velocity component, since the structures are primarily separated in velocity space (upper panels of Fig. 2). Filament detection is restricted to regions with S/N > 12. The images are first preprocessed using percentile-based flattening to suppress bright compact emission and enhance extended low-contrast structures. Filament masks are then generated through adaptive thresholding, combined with minimum area criteria and optional hole filling to remove small-scale noise features. The minimum filament sizes are selected based on the spatial resolution of the data. Table B.1 summarizes the adopted mask parameters prior to filament identification. All spatial thresholds (see below) are larger than the beam minor axis (2.1" ≈ 0.024 pc), ensuring that the detected structures are spatially resolved filaments.

We adopt a smoothing scale of 5 pixels (0.05 pc; ~2 beam sizes) to reduce pixel-scale noise while preserving resolved filamentary structures. A minimum area threshold of ~0.05 pc² is used to reject compact or marginally resolved features. The adaptive threshold window is set to ~0.5 pc, matching the characteristic filament scales. We retain only skeletons longer than 0.8 pc to focus on parsec-scale filaments, and prune branches shorter than 0.2 pc to remove minor substructures while preserving the main filamentary morphology. From the final mask, the algorithm computes a medial skeleton representation, reducing each filament to a one-pixel-wide spine (equidistant from the filament edges). This rigdeline estimation differ of our method to estimate the filament spine in Sect. 4.1 where we use the integrated emission in each filament. The resulting skeleton network is analyzed to identify connected components and branches, and is pruned according to length thresholds in order to keep only the main filamentary paths.

Appendix B.2 Variations in physical parameters

In Fig. B.1 we present the results of the FilFinder detection applied to the two velocity component integrated-intensity maps (see upper panels in Fig. 2). The final skeletons identified in the R1 and R2 regions appear more extended than assumed in our baseline definition. Nevertheless, in R1 we recover a ridge line that closely follows the structure adopted in our analysis. In some cases, FilFinder identifies more complex morphologies, particularly in the R2 filament, where most of the elongated structure follows a single main skeleton, but the lower part splits into approximately four possible branches (see right panel in Fig. B.1). Although one could select only the longest path as the main spine of the filament, we consider that it would not be fully reliable to claim that this represents the principal structure. Higher-resolution data will be required to confirm whether these features correspond to dependent or independent substructures near the filament end. For this reason, and to maintain consistency and simplicity in our main kinematic analysis, we restrict our study to the R1 and R2 regions defined at the beginning of this paper (see Fig. 1, right panel), which we consider our main filamentary structures of interest.

To quantify potential variations in our main analysis, we derived key parameters (mass estimation, mass profiles and associated metrics, SFR, and SFF) within the areas identified by FilFinder for the R1 and R2 regions. Within the yellow contours shown in Fig. B.1, which trace the more extended filamentary areas, we derive an H₂ mass of 2450 M_⊙ for the extended R1 structure and 5000 M_⊙ for the extended R2 filament. These values are approximately a factor of two higher than those obtained using our baseline regions alone, emphasizing the significant contribution of these dominant filaments to the total protocluster mass. Overall, about 60% of the protocluster mass is contained within these main filamentary structures.

The cumulative mass profiles of the two extended filaments preserve the approximately linear behavior previously found for R1 and R2. We therefore evaluate how strongly the mass profiles change when considering the extended structures. To align the data, we adopt the FilFinder skeleton output as the ridge-line reference. For R1 we use the single skeleton returned by the algorithm, while for R2 we select the most extended skeleton as the filament ridge line. On average, we find differences of about 30% relative to the metrics derived using only the R1 and R2 regions (see Table 3 for details). The extended filaments show lower M/L values but larger projected radius, which is a direct consequence of the ridge-line alignment. Despite these lower M/L ratios compared to the baseline R1 and R2 regions, the profiles remain significantly higher than those measured in the comparison regions shown in Fig. 6, reinforcing the conclusion that these filaments correspond to highly dense environments and are dominant mass reservoirs in the G012 protocluster.

Fig. B.2 Example of the R1 filament alignment relative to the integrated-intensity moment map. Left panel: integrated-intensity map of the FVC in the R1 region. The red line represent the ridge line estimated based on the integrated-intensity peak along to the filament. The X and Y axis represent the R1 length and width in units of pc. Right panel: integrated-intensity map aligned respect to the ridge line. The X axis represent the projected radius result of the alignment process.

In addition, we quantify the number of cores, their classification (prestellar or protostellar), and their masses to estimate the SFR and SFE within the extended areas. In the extended R1 filament, we now identify a protostellar contribution (in contrast to the original R1 selection, which contained only prestellar cores). In total, we detect nine cores in the R1 extension (two proto-stellar and seven prestellar) and 15 cores in the R2 extension (11 protostellar and 4 prestellar). Following the same timescale assumptions described in Sect. 4.3, we estimate the SFR for each core population and for the total. For the extended R1 filament, we derive a prestellar SFR of 6.25 M_⊙ Myr⁻¹ and a protostellar SFR of 5.40 M_⊙ Myr⁻¹, resulting in a total SFR of 11.65 M_⊙ Myr⁻¹, corresponding to an increase of approximately 60% relative to the baseline region. For the extended R2 filament, we estimate a prestellar SFR of 6.3 M_⊙ Myr⁻¹ and a protostellar SFR of 47.5 M_⊙ Myr⁻¹, yielding a total SFR of 53.9 M_⊙ Myr⁻¹, which represents only a 2% difference compared to the baseline region.

For the SFE, the extended R1 region has a total value of 0.0041 (with a prestellar contribution of 0.0030 and a protostellar contribution of 0.0011) that differ in around 60% with the SFE in the original R1 region, while the extended R2 region has a total SFE of 0.0062 (with a prestellar contribution of 0.0015 and a protostellar contribution of 0.0047), here the SFE differs only by 20% respect to our previous estimation in Sect. 4.3. This indicates that the R2 region forms cores with approximately twice the efficiency of the R1 filament. Despite the larger changes in SFR and SFE in the R1 filament, the extended R2 filament maintains a broader diversity in its core population, with a slight tendency toward more evolved cores, and a higher SFR than R1 (SFR in R2 is around a factor of five higher). The main difference compared to the parameters estimated within the smaller baseline regions is that the SFE gap between the two filaments becomes smaller (SFE in R2 is higher by about a factor of 1.5 than in R1).

The kinematic patterns of rotation and collapse discussed in the main analysis for R1 and R2 regions not change when considering the extended filaments. However, we highlight a kinematic transition between the upper and lower parts of the R1 filament. Specifically, the upper part exhibits the rotation signatures discussed in Sect. 4.1, whereas the lower part shows more dispersed velocity patterns and V-shaped structures that we discuss in Sect. 6. This transition may reflect the influence of the core population in that region or the impact of the feedback produced by the H ``II``` bubbles in the region.

Appendix C Data alignment

We apply a filament alignment procedure throughout this work. For example, in Sect. 4.1 the filaments are aligned relative to the peak integrated intensity to estimate average velocity gradients, while in Sect. 4.2 the alignment is based on the column density distribution. Below we describe the procedure for the first case (Fig. B.2), although the same method can be applied to other datasets:

1. Map rotation: We first rotate the N₂H⁺ maps (e.g., integrated-intensity and velocity centroid) to align the main filaments along the vertical axis (see Fig. 3, upper-left panel).

2. Ridge-line estimation: Using the N₂H⁺ integrated-intensity map, we estimate the filament ridge line with a custom function that traces the peak intensity along the y-axis, defining the filament spine (red line in Fig. B.2). The procedure is applied independently to the R1 and R2 filaments.

3. Ridge-line smoothing: We smoothed the ridge line using a one-dimensional uniform filter from scipy.ndimage with a window size of 1 pixel. This reduces sharp fluctuations and produces a smoother representation of the filament path.

4. Alignment and analysis: The maps are aligned to the smoothed ridge line (Fig. B.2, right panel), transforming the x-axis into projected radius and enabling a more robust analysis of velocity and mass profiles.

Appendix D Technical considerations to estimate relative abundance

Reprojection of the N₂H⁺ column density map can introduce edge artifacts because the regions with sufficient signal-to-noise are highly irregular. To mitigate this effect, we follow the methodology of Sandoval-Garrido et al. (2025) to identify affected pixels: a) We generate a binary mask from the N₂H⁺ column density map, assigning values of 1 to valid pixels and 0 to pixels below the signal-to-noise threshold. b) The mask is reprojected to the H₂ pixel scale of Dell'Ova et al. (2024). Pixels with intermediate values between 0 and 1 identify reprojection artifacts, primarily affecting regions within 1–3 pixels of the image boundaries. c) We exclude all pixels with values < 1 in the reprojected mask, producing a reprojection-artifact-free N₂H⁺ column density map.

To analyze the relative abundance distribution, we construct a histogram using the Freedman-Diaconis method (Sandoval-Garrido et al. 2025), obtaining an optimal bin width of 0.377 × 10⁻¹⁰ and 40 bins. We adopt the mode as the representative abundance value, since previous studies showed it is robust against changes in binning (Sandoval-Garrido et al. 2025). For the full sample, the mode is 0.84 × 10⁻¹⁰. To evaluate the impact of stricter τ uncertainty constraints on the abundance estimate, we apply progressively stronger masks:

1. $\frac{\tau }{e\left(\tau \right)}>1.5$ (~3% pixels removed)

2. $\frac{\tau }{e\left(\tau \right)}>2$ (~7% pixels removed)

3. $\frac{\tau }{e\left(\tau \right)}>3$ (~15% pixels removed)

The mode remains stable between 0.84 × 10⁻¹⁰ and 1.07 × 10⁻¹⁰. To minimize the impact of large τ uncertainties while preserving most of the filamentary structure, we adopt the τ/e(τ) > 2 mask, which removes only ~7% of the pixels.

Fig. D.1 Mass over length profile of the selected region for the mass validation step. Red line represents the resulting profile considering H2 mass estimated in Sect. 4.3, while blue line shows the final profile considering H2 mass from the column density estimated in Dell'Ova et al. (2024).

Appendix E Mass validation

To validate the mass selection used for the line-mass profiles (Sect. 4.2), we compare the mass distributions derived from the H₂ column density map (Dell'Ova et al. 2024) and the N₂H⁺-based mass estimate from Sect. 3.3. We select a well-detected elongated region within the R2 filament, with a length of 0.4 pc and width of 0.03 pc, and align the filament to compare their cumulative mass distributions. Both distributions follow a similar linear trend, differing only slightly in mass (ΔM ~ 300 M_⊙). We therefore derive the M/L profiles following the method described in Sect. 4.2. The resulting profiles show minor differences in slope and the expected offset in normalization (Fig. D.1). We additionally estimate the normalization constants of the apparent volume density (Fq. 5), gravitational potential (Fq. 6), and gravitational acceleration (Fq. 7) using both mass distributions, finding average differences of ~35%.

Given the overall agreement between the two distributions and the improved mass profile obtained for the R1 filament, we adopt the total H₂ mass map (M_tot) derived from the N₂H⁺ data for the analysis in Sect. 4.2. Using M_tot, we obtain an approximately linear cumulative mass profile in the R1 filament over a filament length of 0.6 pc (Fig. F.1).

Fig. F.1 Fxample of cumulative mass profile in the R1 filament. Color curves highlights cumulative mass measurements at different filament radius. Profiles show a linear trend of the cumulative mass along the filament.

Appendix F Core velocities

We estimate the velocities of the G012 cores identified in the catalogs described in Sect. 2 using N₂H⁺ and DCN spectral lines. Of the 38 cores previously detected in DCN (Cunningham et al. 2023), 36 are also detected in N₂H⁺ with S/N > 12. We additionally use the 96 continuum cores identified by Armante et al. (2024), 58 of which lack previous velocity measurements. To derive velocities, we fit the DCN and N₂H⁺ emission with a single-Gaussian model (see Appendix A). The major and minor axes of each core (Tables F.1 and F.2) define the region used to extract average spectra and estimate velocity centroids and mean fitted parameters.

From the Armante et al. (2024) catalog, we identify 24 cores detected in both tracers (DCN ⋂ N₂H⁺ sample), 20 detected only in N₂H⁺, 2 detected only in DCN (DCN new sample), and 12 undetected in both tracers. We define V_DCN and ${\text{V}}_{{\text{N}}_2{\text{H}}^{\text{+}}}$ as the velocity centroids derived from the DCN and N₂H⁺ fits, respectively. Since DCN traces denser, warmer, and more compact gas, we adopt V_DCN when available. To estimate velocities for the N₂H⁺-only sample, we compute the offset ${\text{V}}_{\text{DCN}}-{\text{V}}_{{\text{N}}_2{\text{H}}^{\text{+}}}$ using 60 cores detected in both tracers: 36 from the previous DCN catalog and 24 from the DCN ⋂ N₂H⁺ sample (similar to the method applied in G353.41 by Álvarez-Gutiérrez et al. 2024).

We find that N₂H⁺ velocities are blueshifted by ~0.16 km s⁻¹ relative to DCN, smaller than the spectral resolutions of both tracers (N₂H⁺: 0.23 km s⁻¹ ; DCN: 0.34 km s⁻¹). Using this correction, we increased the sample of cores with velocity estimates by 48 objects: 24 from the DCN ⋂ N₂H⁺ sample, 20 from the N₂H⁺ sample, and 2 from the DCN new sample.

All Tables

All Figures

Fig. 1 Multiwavelength view and molecular gas distribution of the G012 protocluster. Left panel: Spitzer RGB composite figure of the G012 protocluster at 8 μm (red), 4.5 μm (green), and 3.6 μm (blue). The two cyan contours trace the N2H+ integrated-intensity emission at 25 and 100 K km s−1, respectively. Right panel: N2H+ integrated-intensity map with contours at 25 and 100 K kms−1, corresponding to an S/N > 12, shown in black (same levels as in the left panel). The blue boxes highlight the two main filamentary structures, R1 and R2, with estimated lengths of 0.56 pc for both. We show the average spectra of the two regions in the insets. The black curve represents the data within the black boxes. The red markers represent ionizing regions (crosses), radio sources (plusses), and infrared sources (dots) detected in the region by Haschick & Ho (1983). The green contours show continuum emission from the Atacama Large Millimeter/submillimeter Array Band 6 observations (Ginsburg et al. 2022b). The black ellipse in the bottom right corner represents the beam size of the N2H+ data.

In the text

Fig. 2 N2H+(1−0) integrated intensity (upper panels), velocity centroid (middle panels), and velocity dispersion (bottom panels) of the FVC (left panels) and SVC (right panels). The black contours trace the N2H+ integrated-intensity emission at 25 and 100 K km s−1. The black boxes in all panels display the spatial location of the R1 and R2 filaments. The black ellipse in the bottom right corner represents the beam size of the N2H+ data.

In the text

Fig. 3 Position-position and position-velocity diagrams of the N2H+ observations. Upper left panel: G012 integrated intensity of FVC (blue shades) and SVC (green shades). We display the DCN cores from Cunningham et al. (2023) catalog with black pluses. The Armante et al. (2024) core catalog is divided into four categories: cores only detected with N2H+ (red crosses), cores only detected with DCN data (red circles), cores detected with N2H+ and DCN (red triangles), and cores detected neither in N2H+ nor DCN (black triangles). The symbol sizes are proportional to the estimated core mass. The black ellipse in the bottom right corner represents the beam size of N2H+ data. Upper right panel: PV diagram along the y-axis; the velocity axis is subtracted from the N2H+ systemic VLSR of the protocluster (35.5 km s−1). The black arrow illustrates a slope of 5 pc (km s−1 )−1, which corresponds to a timescale of ~0.2 Myr, that is, the timescales that approximately correspond to some of the extended structures in this PV diagram. At the top of R1, lies a wrapped (or double-helix) type velocity field with spreads of more than ~3 km s−1. On the other hand, R2 appears very compact in velocity along its extent, with small-scale spatial wiggles. Bottom left panel: PV diagram in the perpendicular direction compared to the upper right panel. This shows the emergence, albeit somewhat hidden in the overall velocity field, of an approximately uniform and extended gradient in R1 that is apparent near ΔX ~ 0.60 to 0.75 pc, ΔV ~ −2 to 0 km s−1. In contrast, R2 appears as the compact blob of cores (triangles and x-symbols) on the right-hand-side of the panel, characterized by the absence of an obvious gradient in position and velocity. Bottom right panel: zoomed PV diagrams of the main filaments R1 (left panel) and R2 (right panel), enclosed in the black boxes of the upper left panel. R1 presents a wrapping double-helix signature that is most obvious toward the top of the diagram and which is dominated by the FVC velocities. R2 exhibits comparatively compact velocity variations (ΔVmax~1.5 km s−1) along the filament and contains a high number of massive cores (1–3 M⊙).

In the text

Fig. 4 Left panel: relative abundance map in G012. The black contour highlights the mask of $\frac{\tau }{e\left(\tau \right)}>2$ applied to the N2H+ column density (see Sect. 3.4); most of the pixels removed by this mask do not affect the main filaments significantly. The areas lacking data (top of R1 and bottom of R2) are due to the H2 column density map lack of coverage in the filaments. The representative relative abundance value is 0.93 × 10−10. The black circle in the bottom left corner represents the beam size of the H2 map. Right panel: relative abundance histogram of the values inside the black contour in the left panel. The dashed red line represents the mode of the distribution. In the upper-right corner, the values of the mode and standard deviation are shown, which we consider as the representative value and error of the sample.

In the text

Fig. 5 Average velocity gradients perpendicular to R1 (left panel) and R2 (right panel) filaments. In both diagrams, the ΔV and Δr axes have the same range, allowing for a direct comparison of the slopes between the panels. In R1 the red line represents the linear fit weighted by the integrated-intensity points of the velocity gradient (VG). The associated VG timescale (τ) estimate is displayed in the upper left corner. In contrast, R2 lacks a clear velocity gradient structure as in R1 and is instead comparatively compact in its velocity distribution. In addition, we identified two different velocity distributions that spatially correspond to different regions in the R2 filament: 1) the most compact and central structure (compact feature) is related to the densest region of the filament; and 2) the most elongated feature is spatially related to the eastern edge of the R2 filament, characterized by elongated SiO emission. In the diagram, the lower intensity and scattered points represent regions surrounding the filament.

In the text

Fig. 6 Line-mass profile of the R1 and R2 filaments (red lines). The red shaded area represents the uncertainty in the two estimates, derived from the relative abundance uncertainty. This percentage therefore provides a lower limit to the error of 11% in both distributions. We include profiles of the Orion ISF (Stutz & Gould 2016), the ONC (Stutz 2018), the California L1482 region (Álvarez-Gutiérrez et al. 2021), and the G351.77 protocluster (Reyes-Reyes et al. 2024) with black lines. The main filaments of G012 reveal denser line-mass distributions than those found in other star-forming regions.

In the text

Fig. 7 Ratio of the centrifugal (Fc) and gravitational (Fg) forces in the R1 filament (red line), where θ represents the inclination angle of the filament relative to the POS. We assumed cos(θ) = 1, corresponding to a filament aligned with the POS. The red shaded area represents the lower limit error of 11 %, derived from the relative abundance uncertainty. For N2H+, the gravitational force dominates, with a distribution similar to the internal rotation observed in the California L1482-S filament (Álvarez-Gutiérrez et al. 2021).

In the text

Fig. 8 DCN PV diagrams along the R2 filament and toward the central H ``II``` region. The left panel shows a zoom-in of the N2H+ moment 0 map, overlaid with DCN (green) and H41α (yellow) contours. The middle and right panels present DCN PV diagrams extracted from the red and brown boxes shown in the left panel, respectively. The red box highlights the filament extending toward the central H ``II``` bubble in the protocluster, and the brown box traces the R2 filament, including the northern region without N2H+ emission. The red stars indicate the cores detected in DCN (Cunningham et al. 2023). In the middle panel, the dashed black lines mark the center of the H ``II``` bubble (horizontal line) and the systemic velocity estimated from masers (Immer et al. 2014), corresponding to the approximate velocity of the H ``II``` bubble (37 km s−1 ; vertical line).

In the text

Fig. 9 Comparative PV diagrams of California L1482-south (left panel, Álvarez-Gutiérrez et al. 2021), OMC-2 (right panel, González Lobos & Stutz 2019), and the R1 and R2 filaments (center panels). Similar to California, R1 presents a double-helix feature, low presence of cores, and low temperatures (<25 K, see Dell'Ova et al. 2024). Despite this, line-mass profile of both regions (estimated within the range of the red lines on the central panels) shows that R1 is more massive and dense than California L1482-south. Furthermore, R2 shows PV diagram features comparable to the ISF OMC-2 region. Specifically, R2 presents smooth undulations in velocities characterized with uniform velocities. R2 is also characterized by a high presence of cores.

In the text

Fig. A.1 N2H+ two velocity component fitting example: The panels, from top to bottom, display the raw data (grey curve), the first velocity component (blue curve) overlaid on the raw data, the second velocity component (green curve) overlaid on the raw data, the total model (black curve) fitted to the raw data, and the residuals of the model (grey dashed curve). The PySpecKit parameters for the first and second velocity components are shown in each respective panel.

In the text

	Fig. A.2 Histograms of the velocity distributions for single-(grey) and two-component (blue and green) velocity fits to the spectra (see text). The black dashed line shows the velocity boundary, at 35.6 km s−1, used to distribute one velocity component fits in the merged model cube.
In the text

Fig. A.3 Integrated intensity of the complementary tracers: C18O (J=2−1, upper-left panel; Koley et al. 2025), SiO (J=5−4, upper-right panel; Towner et al. 2024), H41α (bottom left panel; Galván-Madrid et al. 2024), and DCN (J=3−2, bottom right panel; Cunningham et al. 2023). The black contour shows N2H+ integrated intensity at 25 and 100 K km s−1. The scale bar indicates 0.3 pc at the distance of the G012 protocluster. The H41α recombination line and C18O emission are mostly concentrated at the center of G012, where N2H+ is mostly absent and an OB type stars cluster is located (see Fig. 1). In the areas surrounding the main N2H+ filaments, we observe a distribution of lower C18O integrated-intensity emission (~15 K km s−1 ). DCN traces some regions of the main N2H+ structures, and its emission is associated with cores and filaments. Specifically, DCN dominates the top of the R2 filament, central regions of the protocluster, at the edges of R1, and surroundings of the H41α bubbles.

In the text

Fig. B.1 Results from FilFinder filament detection over the FVC (left panel) and SVC (right panel) integrated-intensity maps. In both panels, the yellow contours represent the final masking area where we apply the algorithm (estimated by SNR and spatial extension considerations), green rigdelines represent the skeletons and branches identifying by FilFinder, black rigdelines shows the filament spine identification in Sect. 4.1 based on the integrated-intensity emission in the R1 and R2 regions. We add a zoom in the R1 and R2 regions to compare the differences between FilFinder (green lines) and our previous (black lines) rigdeline identification.

In the text

Fig. B.2 Example of the R1 filament alignment relative to the integrated-intensity moment map. Left panel: integrated-intensity map of the FVC in the R1 region. The red line represent the ridge line estimated based on the integrated-intensity peak along to the filament. The X and Y axis represent the R1 length and width in units of pc. Right panel: integrated-intensity map aligned respect to the ridge line. The X axis represent the projected radius result of the alignment process.

In the text

	Fig. D.1 Mass over length profile of the selected region for the mass validation step. Red line represents the resulting profile considering H2 mass estimated in Sect. 4.3, while blue line shows the final profile considering H2 mass from the column density estimated in Dell'Ova et al. (2024).
In the text

	Fig. F.1 Fxample of cumulative mass profile in the R1 filament. Color curves highlights cumulative mass measurements at different filament radius. Profiles show a linear trend of the cumulative mass along the filament.
In the text