© The Authors 2025

1 Introduction

The study of the formation and evolution of dense filamentary structures in the interstellar medium is critical for understanding star formation (André et al. 2014). All molecular clouds contain filamentary structures, as demonstrated in the Herschel Gould Belt Survey (André et al. 2010). Furthermore, the majority of gravitationally unstable cores, known as pre-stellar cores (Ward-Thompson et al. 2007), which are on the brink of star formation, are located within these filaments (André et al. 2010, Könyves et al. 2015). This suggests that dense cores form within the filamentary structures of molecular clouds. Magnetic fields are widely recognized as playing a key role in the formation of these filaments, the cores within those filaments, and the eventual stars within them (Pattle et al. 2023). The relative contributions of magnetic fields, turbulence, and gravity to filamentary structure formation remain poorly understood (Li et al. 2014, Crutcher 2012). Gravity and turbulence influence both the structure and strength of the magnetic field (Hennebelle & Inutsuka 2019), which in turn shapes the gas dynamics within filaments. Magnetic fields are thought to help organize material into filamentary structures, but their precise role in this process remains unclear; this is especially true for quiescent filaments, which may not always channel material toward denser regions.

Dust polarization observations, notably from the Planck satellite, have revealed ordered magnetic field structures in molecular clouds, particularly in the Gould Belt clouds of the solar neighborhood (Planck Collaboration Int. XXXV 2016). These observations show that in regions of low column density, the gas structure tends to align parallel to the magnetic field, and vice versa This alignment transition, observed at visual extinctions of A_V ~2.7–3.5 mag, underscores the role of magnetic fields in shaping filamentary structures in star-forming regions (Planck Collaboration Int. XXXV 2016, Soler et al. 2017). Recent High-resolution Airborne Wideband Camera-Plus (HAWC+) observations of the Serpens South cloud reveal a second transition in the alignment of gas structures with magnetic fields at A_V ~ 21 mag, where the alignment shifts back to being parallel (Pillai et al. 2020).

In this study, we focused on the Snake filament, also known as Barnard 72, employing both polarization data and molecular line observations to investigate the role of magnetic fields in shaping its filamentary structure. This cloud has received limited attention in the literature; to our knowledge, the only reference to it is found in Nielbock et al. (2012), where it is mentioned but not studied in detail. Molecular line observations were performed using the 30 m IRAM (Institut de Radioastronomie Millimétrique) telescope; we observed the J = 1 → 0 transitions of C¹⁸O and ¹³CO. Additionally, we obtained polarimetric observations from the 1.6 m and 60 cm telescopes at the Observatório do Pico dos Dias (OPD) that span optical and near-infrared (NIR) wavelengths. Through the analysis of polarization properties and gas behavior within the Snake filament, we aim to understand how magnetic fields influence its structure and dynamics, thereby advancing our knowledge of star formation in filamentary regions.

This paper is structured as follows: In Sect. 2 we detail the observational data obtained from multiple telescopes, including submillimeter polarimetry from Planck, molecular line observations from the IRAM 30m telescope, and optical and NIR polarimetric data from the OPD Laboratório Nacional de Astrofísica (LNA) telescopes. The methods for data reduction and calibration are also described. In Sect. 3, we present the analysis and results of data processing, focusing on the filament’s polarimetric and kinematic properties, and examining the relationship between the magnetic field geometry and the filament’s structure. Section 4 investigates the stability of the filament and the polarization efficiency. Finally, in Sect. 5 we summarize the key findings of this work.

2 Data

2.1 Molecular line data

The Snake filament was observed using the 30 m telescope at Pico Veleta (Spain) from IRAM, under project 038-11, in August of 2011. The ground-state transition of C¹⁸O (1-0) and ¹³CO (1-0) at a rest-frame frequency of 110.201354 MHz and 109.782173 MHz was observed, respectively. The observations were made using EMIR E090 receiver, the VErsatile SPectrometer Assembly (VESPA), with a resolution of 20 kHz. The data were calibrated and reduced using CLASS software, GILDAS¹ software. The mean rms noise level in T_mb scale is 0.3 K with 0.05 km s⁻¹ channels.

2.2 Starlight polarization data

The polarimetric observations were carried out using the 1.6 m and the IAG 60 cm telescopes at the OPD² (Brazil) during several missions, 2007, 2008, and 2024 (IAG 60 cm), for the optical data, conducted in the R band (6474 Å) and, 2018 and 2019, for the NIR data, conducted in the H band. The data were obtained using IAGPOL, an imaging polarimeter specifically adapted for polarimetric measurements, consisting of a rotating retarder λ/2-waveplate, a calcite Savart prism, and a filter wheel. For further details on the setup of the polarimeter, refer to Magalhaes et al. (1996). The optical polarimetry was obtained for 20 fields, with a field of view of about 10′ × 10′ each, covering the entire region of the Snake filament. The NIR data were obtained for 11 fields, with a field of view of about 4′ × 4′ each. For each optical field, a set of eight positions of the λ/2-waveplate, positioned at intervals of ${22}_{\cdot }^{\circ }5$, was obtained with exposure of 100 to 120 s at each position. In 2024 we used the IAG 60 cm telescope to obtain additional optical data. For the NIR, 60 images of 10 s each were acquired for each position of the λ/2-waveplate, while the telescope was dithering to remove the thermal background signal, adding up to a total exposure of 600 s per λ/2-waveplate position. The same procedure was repeated for each of the 11 fields of view.

Image processing was performed using the reduction pipeline SOLVEPOL (Ramírez et al. 2017). After standard procedures of image reduction, the pipeline calculates the differential aperture photometry of the two components generated by the calcite prism, at each position of the wave plate. Observations of standard stars of known polarization angle were used to determine the correction of the polarization angle to the equatorial coordinate system, according to which the position angle is measured east from the celestial north pole. Unpolarized standard stars were also observed in order to determine the instrumental polarization, which turned out to be smaller than 0.1%. The transformation of pixel to sky coordinates was performed by matching the pixel of one of the two components of each star found in the reduced images with celestial coordinates of their counterparts in the Gaia DR3 catalog (Gaia Collaboration 2023). The final rms of the transformation solutions was about $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$. In this study, we assumed that the angle of optical and NIR starlight polarization, denoted as ϕ_star, is identical to the orientation of the magnetic field, ψ.

2.3 Polarized dust emission from Planck

We used the 353 GHz Planck³ polarization data at a resolution of 4′.8. From the Stokes Q and U, we studied the polarization angles of dust emission in the sky region containing the Snake filament. To improve the signal-to-noise ratio (S/N) of the extended emission, we applied a Gaussian convolution to the Planck beam, achieving 10′ resolution maps. We derived the polarization angle of dust emission⁴, ϕ = 0.5 tan⁻¹(U/Q). The polarized intensity is computed as $$P=\sqrt{Q^2+U^2}.$$(1)

The Stokes parameters obtained from the Planck database are measured eastward from the north Galactic pole. These measurements were then converted to the equatorial coordinate system, with angles measured eastward from the north celestial pole. We applied the method suggested by Corradi et al. (1998). Using the equation $$\text{\hspace{0.17em}offset angle\hspace{0.17em}}=\arctan \left(\frac{\cos (l-{32.9}^{\circ })}{\cos b\cot {62.9}^{\circ }-\sin b\sin (l-{32.9}^{\circ })}\right),$$(2)

and considering the Snake filament has galactic coordinate l ~ 1.8° and b ~ 6.9°, we calculated the offset of 56°. Since the area covered by the Snake filament is relatively small. The error introduced by a common offset is smaller than 0.25° in the borders of the studied area. We applied the same offset all over the field, and 56° was subtracted from ϕ.

The submillimeter polarization is assumed to be perpendicular to the magnetic field orientation, ψ. We rotated the angles by 90° to obtain the corresponding plane-of-sky magnetic field orientations, ψ_submm. Therefore, throughout this work, the polarization data from optical, NIR, and submillimeter wavelengths consistently show that the polarization vectors align with the magnetic field direction in their respective regions.

2.4 Data handling and cross-match with the StarHorse catalog

For our analysis throughout this paper, we utilized data that correspond to entries in the StarHorse catalog (Anders et al. 2022). The StarHorse catalog provides stellar parameters, for example distances and extinctions, for 362 million stars, by combining Gaia Early Data Release 3 data (Gaia Collaboration 2016, 2021) with photometric surveys such as Pan-STARRS1 (Chambers et al. 2016), SkyMapper (Onken et al. 2019), 2MASS (Two Micron All Sky Survey; Cutri et al. 2003), and AllWISE (All-Sky Infrared Survey Explorer; Cutri et al. 2013). Its improved precision and broad wavelength coverage enable distance accuracies of approximately 3% for stars with a magnitude of 14 and 15% for magnitudes of 17 (Anders et al. 2022). As discussed in Sect. 3, the Snake Nebula appears to be the dominant structure toward the studied line of sight, up to a distance of 1.2 kpc; in addition, at that distance, the line of sight is at a height of ~150 pc above the Galactic plane, avoiding the densest regions of the Galactic disk, as indicated by the lack of significant cloud material beyond 1 kpc in Vergely et al. (2022). Because of that, we restricted our analysis to stars located within 2 kpc. This selection minimizes uncertainties associated with both distance estimates and polarization measurements, which tend to increase for more distant and fainter stars. Consequently, all analyses and visualizations presented herein focus solely on stars within a distance of 2 kpc. In the NIR dataset, we identified 1244 stars, with 491 of these entries corresponding to the StarHorse catalog. Notably, 57 stars are situated within 2 kpc of the Sun. Conversely, the optical dataset includes 19627 stars, of which 13252 match entries in the StarHorse catalog. Among these, 2037 stars are located within 2 kpc of the Sun.

2.5 Archive data

We employed archival data from the Gould Belt Survey, collected with the Herschel Space Observatory, to derive the dust temperature map (André et al. 2010; Roy et al. 2014). A visual extinction map was constructed from a deep NIR imaging survey, combining data from the New Technology Telescope (NTT), the Very Large Telescope (VLT), the 3.5 meter telescope at the Centro Astronómico Hispano Alemán (CAHA), and the 2MASS survey (for more details see Román-Zúñiga et al. 2010).

3 Analysis and results

3.1 Distance of the cloud and trends with stellar distance

We used the 3D dust extinction map provided by Edenhofer et al. (2024) to obtain the dust distribution within the Snake filament region. This map, which extends up to 1.25 kpc from the Sun with a 14′ angular resolution, is based on data from 54 million stars as analyzed by Zhang et al. (2023). They derived the stars’ atmospheric parameters, distances, and extinctions by forward-modeling the low-resolution Gaia BP/RP spectra (Carrasco et al. 2021). We employed the publicly accessible version of the map⁵. Figure 1 shows the dust extinction distribution as a function of distance along the line of sight toward the Snake filament. To obtain a representative profile, we computed the average extinction over a grid of Galactic coordinates covering the main body of the filament (l = 1.60–2.09, b = 6.58–7.60). We scaled the map values by 2.8 to convert them to A_V in magnitudes based on the recommendation of Edenhofer et al. (2024). This averaging reveals a well-defined extinction peak at a distance of 154 pc, which we interpret as the dominant contribution from the Snake filament (Fig. 1). This distance agrees with the estimated distance of the entire Pipe Nebula (145 ± 16 pc, Alves & Franco 2007; 163 ± 5 pc, Dzib et al. 2018).

Figure 2 (upper panel) illustrates the trend of increasing visual extinction as a function of distance, characterized by a gradual slope for stars with distances up to 2 kpc. The plot also includes the mean visual extinction (A_V) calculated in bins of 100 pc. This line offers a smoothed representation of the overall trend of increasing extinction with distance, helping to highlight the general pattern across the dataset. We also show the ±1σ variation from the mean, providing an uncertainty in the extinction measurements over distance. The optical and NIR datasets show similar trends, while fewer data points are in the NIR dataset. The bottom panel of Fig. 2 presents the polarization percentage as a function of distance using the same dataset. The mean polarization percentage, calculated in 100 pc bins, reveals an initial increase in polarization, followed by a relatively steady trend extending up to 2 kpc. This steady trend indicates the absence of any dense clouds located beyond the Snake filament, up to the 2 kpc threshold. For Fig. 2, we used S/N > 1 for the optical data to include low-polarization stars, which are likely located in the foreground of the Snake filament. For NIR data, we applied S/N > 3, as the observing process for NIR is inherently noisier, resulting in lower S/N values. For the remaining analysis, we used a more stringent criterion of S/N ≥ 5 for the optical data and S/N ≥ 3 for the NIR data, to ensure data reliability while retaining sufficient sources.

Fig. 1 2D dust extinction map illustrating the dust distribution as a function of distance up to 1.25 kpc from the Sun (Edenhofer et al. 2024). The peak at 154 pc indicates a distance corresponding to the Snake region.

Fig. 2 Top panel: visual extinction (AV) as a function of the distance of stars within 2 kpc of the Sun. The extinction and distance values are derived from the StarHorse catalog. Lower panel: polarization percentage as a function of the distance. A vertical dashed line at 154 ± 15 pc marks the distance of the Snake filament. The error bars represent the uncertainties in the extinction and polarization measurements, scaled by 50%. The purple line represents the mean value in bins of 100 pc, with the light blue shaded area indicating the ±1σ uncertainty.

Fig. 3 Dust extinction map of the Snake region in color scale at a spatial resolution of 20″ from Román-Zúñiga et al. (2010). All segments represent the B-field orientation. Blue and red segments are inferred from the optical and NIR data, respectively, with S/N≥3. The green segments show the magnetic field orientation from Planck data at the resolution of 10′. The vectors are scaled to the same size to enhance visualization clarity. The cyan contours represent the C18O (1-0) integrated intensity at levels of 0.8 K km/s and 1.6 K km/s.

3.2 Spatial distribution of the magnetic field

Figure 3 shows the dust extinction map of the Snake filament region at a spatial resolution of 20″ (Román-Zúñiga et al. 2010). Polarization vectors in optical, NIR, and submillimeter are plotted over the extinction map.

The submillimeter vectors are rotated by 90° to show the B-field orientation. The NIR data probe the denser part of the Snake filament, whereas the optical data trace the lower visual extinction of the observed field.

Figure 4 reveals two prominent peaks in the field lines orientation angles. These two peaks are visible in the optical, NIR, and submillimeter data. Specifically, the optical and NIR data exhibit their first peaks within the range of approximately 85° to 102. The optical data also show a second prominent peak at around 180°. In contrast, the submillimeter data reveal two distinct peaks at 66° and 169°, with the second peak of the NIR data aligning closely with the second peak observed in the submillimeter data.

In contrast, the infrared data have a random distribution in Fig. 4. The overall polarization distribution across the entire region does not appear to follow any specific pattern. However, it is important to note that the number of reliable NIR polarization detections is limited compared to the optical data, which reduces the statistical significance of any observed trends in the angle distribution. Despite this limitation, a peak in the distribution can be observed around 80° to 100°.

Fig. 4 Histogram of the distribution of magnetic polarization angles for the three datasets (optical, NIR, and submillimeter). Optical data peaks are highlighted by dashed blue lines and the submillimeter dataset distribution peaks by dashed green lines.

3.3 The relative orientation of the filament with the magnetic field

To further characterize the Snake’s filamentary structure and its relationship with polarization data, we used the Python package radfill (Zucker & Chen 2018). The spine of the filament is defined by radfil using the fil_finder package (Koch & Rosolowsky 2015) based on the Av map and mask map provided. We extracted evenly spaced cuts perpendicular to the spine at intervals corresponding to every 7 pixels along its length. Consequently, the center of the resulting profiles aligns with the peak value of Av. Using this spine, we determined a filament length of 2.3 pc. Figure 5 illustrates the spine and the perpendicular cuts. We selected a 0.04 degree distance from the spine to make sure we covered the whole snake cloud in our sample (see Fig. 5). The plot also includes polarization data from optical and infrared wavelengths. The projections of the polarization vectors onto the spine was done by finding each polarization vector’s closest point on the spine and shifting it accordingly.

Our goal was to analyze how the polarization angles change as we move along the filament’s spine. The behavior of polarization vectors along the spine provides critical insights into the magnetic field geometry. Specifically, analyzing the relationship between the polarization vectors and the filament’s curvature allows us to evaluate whether the vectors are predominantly aligned with the filament (parallel) or if there are regions where they deviate significantly (perpendicular).

To quantify the alignment between the magnetic field vectors and the filament orientation, we calculated the projected Rayleigh statistic (PRS; Jow et al. 2018). This method provides a statistical measure of whether the polarization angles are preferentially aligned parallel or perpendicular to the filament’s local orientation (tangent to the spine). The filament spine is divided into 44 positions along its length. For each of these positions, we computed the PRS value using all polarization vectors projected onto that spine location (cyan vectors in Fig. 5), determined based on the closest spatial distance. Following Jow et al. (2018), the PRS was defined as $$PRS=\frac{\sum _i^n\cos {\theta }_i}{\sqrt{n/2}},$$(3)

where θ is twice the difference between the filament spine angle and the polarization angle. PRS >> 0 indicates strong parallel alignment, and PRS << 0 indicates strong perpendicular alignment. The uncertainty associated with the PRS value is from Jow et al. (2018): $${\sigma }^2=\frac{2\sum _i^n{(\cos {\theta }_i)}^2-{\mathrm{P}\mathrm{R}\mathrm{S}}^2}{n}.$$(4)

Figure 6 shows the PRS values as a function of right ascension (RA), with error bars representing the uncertainties at each position along the spine. At both ends of the Snake filament, the polarization vectors are predominantly aligned with the spine, as indicated by PRS values greater than zero, above the dotted line. In contrast, near the central, denser region of the filament, the polarization vectors tend to become more perpendicular to the spine, reflected in the negative PRS values. The choice of maximum radius from the filament spine affects the PRS values, as it determines which polarization vectors are included in the analysis. To assess this effect, we performed the PRS calculation using a smaller radius of 0.02° (compared to a distance of 0.04° from the spine). The comparison reveals that approximately 9% of data show a change in the sign of their PRS value (from positive to negative or vice versa). This suggests that while the exact numerical values of PRS may vary slightly with radius, the overall conclusion regarding the alignment of the magnetic field with the filament spine remains robust within this tested range. We are committed to using a radius of 0.04° in our analysis, as this selection ensures that the entire Snake filament is encompassed, as defined by the boundaries of the dust extinction map.

Fig. 5 Dust extinction map illustrating the masked filamentary structure. The spine of the filament is outlined by a prominent red curve, and the thin red lines indicate the locations of perpendicular cuts. Blue circles represent the peak pixel intensity along each cut. The orange vectors represent optical (S/N≥5) and infrared (S/N≥3) polarization segments. The cyan vectors are projections of all the orange polarization vectors onto the spine. Each polarization vector was projected by finding its closest point on the spine and shifting it accordingly.

Fig. 6 PRS values as a function of RA. The PRS is calculated at each position along the filament spine using the group of polarization angles projected onto that location. Each point represents the PRS for a given spine position. The dotted line indicates PRS = 0, while the error bars represent the uncertainties.

3.4 Spectral line fitting

CO isotopologues are widely utilized to probe gas in molecular clouds at varying densities. The ¹³CO (1-0) isotopologue is typically associated with tracing the diffuse outer regions of molecular clouds. In contrast, the rarer C¹⁸O (1-0) line, with a critical density of ~10³ cm⁻³, serves as an excellent tracer for regions with higher column and volume densities (relative to ¹³CO (1-0)), effectively avoiding saturation. In this section we compare the polarization data with the gas kinematics using the complete set of NIR observations rather than limiting the analysis to data cross-matched with the StarHorse catalog (see Sect. 2.4). We established that no significant structures exist behind the Snake filament; therefore, all observed NIR data points are considered valid for this analysis.

Figure 7 displays the ¹³CO (1-0) and C¹⁸O (1-0) channel maps in 0.13 and 0.07 km s⁻¹ steps, respectively. For both lines, the data presented have an S/N>3. Polarization vectors for optical (white) and NIR (cyan) wavelengths are overplotted on each panel. The ¹³CO (1-0) emission at velocities lower than 4.4 km s⁻¹ is confined to the east region of the filament, while higher velocities show more concentrated emission at the center of the filamentary structure. At the filament’s center, velocities range from 4.53 km s⁻¹ to 5.18 km s⁻¹, transitioning from west to east (see the second row). While the C¹⁸O (1-0) emission is only concentrated at the center of the filament, it appears at lower velocities on the western side of the map and shifts to higher velocities toward the east side.

The kinematic properties of the observed C¹⁸O (1-0) molecular line were derived by fitting Gaussian profiles to each spectrum in the dataset using the Python package pyspeckit (Ginsburg & Mirocha 2011). A threshold of S/N > 5 was applied, ensuring the fits were constrained to regions with sufficient signal. A single-Gaussian fit was performed across the unmasked pixels with initial parameter guesses refined iteratively using neighboring pixel values. Fits with residuals exceeding twice the RMS or velocity dispersions (σ_v) broader than 0.35 km s⁻¹ were masked and excluded. The resulting C¹⁸O (1-0) model cube showed only 0.1% of the pixels had a bad fit, demonstrating the robustness of the one-Gaussian fitting approach. The top panel of Fig. 8 displays the centroid velocity map of the C¹⁸O (1-0) line. The C¹⁸O (1-0) line exhibits a narrow linewidth, ranging from ~0.1 to 0.2 km s⁻¹, indicative of quiescent C¹⁸O (1-0) gas dynamics in this region. Given that the sound speed in the region is c_c = 0.2 km s⁻¹ (for T = 15 K, from the Herschel map), these velocity dispersions suggest that the motions are transonic (see Appendix A for dispersion velocity map and sound speed calculation).

In contrast, the ¹³CO (1-0) line exhibits a more complex kinematics, particularly in the central regions of the Snake filament, where two velocity components are observed. Similar to C¹⁸O (1-0), a single-Gaussian fit was first applied to pixels with S/N > 4, followed by a dual-Gaussian fit for regions with significant residuals or broad line widths. Through visual inspection, we find that lines broader than 0.7 km s⁻¹ exhibit profiles consistent with two velocity components. Based on this observation, we established a masking criterion that excludes pixels where the residuals exceed twice the RMS and velocity dispersions broader than 0.7 km s⁻¹ in the single-Gaussian fit. Accurately fitting the ¹³CO (1-0) line with two Gaussian components in the central region of the Snake filament is challenging due to the line’s high opacity, which exceeds a value of 1 (see Sect. 3.5). The middle panel of Fig. 8 presents the centroid velocity map of the ¹³CO (1-0) line, derived using a single-Gaussian fitting. The central velocities (V_LSR) of the ¹³CO (1-0) line range approximately from ~4.7km s⁻¹ to ~5.2 km s⁻¹, increasing from west to east.

A similar velocity gradient is observed in the C¹⁸O (1-0) line, with velocities increasing from lower to higher values from the west to the east side. The bottom panel of Fig. 8 presents six spectral profiles at different random positions along the Snake filament for both ¹³CO (1-0) and C¹⁸O (1-0) lines. These profiles highlight the saturation of the ¹³CO (1-0) line, which prevents an accurate Gaussian fit. For positions 1 to 3 (first row of the plot), the C¹⁸O (1-0) line shows no detectable signal, whereas the ¹³CO (1-0) line exhibits a strong signal with two distinct velocity components. In positions 4 to 6 (second row of the plot), both lines show clear signals, each with a single velocity component. The peak velocities of the ¹³CO (1-0) and C¹⁸O (1-0) lines are in good agreement.

Figure 9 presents a 3D position-position-velocity diagram of the C¹⁸O (1-0) line and the two velocity components of the ¹³CO (1-0) line. The higher velocity component of the ¹³CO (1-0) line overlaps significantly with the velocity of the C¹⁸O (1-0) line. From Figs. 8 and 9, we conclude that the two gas tracers follow the same cloud and with similar velocities. The gas shows lower velocities at both ends of the Snake filament, with velocities increasing toward the central region of the filamentary structure, which corresponds to the denser part of this region. This pattern suggests that gas is accreting toward the central region, a behavior also observed in other filaments (e.g., Feng et al. 2024). As discussed in Sect. 3.3, the magnetic field lines are aligned with the filament’s spine at both ends. We infer from this that the magnetic field plays a role in guiding material toward the central region of the filament, where the gas density increases. In this denser region, the magnetic field morphology becomes more complex, leading to the observed orientation of magnetic field lines that appear perpendicular to the filament’s spine.

Fig. 7 (a) Channel maps between 4.01 and 5.44 km s−1 for the 13CO (1-0) spectral cube, in 0.13 km s−1 steps. White and cyan segments illustrate the B-field orientation in the optical (S/N≥5) and infrared (S/N≥3) bands, respectively. The velocity of each channel is shown on the top left of the panels. (b) Same as (a) but for the C18O (1-0) line. Channels are between 4.60 and 5.18 km s−1, in 0.07 km s−1 steps.

Fig. 8 Top panel: centroid velocity map of the C18O (1-0) line. Middle panel: Centroid velocity map of the 13CO (1-0) line. In both maps, we highlight six selected profiles at specific positions along the Snake, marked by white stars. The cyan contours represent the C18O (1-0) integrated intensity at levels of 0.8 K and 1.6 K. The beam size is shown in the bottom-left corner, and the scale bar in the bottom-right. Bottom panel: comparison of the six spectral profiles corresponding to the positions marked by white stars in both maps. The dashed orange profiles represent the C18O (1-0) line multiplied by 5, and the solid purple profiles correspond to the 13CO (1-0) line. The dashed black lines indicate the peak velocity of the 13CO (1-0) line, with the velocity values labeled next to them.

Fig. 9 Position-position-velocity plot illustrating the kinematics of the 13CO (1-0) and C18O (1-0) lines. Each data point represents the spatial location and centroid velocity of a Gaussian component, with colors distinguishing different Gaussian fits. The centroid velocity of the C18O (1-0) line is shown in orange, while the purple and red points correspond to the brighter and fainter intensity components of the 13CO (1-0) line, respectively.

3.5 Line opacity

We estimated the opacity of ¹³CO (1-0) and C¹⁸O (1-0) using the following expression: $$\tau =-\ln [1-\frac{T_{\mathrm{m}\mathrm{b}}}{J_{\nu }(T_{\mathrm{e}\mathrm{x}})-J_{\nu }(T_{\mathrm{b}\mathrm{g}})}],$$(5)

where T_ex = 16 K is the excitation temperature derived from the Herschel dust temperature map, T_mb is the peak main beam temperature of the line, J_ν represents the equivalent Rayleigh-Jeans temperature, and T_bg = 2.73 K is the temperature of the cosmic microwave background. To estimate uncertainties, we varied T_ex within a physically meaningful range, ensuring that the lower limit does not exceed the minimum excitation temperature, T_ex = 15 K. The optical depth of C¹⁸O (1-0) is determined to be $\tau ={0.6}_{-0.1}^{+0.1}$, calculated using T_mb = 5.5 K. The uncertainties were estimated by considering a ±1 K variation in T_ex. Similarly, for the ¹³CO (1-0) line, we obtain $\tau ={2.1}_{-0.9}^{+0.4}$ using T_mb = 11 K. These results confirm that the C¹⁸O (1-0) gas is optically thin, whereas the ¹³CO (1-0) gas is optically thick. The C¹⁸O (1-0) line has a relatively low critical density (n_c ~ 10³ cm⁻³), making it reasonable to consider that the excitation temperature (T_ex) is approximately equal to the kinetic temperature (T_k), considering dust-gas coupling. Although the Herschel dust temperature is typically higher than the excitation temperature, as mentioned earlier, we also calculated the line opacity for a variation in T_ex as an uncertainty.

4 Discussion

4.1 Filament stability and absence of star formation

Filament stability can be assessed by comparing their mass-to-length ratio to a critical value (Ostriker 1964; Fiege & Pudritz 2000). To incorporate the effect of nonthermal motions, we followed the formulation introduced by Hacar et al. (2023): $${\left(\frac{M}{L}\right)}_{\mathrm{c}\mathrm{r}}=m_{\mathrm{v}\mathrm{i}\mathrm{r}}({\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t}})=\frac{2{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t}}^2}{G}\sim 465{\left(\frac{{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t}}}{1\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)}^2{\mathrm{M}}_{\odot }{\mathrm{p}\mathrm{c}}^{-1},$$(6)

where G is the gravitational constant and M_⊙ is the solar mass. σ_tot is the total gas velocity dispersion, ${\sigma }_{\text{tot\hspace{0.17em}}}^2={\sigma }_{\text{nt\hspace{0.17em}}}^2+c_{\mathrm{s}}^2$. For the detail calculation of the σ_tot, see Appendix A.

If the actual mass-to-length ratio of the filament, ${\left(\frac{M}{L}\right)}_{\text{fil\hspace{0.17em}}}$, exceeds the critical value, the filament will collapse radially. To estimate the mass of the filament, we used the following equation: $$M_{\mathrm{f}\mathrm{i}\mathrm{l}}=1.13\times {10}^{-7}\times \sum _{i=1}^n{\left(\frac{N\left({\mathrm{H}}_2\right)}{{\mathrm{c}\mathrm{m}}^{-2}}\right)}_i{\left(\frac{\theta }{\mathrm{\prime }\mathrm{\prime }}\right)}^2{\left(\frac{d_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d}}}{\mathrm{p}\mathrm{c}}\right)}^2{\mu }_{{\mathrm{H}}_2}{m}_{\mathrm{H}}{M}_{\odot },$$(7)

where the summation ${\sum }_{i=1}^n\left(N{\left({\mathrm{H}}_2\right)}_i\right)$ includes contributions from n pixels within the filament area as introduced in Sect. 3.3, m_H = 1.67 × 10⁻²⁴ g is the mass of a hydrogen atom, μ_H2 = 2.8 (from Kauffmann et al. 2008) is the mean molecular weight per hydrogen molecule, and d_cloud = 154 pc is the distance to the cloud, θ = 10″ is the pixel size. We calculated the gas column density of the filament, N(H₂), from the visual extinction, A_V (Fig. 3). Using the relation provided by Bohlin et al. (1978), N(H₂) = 9.4 × 10²⁰ × A_V cm⁻² mag⁻¹, we converted A_V into the molecular hydrogen column density. The filament length is estimated to be L = 2.4 pc, with a detailed description of the measurement method provided in Sect. 3.3. The mass per unit length of the filament is determined to be 14.4 M_⊙ pc⁻¹, which is below the critical mass per unit length, ${\left(\frac{M}{L}\right)}_{\text{cr\hspace{0.17em}}}=31.3M_{\odot }{\mathrm{p}\mathrm{c}}^{-1}$. This indicates that the Snake filament is stable against gravitational collapse. This finding is valid assuming negligible external pressure onto the filament. The magnetic field would only increase the support against collapse. Consequently, this is likely the primary reason for the absence of star formation in this region.

We checked for the presence of recent star formation activity using version 7.1 of the SPHEREx Target List of ICE Sources (SPLICES), which is part of the SPHEREx all-sky survey mission (Ashby et al. 2023). The catalog contains 8.6 × 10⁶ objects that are brighter than approximately 12 Vega magnitudes in the W2 band. To identify young stellar objects (YSOs), we used the W1[3.6 mag]–W2[4.5 mag] versus W2[4.5 mag]–W3[5.8 mag] color-color diagram of asymptotic giant branch candidates from SPLICES (see Ashby et al. 2023). The regions defined by Koenig & Leisawitz (2014) for the selection of Class I and Class II YSOs. No YSOs have been detected in the area using the Spitzer GLIMPSE color-color diagram. We utilized the mid-infrared colors [3.6–4.5] versus [4.5–5.8] to identify YSOs in the region. We considered a region of size 2 pc that encompasses the filamentary features associated with the Snake region, shown in Fig. 3, and find 822 objects that are detected in the region (black dots in Fig. 10). We applied a distance constraint, selecting objects that are 1 kpc away, and we only find 19 objects within 1 kpc. Among these, we identify only one object as Class II, which is situated far from the Snake region.

Fig. 10 W1–W2 versus W2–W3 color–color diagram of asymptotic giant branch candidates from the SPLICES survey. The areas within the dashed lines are those used by Koenig & Leisawitz (2014) to identify Class I and Class II YSOs. Black dots represent all the data from SPLICES, and purple diamonds indicate objects within a distance of 1 kpc.

4.2 Polarization efficiency

Figure 11 presents the scatter plots of polarization efficiency (p_pol/A_V) as a function of visual extinction (A_V) in log-log space. The optical and NIR datasets are organized into 30 bins. The error bars on each data point are calculated as the standard deviation of the polarization measurements within each bin, reflecting the measurement uncertainties in the analysis. The optical and NIR data show a peak density within the visual extinction range of 1–3 magnitudes.

The plot includes two fitted lines corresponding to two datasets. The optical fit is represented by the equation: p_pol/A_V = −(0.73 ± 0.09) log A_V + 0.04. Conversely, the NIR fit is described by p_pol/A_V = −(0.75 ± 0.27) log A_V + 0.02, which shows a slightly steeper decline in polarization efficiency with increasing visual extinction; however, the two follow the same overall trend, as they span very similar ranges of extinction values. The optical polarization is primarily associated with the diffuse cloud surrounding the Snake filament, while the infrared data predominantly trace slightly denser regions. However, both datasets in our analysis correspond to relatively low-density regions within the Snake filament, spanning a visual extinction range of 0.5 to 4.5 magnitudes (see Fig. 3). The polarization efficiency exhibits a decrease at both optical and infrared wavelengths. The slopes observed for both datasets indicate depolarization in this region, consistent with the radiative alignment theory (Lazarian & Hoang 2007). As visual extinction increases, reduced radiation penetration within the cloud diminishes the effectiveness of radiative torques in aligning dust grains. As reported by Alves et al. (2014), NIR polarization measurements of the starless core Pipe-109 in the Pipe Nebula yield α (which is the slope of the fit) values of 1.0 for A_V < 9.5 mag and 0.34 for A_V > 9.5 mag. Redaelli et al. (2019) determined a steep slope of α = 1.21 for the FIR polarization data associated with the protostellar core IRAS 15398-3359. Similarly, Tabatabaei et al. (2024) derived a slope of 0.75 from NIR data of a filamentary structure in Barnard 59. However, other mechanisms may also contribute to this depolarization trend. Seifried et al. (2019) conducted radiative transfer modeling using the POLARIS code and show that dust grains remain well aligned even at high densities (n > 10³ cm⁻³) and extinctions (A_V > 1), provided the radiative alignment theory alignment is active. They found that the frequently observed decrease in polarization degree is primarily caused by large variations in the magnetic field orientation along the line of sight, rather than a breakdown of grain alignment. Wang et al. (2019) reported a power-law slope of α = 0.56 from 850 μm polarization observations of the IC 5146 filament, which implies that dust grains in this A_V ~ 20–300 mag range can still be aligned with magnetic fields.

Fig. 11 Scatter plots of the polarization efficiency (ppol/AV) as a function of the visual extinction (AV) in magnitudes in logarithmic scale. Measurement uncertainties are displayed as an error bar for each data point. Solid lines are the best fit for the datasets, as explained in the main text. Contours represent levels of constant probability density estimated using a 2D kernel density estimation of the data (using the seaborn Python library; Waskom 2021). These contours indicate where the data points are most densely concentrated. Both sets of contours span density levels from 0 to 0.8, with steps of 0.1, effectively highlighting regions of high and low data concentration. The best-fit parameters are shown in the top-right corner.

5 Summary

In this work, we investigated the magnetic field geometry and kinematic properties of the Snake filament by combining archival submillimeter polarimetry from Planck and molecular line data from IRAM with new optical and NIR polarization measurements. Our analysis confirms that the Snake is the dominant dust structure along the line of sight within 2 kpc of the Sun, with no significant molecular structures located beyond this distance. As a result, the observed polarization angles of background stars trace the magnetic field associated with the Snake filament. Additionally, we determined the distance of the source to be 154 pc using Gaia data (see Sect. 3.1). This distance puts the Snake filament at the same distance as the Pipe Nebula.

The Snake filament is gravitationally stable. The estimated mass-to-length ratio is below the critical threshold, explaining the lack of star formation in the region. The polarization efficiency (p_pol/A_V) decreases with increasing visual extinction, with fitted slopes of 0.73 for optical data and 0.75 for NIR data. This behavior suggests depolarization in denser regions, where reduced radiation penetration weakens the dust–grain alignment.

Kinematic analysis reveals a velocity gradient along the filament in both ¹³CO (1-0) and C¹⁸O (1-0) data. The C¹⁸O (1-0) line, tracing denser regions, exhibits well-constrained single-Gaussian fits, while the ¹³CO (1-0) line shows more complex kinematics, including multiple velocity components in the denser central region. Furthermore, the polarization vectors along the filament’s spine exhibit a strong alignment with the filament’s spine tangent at both ends. Consequently, it appears that the field may help channel material toward the central region, where the gas density increases. In this denser central area, the magnetic field morphology appears more complex, resulting in field lines that are oriented perpendicular to the filament’s spine.

These results demonstrate that the Snake filament is immersed in a well-ordered magnetic field, which plays a critical role in maintaining its structural coherence and stability. The observed polarization and kinematics suggest that magnetic fields significantly influence the dynamics of the filament, shaping its evolution. This study contributes to our understanding of the role of magnetic fields in filamentary structures and their connection to star formation processes.

Acknowledgements

Elena Redaelli acknowledges the support from the Minerva Fast Track Program of the Max Planck Society. G. Franco acknowledges the partial support from the Brazilian agencies CNPq and FAPEMIG. We thank the staff of OPD/LNA (Brazil) for their hospitality and invaluable help during our observing runs. This research has used data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome, and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The authors express their gratitude to Jaime Pineda for his valuable support with the data reduction process of the IRAM data.

Appendix A Velocity dispersions

The linewidth of the observed C¹⁸O (1-0) molecular line was determined by fitting Gaussian profiles to each spectrum in the data cube. Figure A.1 illustrates the derived velocity dispersion map, σ_v, highlighting the narrow linewidths of the C¹⁸O (1-0) emission, with values typically around 0.1 km s⁻¹ to 0.2 km s⁻¹. The isothermal sound speed was obtained from $$c_{\mathrm{s}}=\sqrt{\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}}}=0.2{\mathrm{k}\mathrm{m}\mathrm{s}}^{-1},$$(A.1)

where μ = 2.37 (Kauffmann et al. 2008) is the mean molecular weight per free particle and m_H is the mass of hydrogen atom. We used the dust temperature, T_d = 15 K from Herschel map as a proxy for the gas kinetic temperature.

The nonthermal velocity dispersion of the gas can be derived from $${\sigma }_{\mathrm{n}\mathrm{t}}=\sqrt{{\sigma }_{\mathrm{v}}^2-{\sigma }_{\mathrm{t}\mathrm{h}}^2}=0.14{\mathrm{k}\mathrm{m}\mathrm{s}}^{-1},$$(A.2)

with ${\sigma }_{\text{th\hspace{0.17em}}}=\sqrt{k_{\mathrm{B}}T/m_{\text{obs\hspace{0.17em}}}}=0.06$ km s⁻¹ and σ_v = 0.15 km s⁻¹, where m_obs is the mass of observed molecule line, C¹⁸O (1-0). As mentioned in the main text, the total gas velocity dispersion is calculated using $${\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t}}=\sqrt{{\sigma }_{\mathrm{n}\mathrm{t}}^2+c_{\mathrm{s}}^2}=0.24{\mathrm{k}\mathrm{m}\mathrm{s}}^{-1}.$$(A.3)

Fig. A.1 Velocity dispersion (σV) map of the C18O (1-0) line. The cyan contours represent the C18O (1-0) integrated intensity at levels 0.8 and 1.6 K. The beam size is indicated in the bottom-left corner, and the scale bar in the bottom-right.

References

All Figures

	Fig. 1 2D dust extinction map illustrating the dust distribution as a function of distance up to 1.25 kpc from the Sun (Edenhofer et al. 2024). The peak at 154 pc indicates a distance corresponding to the Snake region.
In the text

Fig. 2 Top panel: visual extinction (AV) as a function of the distance of stars within 2 kpc of the Sun. The extinction and distance values are derived from the StarHorse catalog. Lower panel: polarization percentage as a function of the distance. A vertical dashed line at 154 ± 15 pc marks the distance of the Snake filament. The error bars represent the uncertainties in the extinction and polarization measurements, scaled by 50%. The purple line represents the mean value in bins of 100 pc, with the light blue shaded area indicating the ±1σ uncertainty.

In the text

Fig. 3 Dust extinction map of the Snake region in color scale at a spatial resolution of 20″ from Román-Zúñiga et al. (2010). All segments represent the B-field orientation. Blue and red segments are inferred from the optical and NIR data, respectively, with S/N≥3. The green segments show the magnetic field orientation from Planck data at the resolution of 10′. The vectors are scaled to the same size to enhance visualization clarity. The cyan contours represent the C18O (1-0) integrated intensity at levels of 0.8 K km/s and 1.6 K km/s.

In the text

	Fig. 4 Histogram of the distribution of magnetic polarization angles for the three datasets (optical, NIR, and submillimeter). Optical data peaks are highlighted by dashed blue lines and the submillimeter dataset distribution peaks by dashed green lines.
In the text

Fig. 5 Dust extinction map illustrating the masked filamentary structure. The spine of the filament is outlined by a prominent red curve, and the thin red lines indicate the locations of perpendicular cuts. Blue circles represent the peak pixel intensity along each cut. The orange vectors represent optical (S/N≥5) and infrared (S/N≥3) polarization segments. The cyan vectors are projections of all the orange polarization vectors onto the spine. Each polarization vector was projected by finding its closest point on the spine and shifting it accordingly.

In the text

	Fig. 6 PRS values as a function of RA. The PRS is calculated at each position along the filament spine using the group of polarization angles projected onto that location. Each point represents the PRS for a given spine position. The dotted line indicates PRS = 0, while the error bars represent the uncertainties.
In the text

Fig. 7 (a) Channel maps between 4.01 and 5.44 km s−1 for the 13CO (1-0) spectral cube, in 0.13 km s−1 steps. White and cyan segments illustrate the B-field orientation in the optical (S/N≥5) and infrared (S/N≥3) bands, respectively. The velocity of each channel is shown on the top left of the panels. (b) Same as (a) but for the C18O (1-0) line. Channels are between 4.60 and 5.18 km s−1, in 0.07 km s−1 steps.

In the text

Fig. 8 Top panel: centroid velocity map of the C18O (1-0) line. Middle panel: Centroid velocity map of the 13CO (1-0) line. In both maps, we highlight six selected profiles at specific positions along the Snake, marked by white stars. The cyan contours represent the C18O (1-0) integrated intensity at levels of 0.8 K and 1.6 K. The beam size is shown in the bottom-left corner, and the scale bar in the bottom-right. Bottom panel: comparison of the six spectral profiles corresponding to the positions marked by white stars in both maps. The dashed orange profiles represent the C18O (1-0) line multiplied by 5, and the solid purple profiles correspond to the 13CO (1-0) line. The dashed black lines indicate the peak velocity of the 13CO (1-0) line, with the velocity values labeled next to them.

In the text

Fig. 9 Position-position-velocity plot illustrating the kinematics of the 13CO (1-0) and C18O (1-0) lines. Each data point represents the spatial location and centroid velocity of a Gaussian component, with colors distinguishing different Gaussian fits. The centroid velocity of the C18O (1-0) line is shown in orange, while the purple and red points correspond to the brighter and fainter intensity components of the 13CO (1-0) line, respectively.

In the text

	Fig. 10 W1–W2 versus W2–W3 color–color diagram of asymptotic giant branch candidates from the SPLICES survey. The areas within the dashed lines are those used by Koenig & Leisawitz (2014) to identify Class I and Class II YSOs. Black dots represent all the data from SPLICES, and purple diamonds indicate objects within a distance of 1 kpc.
In the text

Fig. 11 Scatter plots of the polarization efficiency (ppol/AV) as a function of the visual extinction (AV) in magnitudes in logarithmic scale. Measurement uncertainties are displayed as an error bar for each data point. Solid lines are the best fit for the datasets, as explained in the main text. Contours represent levels of constant probability density estimated using a 2D kernel density estimation of the data (using the seaborn Python library; Waskom 2021). These contours indicate where the data points are most densely concentrated. Both sets of contours span density levels from 0 to 0.8, with steps of 0.1, effectively highlighting regions of high and low data concentration. The best-fit parameters are shown in the top-right corner.

In the text

	Fig. A.1 Velocity dispersion (σV) map of the C18O (1-0) line. The cyan contours represent the C18O (1-0) integrated intensity at levels 0.8 and 1.6 K. The beam size is indicated in the bottom-left corner, and the scale bar in the bottom-right.
In the text