© The Authors 2024

1 Introduction

Stars form in cold dense cores embedded in molecular clouds. These cores are characterised by higher densities and lower temperatures, compared to the parental cloud (Myers 1983; Myers & Benson 1983; Caselli et al. 2002). Studies with high-density (n(H₂) > 10⁴ cm⁻³) tracers have revealed subsonic levels of turbulence inside cores (Barranco & Goodman 1998; Rosolowsky et al. 2008), in contrast to the supersonic linewidths observed in the ambient cloud (traced by lower density tracers, like CO; Larson 1981).

A useful tracer for observing dense cores in molecular clouds is NH₃. Even though it has a relatively low critical density (a few times 10³ cm⁻³), the hyperflne structure of the NH₃ inversion transitions allows for the individual hyperflnes to remain optically thin even at high column densities (Caselli et al. 2017). Fitting the multiple hyperflnes of NH₃ simultaneously also gives a very precise constraint on the kinematical information. Moreover, if the NH₃ (2,2) line is reliably detected along with the (1,1), it allows for a direct measurement of the gas temperature (e.g. Friesen et al. 2017).

Previous observations have shown filaments to be widely present in molecular clouds in star-forming regions and connected to subsonic dense cores (Goldsmith et al. 2008; André et al. 2010; Polychroni et al. 2013). These filamentary structures have been extensively observed using Herschel (Ward-Thompson et al. 2010; André et al. 2014; Arzoumanian et al. 2019) and studied with respect to a range of physical properties, including density and temperature. Molecular line observations of filaments have revealed complex kinematical structures (Kirk et al. 2007; Henshaw et al. 2013; Hacar et al. 2013). Using molecular lines, filaments have been observed to have much smaller lengths (Hacar et al. 2017; Monsch et al. 2018; Suri et al. 2019) than the filaments identified with continuum observations. These filaments have been observed with lengths ranging from 0.01 pc to 0.3 pc; therefore, they appear to be scale-free.

Chen et al. (2020a) suggested the use of a dimensionless parameter, C_v, to distinguish between two scenarios resulting in the formation of filaments, namely, turbulent-driven and gravity-dominated compression. This parameter, which depends on the velocity gradient orthogonal to the major axis of the filament and the mass per unit length of the filament, has been used to differentiate between the formation scenarios of filaments in observations of star-forming regions (e.g. Zhang et al. 2020; Gong et al. 2021; Hsieh et al. 2021). However, in these studies, a single value is calculated representing the entire filament and no distribution along the filament is reported. Analyses of velocity gradients along the filament could highlight the robustness of this parameter, as well as to study the overall structure of the filaments.

We study the subsonic core Barnard 5 (hereafter B5) in the Perseus molecular cloud with combined Very Large Array (VLA) and Green Bank Telescope (GBT) observations using NH₃ (1,1) and (2,2) transitions. B5 is located at a distance of 302 ± 21 pc (Zucker et al. 2018), and hosts a young stellar object (YSO), the Class-I protostar B5-IRS1 (Fuller et al. 1991). Using NH₃ hyperfine transitions, Pineda et al. (2010) revealed the turbulence inside the core to be subsonic and nearly constant, with a sharp transition to coherence at the core boundary. Further interferometric observations with high spatial resolution (6", ≈1800 au) revealed two filaments embedded in the subsonic core (Pineda et al. 2011) and three gravitationally bound condensations (Pineda et al. 2015; Schmiedeke et al. 2021). These filaments are similar in length to those presented in Suri et al. (2019). However, it is not yet well understood how such subsonic filaments within dense cores relate to their large-scale counterparts. In this project, we use this high-resolution data, smoothed to a slightly larger beam (8"), to study the detailed velocity structure within B5 and the embedded filaments with multi-component analysis.

Fig. 1 Integrated intensity of NH3 (1,1) in Barnard 5. The black-solid and black-dashed contours show the two filaments and the three condensations in B5, respectively (contours and nomenclature adapted from Schmiedeke et al. 2021). The beam size and scale-bar are shown in the bottom-left and the top-right corners, respectively. The white-solid contour shows the core boundary calculated using single-dish observations by Pineda et al. (2010).

2 Ammonia maps

We used single-dish observations of B5 with the GBT combined with interferometric observations with the VLA of NH₃ inversion transitions (1,1) and (2,2). The GBT observations were performed under project 08C-088 with on-the-fly (OTF) mapping (Mangum et al. 2007) in frequency-switching mode. The data have a spectral resolution of about 3 kHz (or 0.04 km s⁻¹). More details of the data processing can be found in Pineda et al. (2010), where the data were first published. The interferometric data were obtained with the VLA in the D-array configuration on 16–17 October 2011 and in the DnC-array configuration on 13– 14 January 2012 (project 11B-101). The spectral resolution of the data was 0.049 km s⁻¹. A detailed description of the data can be found in Pineda et al. (2015). The resolution of the combined data is 6″.

We smoothed the data to a beam of 8″ to increase the sensitivity, which enables us to identify and fit multiple components towards the filaments. To avoid oversampling, we fixed the pixel-to-beam ratio to 3. In this work, we follow the nomenclature of Pineda et al. (2015) and Schmiedeke et al. (2021) regarding the different regions in B5 clump 1, filaments 1 and 2, and condensations 1, 2, and 3). Figure 1 shows the integrated intensity map of NH₃ (1,1) in B5, along with the locations of the two filaments and the three condensations inside the subsonic core.

3 Analysis

3.1 Line fitting

We used the cold_ammonia model in the python package pyspeckit (Ginsburg & Mirocha 2011; Ginsburg et al. 2022) to fit the NH₃ (1,1) and (2,2) data¹. This model is appropriate for regions with temperatures <40 K (see Friesen et al. 2017). The model fits the NH₃ (1,1) and (2,2) spectra simultaneously with kinetic and excitation temperatures (T_K and T_ex), NH₃ column density (N(NH₃)), velocity dispersion (σ_v), centroid velocity (v_LSR), and the ortho-NH₃ fraction of the total NH₃column density (f_ortho) as parameters.

We followed the Bayesian approach detailed in Sokolov et al. (2020) to fit multiple components and selected the optimum number of components in the model. In classical line fitting procedures, the best-fit model is obtained by minimising the χ² parameter: $${\chi }^2(\theta )={\displaystyle \sum _i\frac{{\left(I_i-{ℳ}_i(\theta )\right)}^2}{{\sigma }^2}},$$(1)

where I_i and 𝓜_i(θ) are the observed and model intensity in channel i, respectively, and σ is the noise in the data, while θ represents the parameters (see above) that define the model. χ² is related to the likelihood function, 𝓛(θ), as follows: $$\ln ℒ(\theta )=\frac{{\chi }^2}{2}+\text{ const}.$$(2)

In the model selection employed in our work, the Bayes factor $K_b^a$ gives the estimate of whether model 𝓜_a is favoured over model 𝓜_b: $$K_b^a=\frac{P\left({ℳ}_b\right)Z_a}{P\left({ℳ}_a\right)Z_b}$$(3)

where P(𝓜_l)s are the posterior probability densities and Z_ls are the likelihood integrals over the parameter space. The likelihood integral, Z, is given by: $$Z={\displaystyle {\int }_{\theta }p}(\theta )ℒ(\theta )\text{d}\theta ,$$(4)

where p(θ) are the prior probability densities of the parameters represented by θ. The posterior probability densities can be calculated using the likelihood function as: $$P(\theta )\propto p(\theta )\times ℒ(\theta )$$(5)

We adopted a threshold of $K_b^a=5$ to indicate whether model 𝓜_a is preferred over 𝓜_b (following Sokolov et al. 2020). Figure A.1 shows the detailed maps of the factors $K_0^1,K_1^2$, and $K_2^3$

The initial priors were estimated to accommodate the ranges of these parameters in B5 known from previous works with single-component fits (Pineda et al. 2011, 2021; Schmiedeke et al. 2021). The ranges of the priors were suitably increased if they were found to be inadequate for the multi-component fits. As the priors in the model, we use uniformly distributed ranges of 6–15 K for T_K, 4–12 K for T_ex, 10¹³·²–10¹⁵ cm⁻² for $N_{{\text{NH}}_3},$, 0.04–0.3 km s⁻¹ for σ_v and 9–11 km s⁻¹ for v_LSR. Since we only have detections of the ortho-NH₃(1,1) and (2,2) lines, we cannot fit for ƒ_ortho. Therefore, for each fit, the ortho-ammonia fraction was fixed at 0.5 (for an assumed ortho-para ratio of unity). Additionally, in order to reduce the likelihood of spurious components, the minimum separation between two components was required to be 0.1 km s⁻¹ (approx. twice the spectral resolution); while the maximum separation between two components was set at 2 km s⁻¹ (more than double the maximum range in velocity across B5 from single component fit results, see Pineda et al. 2010; Schmiedeke et al. 2021) to reduce the parameter space to be covered in the likelihood calculations (this was achieved by suitable priors for the additional components).

After the initial runs, we identified that all pixels with a good 1-, 2- and 3-component fits are within a signal-to-noise ratio (S/N) of 4, 7, and 10, respectively², in the peak intensity of the NH₃ (1,1) line. We therefore apply these thresholds in the parameter maps presented here. We attempted to fit up to four components to the data, and found that three was the maximum number of components in the best-fit model for regions larger than a beam. Figure 2 shows the number of components in the best-fit model in different parts of B5. Throughout the core, we are able to fit at least one component (see Fig. 2). Additionally, towards the two filaments and clump 1, we could fit two components and in some regions inside the filaments, we find that a three-component fit is the optimal model.

Fig. 2 Number of components detected towards different parts of B5, as decided by the Bayesian approach described in Sect. 3.1 (left). These components are then subdivided into 5 separate components according to their velocities (see Sect. 3.2). The black-solid and black-dashed contours show the two filaments and the three condensations in B5, respectively (contours adapted from Schmiedeke et al. 2021). The black-dotted contour shows the boundary of the coherent core in Pineda et al. (2010). The beam size and scale-bar are shown in the bottom-left and the top-right corners, respectively. Examples of NH3 (1,1) spectra from regions with one-, two-, and three-component fits (right, with the positions of the spectra indicated in the left panel). Only a small range in velocity is shown here to clearly highlight the individual components, which are illustrated in blue, red, and purple. In each panel, the fitted model is shown in green, while the residuals are shown in grey. Note: there are multiple, closely-spaced hyperfines in the velocity range shown here. The strongest two of these hyperfines can be seen in the narrow component (blue and purple in panels 2 and 3). However, these hyperfines are blended together in the broad component (green in panel 1, red in panels 2 and 3), and therefore the hyperfine structure is not visible.

Fig. 3 Velocity map of ‘global component’ (left). Velocity map of the four additional components (right). Each of these four components show a very narrow range of velocities. The median velocities of these components are 9.61 km s−1 (far-blue), 9.99 km s−1 (mid-blue), 10.32 km s−1 (mid-red), and 10.58 km s−1 (far-red). The black-solid and black-dashed contours show the two filaments and the three condensations in B5, respectively. The common colour scale is shown in the left panel. The black-dotted contours show the boundary of the coherent core in Pineda et al. (2010). The scale-bar and beam size are shown in the top-right and bottom-left corners, respectively.

3.2 Component assignment

Figure 2 shows the number of components fit towards different parts of B5. Since each pixel is fitted independently, and the number of components varies across the map, the identification of spatially coherent features is non-trivial. With the exception of the two filaments and clump 1 (see Fig. 1), we detected a single continuous component in the majority of the coherent core. We first isolated this component with an approach similar to that of Chen et al. (2022). In regions with two or three components, we searched for the fitted component which provides the smoothest map of velocity and the brightness temperature, when combined with the single component in the outer region (see Appendix B for details). This gives us a map of the component that is kine-matically coherent over the largest spatial extent³. Hereafter, this component is referred to as the ‘global component’, as it is present throughout all the regions discussed in this work. The velocity and velocity dispersion of this component are shown in the left panels in Figs. 3 and 4. This component shows very similar properties as those presented in the single dish analysis of the region, which is more sensitive to the extended emission, (see Fig. 2 in Pineda et al. 2010).

We are then left with one or two more components in positions close to or towards the filaments. Upon visual inspection of these two components, we find that they are either red- or blue-shifted with respect to the median velocity of the global component (v_med = 10.15 km s^–1). Therefore, we grouped the remaining two components into red-shifted and blue-shifted components. Each of these two components shows indications of further subdivision in velocity, and therefore, they are split into two further components each to give four kinematically coherent maps. The velocities and velocity dispersions of these four components are shown in the right panels of Figs. 3 and 4. Due to insufficient sensitivities, we were not able to detect all of these components continuously, hence, these four maps appear patchy.

We note that these different velocity components are not shown separately in Fig. 2, which only shows the number of components in the best-fit model at each pixel.

Chen et al. (2022) presented a kinematical analysis of B5 using the same data set at the original resolution of 6". However, they take a different statistical approach for the multi-component fit to the data using the MUFASA software (Chen et al. 2020b). Moreover, Chen et al. (2022) attempted to fit up to only two velocity components in the region and focus their analysis on the one component that is kinematically coherent over the core region (very similar to ‘global component’ in this work), while the other component is excluded. In contrast, we focus our discussion on the four additional velocity components that we extract with the Bayesian approach to multi-component fit. Therefore, a direct comparison between the two works is not straightforward, even though both present an analysis of the signatures that hint at infalling material in the B5 filaments⁴.

Fig. 4 Velocity dispersion of the ‘global component’ (left). Velocity dispersions of the four additional components (right). The white-solid and white-dashed contours show the two filaments and the three condensations in B5, respectively. The white/black-dotted contours show the boundary of the coherent core in Pineda et al. (2010). The scale-bar is shown in the top-right corner and the beam size is shown in the bottom-left corner.

4 Results and discussions

4.1 Identified velocity components

As can be seen from the left panel of Fig. 3, the ‘global component’ shows a smooth velocity field, fluctuating in the north-south direction. As this is also the brightest component, it dominates the velocity obtained from a single component fit to the data. Therefore, this velocity map is very similar to those obtained with a single component fit (Pineda et al. 2011; Schmiedeke et al. 2021). Inside the core boundary (from Pineda et al. 2010, dotted line in Fig. 3), we attribute this component to the dense core material, and outside the core boundary, to the ambient cloud. This can be clearly seen from the velocity dispersion in the component (Fig. 4). Inside the core boundary, the dispersion is subsonic (≲0.2 km s^–1) and outside the core, it is supersonic.

The four panels on the right of Fig. 3 show the velocities of the additional components we identify in B5. We refer to these components as ‘far-reď, ‘mid-reď, ‘mid-blue’, and ‘far-blue’ according to their shift in velocities (as shown in the figure). We are not able to detect all four components towards all positions, with a maximum of two (out of these four) components being detected at a given location. This is most likely due to insufficient sensitivity in the data.

The mid-blue and the mid-red components seem to trace roughly the two filaments. Each component shows almost constant velocities (median velocities of 9.99 km s^–1 and 10.32 km s^–1, respectively) with a standard deviation less than 0.1 km s^–1. This suggests that each component is tracing a single, velocity coherent structure. The velocity dispersion of these two components is subsonic inside the two filaments (Fig. 4), which further shows that they are tracing material inside a subsonic core. Given the spatial extent and constant separation in velocity of these two components, their origin is most likely the two opposite sides of material undergoing infall onto the filaments (see Fig. 5).

The far-blue component shows an almost constant velocity around 9.6 km s^–1 (standard deviation of 0.16 km s^–1), which indicates that the component is tracing a velocity coherent structure. The velocity dispersion in this component is largely supersonic (Fig. 4), which suggests that this component traces gas in the ambient cloud, rather than the core. The velocity and the velocity dispersion of this component are similar to those of the global component in the north. Therefore, this component likely shows material flowing from the ambient molecular cloud towards B5 from the north and north-west. We explore this possible mass infall in Sect. 4.4. The location and velocity of this component is consistent with those of the methanol peak (Wirström et al. 2014; Taquet et al. 2017, observations with Herschel, Onsala, and IRAM 30 m observatories), which is indicative of gas collision. We note that there is a significant apparent gap in the map of this component between clump 1 and the central region. This could be due to some of the material in this component merging with the core material, making it difficult to distinguish from the core (global component). The lower velocity dispersions in the far-blue component towards the middle of the core, as compared to the northern part, could point to the turbulence in the component being dissipated in the process of merging with the core material.

The far-red component, which is red-shifted and at a velocity of approximately 10.6 km s⁻¹, could be largely divided into two parts: the arc-like structure just south of the boundary of filament-2, and region near condensation-1. The arc-like structure might show material falling into B5 from that direction similar to the aforementioned scenario towards north, as the velocities in this component are systematically different from that in the core and filament components (Fig. 3). The region just south of condensation-1 might trace material falling towards the condensation due to gravitational collapse (see also Chen et al. 2022). A corresponding region is observed in the far-blue component north of condensation 1, which could be tracing the opposite side of this infall. Figure 5 shows the arrangement of the far-blue, mid-blue, mid-red and far-red components in the line of sight, as discussed in this section.

Fig. 5 Schematic showing the spatial distribution of the different components in the line of sight. The coloured arrows show the direction of the velocity of the corresponding component. The ellipses represent the core and a filament inside the core, as indicated. The grey arrows indicate the north and south directions in the maps presented in this work.

4.2 Filament formation scenarios

Chen et al. (2020a) proposed an analysis to determine whether filament formation is dominated by self-gravitational collapse or via direct shock compression. They define the dimensionless parameter C_v to differentiate between the two cases as $$C_{\text{v}}\equiv \frac{\text{Δ}{\text{V}}_{\text{h}}^2}{GM(r)/L}$$(6)

where Δv_h is half of the velocity gradient across the filament out to a distance, r, from the centre, G is the gravitational constant, and M(r)/L is the mass per unit length of the filament at a location, r. Here, $\text{Δ}{\text{v}}_{\text{h}}^2$ is representative of the kinetic energy of the flow transverse to the filament, while the denominator in Eq. (6) represents the gravitational potential energy in the filament. Therefore, the case where C_v ≤ 1 and C_v ≫ 1 represents scenarios where gravity or turbulence and flow dynamics, respectively, dominates the filament formation.

We determine the boundaries and the spine of the filaments following the analysis of Schmiedeke et al. (2021). We define orthogonal cuts to each of the filament spines with a separation of one beam between each cut. These cuts are shown in Fig. 6. For each of these cuts we calculated the difference in the velocity at the two ends (at the filament boundaries) and the mass per unit length. As described in Sect. 4.1, the mid-red and mid-blue components likely trace infall onto the filaments in the line of sight.

For this analysis, we needed to measure a velocity gradient in the plane of sky, both inside and in the neighbourhood of the filaments. Therefore, we used the map of the global component (left panel in Fig. 3) here. We calculated the mass using an NH₃ intensity to mass conversion factor of 1.5 ± 0.7 M_⊙ (Jy beam⁻¹)⁻¹, derived in Schmiedeke et al. (2021) by comparing the NH₃ (1,1) emission and the 450 µm continuum flux from SCUBA-2⁵.

Figure 7 shows six selected examples of the gradients in velocity across the orthogonal cuts to the filaments 1 and 2. The corresponding cuts are shown with thick lines in Fig. 6. The gradients in velocity across all the cuts are shown in Appendix C. Overall, we find C_v values less than 0.03, and so, across these two filaments, C_v ≪ 1 (Fig. C.3).

The average Δv_h for the orthogonal cuts along the two filaments is 0.04 km s⁻¹. For the C_v values to be near unity, Δv_h, needs to be ≈O.5 km s^{– 1}, on average, across the cuts. Therefore, the observed Δv_h values are more than ten times lower than what is needed for C_v = 1 at the mass-per-unit-length across the filaments. For comparison, the sound speed at the temperatures around the filaments is smaller than 0.2 km s^–1. The gas in that region is subsonic, therefore, the observed velocity dispersions are even smaller. By using the components that are likely tracing infall onto the filaments in the line of sight (mid-red and mid-blue), we get Δv_h ≈ 0.15 km s⁻¹, which is also much lower than what is required to have C_v = 1.

As mentioned by Chen et al. (2020a), filaments formed in a planar-geometry (preferred direction of accretion) could have C_v ≪ 1 due to projection effects or if the filaments are at an early formation stage. The latter scenario is consistent with the prestellar condensations, but is not consistent with the presence of the Class 1 protostar in B5, nor the highly supercritical M/L of the filaments. In the former scenario, if the slab containing the filament is inclined at an angle, θ, with the plane of the sky, the observed velocity difference across the filament, Δv_h,obs will be $$\text{Δ}{\text{v}}_{\text{h},\text{obs}}=\text{Δ}{\text{v}}_{\text{h}}\sin \theta ,$$(7)

where Δv_h is the true velocity difference. The true C_v value would then be C_v,_obs/sin²θ, where C_v,_obs is what we report here. Considering the C_v values we obtain across the filaments, the inclination angle needs to be |θ| ≲ 7°, for them to be comparable to the values obtained by Chen et al. (2020a), considering that |θ| = 45° in their results. Therefore, it is very unlikely for the estimated low values of C_v to be entirely due to inclination effects.

Since we observe C_v ≪ 1 across the two filaments, we can rule out the turbulence/flow formation scenario. Schmiedeke et al. (2021) find mass-per-unit-length (M/L) in these two filaments to be much higher (~80 M_⊙ pc^–1) than what can be thermally supported (~17M_⊙ pc^–1). Therefore, both the filaments are strongly dominated by gravity. However, we do not see a corresponding strong gradient in the velocity field, which suggests the presence of additional support against gravity. A strong poloidal magnetic field that is parallel to the spine of the filament could provide this necessary support. Schmiedeke et al. (2021) estimated the required magnetic field strength to be ~500 µG. The estimated Alfven velocity, $V_{\text{A}}(=B/\sqrt{4\pi \rho },$, where B is the magnetic field strength, and ρ is the gas mass density) in the filaments is ~0.5 km s⁻¹, which is almost a factor of three higher than the observed velocity dispersions in the filaments (see Fig. 4). This estimated value of V_A is similar, however, to the Δv_h required to have C_v = 1; that is, the scenario where gravitational potential energy is comparable to the kinetic energy in the filaments. Therefore, it supports the additional support against gravity being provided by the estimated magnetic field.

We note that the values for C_v in our work are much smaller than those for the models considered by Chen et al. (2020a), who reported C_v ~ 0.3. This suggests that their simulations are not good representations of the filaments in B5. We discuss three possible factors behind the significant difference in the C_v values between our work and the simulations of Chen et al. (2020a) in Appendix D.

Fig. 6 Orthogonal cuts to the major axes of the two filaments. The cuts corresponding to the velocity gradients shown in Fig. 7 are shown with thick lines. The background greyscale shows the integrated intensity of NH3(1,1). The cuts are spaced one beam apart. The beam is shown in the bottom-left corner.

4.3 Infall in the filaments

Assuming that both mid-blue and mid-red velocity components (see Sect. 4.1) trace the two opposite sides of infalling material, we can estimate the rate of infall onto the two filaments. The infall rate is then given by: $$\dot{m}=S\times \rho \times v_{\text{infall }},$$(8)

where S is the surface area of the filament, ρ is the volume density inside the core, and υ_infall is the velocity of the infall, taken to be the mean difference between the red- and blue-shifted filament components (mid-red and mid-blue). As both the red-and blue-shifted filament components (mid-red and mid-blue components, respectively) have roughly constant velocities of 10.3 km s^–1 and 9.9 km s⁻¹, respectively, we take $${\text{v}}_{\text{infall\hspace{0.17em}}}=\frac{\langle {\text{v}}_{\text{red\hspace{0.17em}}}-{\text{v}}_{\text{blue\hspace{0.17em}}}\rangle }{2}=0.2\hspace{0.17em}\text{km}\hspace{0.17em}{\text{s}}^{-1}$$(9)

Figure 8 shows a schematic of the geometry considered to calculate the approximate surface area of the two filaments. In the case of filament-1, we estimate the surface area assuming the filament to be roughly cylindrical with a diameter d = 0.05 pc and a length l = 0.2pc. The surface area then is $S=2\pi \left(\frac{d}{2}\right)l$. For filament-2, we estimate the surface area by dividing the filament into two roughly cylindrical regions with diameters d₁ = 0.07 pc and d₂ = 0.05 pc, and with lengths l₁ = 0.14 pc and l₂ = 0.13 pc, respectively. The total surface area is $S=S1+S2=2\pi \left(\frac{d_1}{2}\right)l_1+2\pi \left(\frac{d_2}{2}\right)l_2$.

Using the filament masses from Schmiedeke et al. (2021) and the geometry described above, we obtain a volume density of $n_{{\text{H}}_2}=3.5\times {10}^5 {\text{cm}}^{-3}$ for filament-1 and $n_{{\text{H}}_2}=2.5\times {10}^5 {\text{cm}}^{-3}$ for filament-2. As the different components are in the same line of sight, the exact densities associated with the components are difficult to distinguish. Therefore, we use the average density of the filament to calculate the infall rate. With these estimates, we get infall rates of 1.6 × 10^–4 M_⊙ yr^–1 for filament-1 and 1.8 × 10⁻⁴ M_⊙ yr⁻¹ for filament-2. The total masses of the two filaments are 9.4 M_⊙ and 13.4 M_⊙, respectively (Schmiedeke et al. 2021). At these rates, the formation timescale of the two filaments, as well as the mass doubling time, is 6–7 × 10⁴ yr, which is an order of magnitude less than typical core lifetimes of 0.5-1 Myr (Offner et al. 2022). This implies that at these infall rates, the filament masses can change significantly during the evolution of the core. As the material likely undergoing infall is at the surface, the density is expected to be lower than the average density, and hence, these estimates are upper-limits for the infall rates. Using observations with HC₃N (which traces chemically fresh gas, as compared to NH₃; Suzuki et al. 1992) and comparison with the NH₃data, Valdivia-Mena et al. (2023) suggest infall of fresh material from the B5 core onto the filaments. This could correlate to the infall traced by the mid-red and mid-blue components, as discussed above.

Fig. 7 Selected examples of the line of sight velocity profiles in the orthogonal cuts in filament-1 (top) and filament-2 (bottom). The red-dotted vertical lines show the position of the filament spine and the black-dotted vertical lines show the boundary of the filament for the respective orthogonal cuts shown. The Δv parameter (Δvh in Eq. (6)) for the cut is shown in the upper-left corner. The Δvh value needed to have Cv = 1 $\left(\text{Δ}{\text{v}}_{C_V=1}\right)$ is shown in the bottom-right corner for the upper panels, and below the Δvh in the lower panels.

Fig. 8 Geometry considered to calculate the mass inſall onto filament-1 (top) and filament-2 (bottom).

4.4 Mass infall onto B5

Assuming that the far-blue component (top left panel on the right of Fig. 3) traces material falling onto B5, we can similarly estimate its rate of infall. We assume that this material is flowing through a cylindrical region shown in Fig. 9, although only one-third of it (approximately) is detected⁶.

We estimated the rate of infall of material onto B5, using Eq. (8). Again, as this component is in the same line of sight as the other components, it is difficult to obtain the exact density associated to the material. So, we assumed a volume density of $n_{{\text{H}}_2}~{10}^4 {\text{cm}}^{-3},$, characteristic for core material, for the infall rate calculation. As the mean velocity in the core around this region is ≈ 10.2 km s^{– 1} (Fig. 3), we take the infall velocity to be 0.6 km s^–1 (as the far-blue material is at υ = 9.6 km s⁻¹). We then obtain a mass infall rate of 3.5 × 10⁻⁵ M_⊙ yr⁻¹ onto B5. We note that as this component is sparsely sampled, this is an upper limit of the mass infall rate.

Using observations with HC₇N, Friesen et al. (2013) calculated a lower limit of 2–5 × 10⁻⁶ M_⊙ yr^–1 for accretion onto dense filaments in the Serpens South, which is a cluster-forming region. As HC₇N traces a lower density compared to NH₃, the infall observed by Friesen et al. (2013) is likely in a different, more external layer in the cloud, than what would be traced by NH₃. The observation of different regions of mass infall with comparable rates might point to the existence of a mechanism through which new material is being delivered to coherent cores. In Choudhury et al. (2021), we also discuss possible material accretion onto coherent cores. These results highlight the possibility of observing infall onto cores in different star-forming regions, using suitable tracers of different densities. If accretion of new material onto the core is observed towards more regions, there could be a need for a significant update to the current core evolution models, where cores are assumed to evolve in isolation.

Fig. 9 Geometry considered for calculating the mass infall onto B5.

5 Conclusions

We used combined GBT and VLA observations of the star-forming core Barnard 5 in Perseus to study the multi-component structure of the core. Our results can be summarised as follows:

• 1.

  By smoothing the data to a beam of 8" (the original beam was 6"), we are able to identify different velocity components in the line of sight. Apart from the brightest component (which has properties similar to those from a single component fit), we detected four separate velocity components. We label these four components ‘far-red’, ‘mid-red’, ‘mid-blue’, and ‘far-blue’, alluding to their redshift or blueshift relative to the average velocity of the cloud;

• 2.

  The mid-red and mid-blue components are separated by 0.4 km s^–1 and they are detected primarily towards both filaments in B5. We postulate that these two components likely trace the two opposite sides in the line of sight of material contracting within the filaments;

• 3.

  We find that the simulations presented in Chen et al. (2020a) exploring filament formation scenarios do not accurately represent the filaments in B5. For both filaments in B5, we find C_v ≪ 1 across all the orthogonal cuts, and therefore, we rule out the turbulence- and flow-dominated formation scenarios for the two filaments. Both filaments have highly supercritical masses per unit of length, but do not show any strong orthogonal velocity gradient. This indicates a presence of additional support against gravity, likely via a poloidal magnetic field that is parallel to the spine of the filaments;

• 4.

  Assuming simple geometry and uniform density inside the core, we estimated the rate of infall of material onto the two filaments (seen with mid-red and mid-blue components) to be 1.6 × 10^–4 M_⊙ yr^–1 and 1.8 × 10^–4 M_⊙ yr^–1, respectively, for filament-1 and filament-2. Since the density at the surface, where the accretion likely takes place, is lower than the average density in the core, we note that these rates are upper-limits for the infall. The mass doubling time at these infall rates is 6–7 × 10⁴ yr. Therefore, the filament masses can change significantly during the typical core lifetime of 0.5–1 Myr;

• 5.

  We also obtain an upper limit of 3.5 × 10^–5 M_⊙ yr^–1 for the rate of infall of mass onto B5 (using the far-blue component) from the north. This rate is comparable to the rate of infall onto dense filaments in Serpens South. If more cores are observed to accrete new material after their formation, that would indicate that a significant update is required in the current core evolution models, where the cores are treated as evolving in isolation.

Our results highlight that infall at different scales within dense cores can be traced with high-sensitivity observation of suitable molecular lines using multi-component analysis. Observations of accretion of fresh new material onto dense cores during its evolution suggest that the current physical and chemical models of core evolution might need a significant update.

Acknowledgements

S.C., J.E.P., P.C. and M.T.V. acknowledge the support by the Max Planck Society.

Appendix A Maps of Bayes factor $K_b^a$

As mentioned in Sect. 3.1, we use the Bayes factor, $K_b^a$ (given by Eq. 3), to decide whether model a is favoured over model b. Figure A.1 show the maps of the Bayes K-factors calculated between no fit and a 1-component fit $\left(K_0^1\right)$, between 1- and 2-component fits $\left(K_1^2\right)$ and between 2- and 3-components $\left(K_2^3\right)$. Following the work of Sokolov et al. (2020), we consider a threshold of $\mathrm{ln\hspace{0.17em}}{K}_b^a\ge 5$ to select model a over model b.

Fig. A.1 Maps of the Bayes factors, $K_0^1$ (top left), $K_1^2$ (top right), and $K_2^3$ (bottom). The black-solid or white-solid contours show the two filaments and black-dashed or white-dashed contours show the three condensations in B5, respectively. The scale-bar and beam size are shown in the top-right and bottom-left corners, respectively.

Appendix B Clustering method used in component assignment

As mentioned in Sect. 3.2, in the pixels with more than one velocity component in the best fit model, we select the component that results in smooth maps of centroid velocity (v_LSR) and the brightness temperature (T_MB). In other words, our grouping is done to minimise sudden and sharp changes in these two parameters between nearby pixels. For this purpose, we define a weighted dimensionless parameter distance, d_p, between two pixels as $$d_p=\sqrt{{\left(w_v\times \frac{\text{Δ}{\text{v}}_{\text{LSR}}}{{\text{v}}_{\text{norm}}}\right)}^2+{\left(\left(1-w_v\right)\times \frac{\text{Δ}{\text{T}}_{\text{MB}}}{{\text{T}}_{\text{MB},\text{norm}}}\right)}^2},$$(B.1)

where Δv_LSR and ΔT_MB are the differences in the two parameters at the two pixels, and w_v is the fractional weight (0 < w_v < 1) applied to Δv_LSR. The separations in v_LSR and T_MB are normalised by v_norm and T_MB,norm, respectively, which are the maximum differences in the respective parameters in the entire region. We use v_norm = 1 km s^–1 and T_MB,norm = 7.5 K (calculated from one-component fit results). Starting with a pixel in the region with a one-component fit, we look at each of the neighbouring pixels, search for the component with the lowest d_p, and group that component to the one-component maps. We then move radially outward, and repeat this process for the neighbours of each new pixel, until the entire region is covered. Thus, we obtain the extended component, which we refer to as ‘global component’. We explore clustering with different relative weights in the parameters (w_v = 0.2–0.8) and different starting points, but find that it does not significantly change the maps of global component.

We then have one or two additional components in some pixels. Using different clustering methods as mentioned above, we are unable to get smooth parameter maps for these two components, which indicates the need for further subgroups. The additional components are then grouped according to their velocity, as described in Sect. 3.2.

Appendix C Line of sight velocities along orthogonal cuts through major axes of the filaments

Figs. C.1 and C.2 show the gradient in velocity across cuts orthogonal to the major axes of the filaments 1 and 2, respectively. These cuts are spaced 1 beam apart, and are perpendicular to the filament spine at that location.

Fig. C.1 Line of sight velocity profiles in the orthogonal cuts in filament 1. The red-dotted vertical lines show the position of the filament spine and the black-dotted vertical lines show the boundary of the filament for the respective orthogonal cuts shown. The Δv parameter (Δvh in Eq. 6) for the cut is shown in the upper-left corner. The Δvh value needed to have Cv = 1 $\left(\text{Δ}{\text{v}}_{{\text{C}}_{\text{V}}=1}\right)$ is shown in the bottom-right corner.

Fig. C.2 Line of sight velocity profiles in the orthogonal cuts in filament-2. The red-dotted vertical lines show the position of the filament spine and the black-dotted vertical lines show the boundary of the filament for the respective orthogonal cuts shown. The Δv parameter (Δvh in Eq. 6) for the cut is shown in the upper-left corner. The Δvh value needed to have Cv = 1 $\left(\text{Δ}{\text{v}}_{{\text{C}}_{\text{V}}=1}\right)$ is also shown.

Figure C.3 shows the Kernel Density Estimate (KDE) distribution of C_v values across the two filaments. This distribution highlights that the values have a strong distribution peak below C_v = 0.01, and a smaller secondary peak at ≈ 0.035.

Fig. C.3 Kernel density estimate of the distribution of the Cv values across both filaments.

Appendix D Comparison with filament formation scenarios considered in Chen et al. (2020a)

As mentioned in Sect. 4.2, there is a systematic and significant difference between the C_v values presented in this work (C_v ≪ 1) and those from Chen et al. (2020a) simulations (C_v ~ 0.3).

One key difference between the two works is in the mass-per-unit length (M/L) of the filaments. Chen et al. (2020a) considered M/L ~ 10 M_⊙ pc^–1 in their simulations, whereas for the filaments in B5, M/L ≈ 80 M_⊙ pc^–1 (Schmiedeke et al. 2021). The C_v ≪ 1 values for the B5 filaments could also imply that the filament formation in B5 was rather cylindrically symmetric (closer to the scenario depicted in the right part of Fig. 1 in Chen et al. 2020a), as opposed to being formed in a planar geometry, which was explored by their simulations (left part in the same figure).

Another possible difference could be the magnetic fields present in the filaments. As mentioned in Sect. 4.2, for the required support against gravity in the B5 filaments, the estimated magnetic field is poloidal and parallel to the spine of the filament, with a strength of ~ 500 µG. This would suggest that the strength and structure of the magnetic field is vastly different from that in the simulations of Chen et al. (2020a), where the magnetic fields are ~ 10 µG, and roughly perpendicular to the filament spine.

While some velocity profiles in the orthogonal cuts to the filaments (Figs. 7, C.1, and C.2) have similar shapes to those in simulations of Chen et al. (2020a), the absolute difference in velocities at the boundary in these cuts is still relatively small. Only 6 out of the 34 orthogonal cuts show a velocity gradient ≥ 0.1 km s^–1 at the filament boundary (black-dotted lines in Figs. 7, C.1, and C.2), and none of them exceed 0.15 km s^–1. In comparison, the velocity gradient seen in the filaments in the simulations of Chen et al. (2020a) is ≳ 0.15 km s^–1, except in the early stages. This could be another potential factor behind the difference in the C_v values between the two works.

References

All Figures

Fig. 1 Integrated intensity of NH3 (1,1) in Barnard 5. The black-solid and black-dashed contours show the two filaments and the three condensations in B5, respectively (contours and nomenclature adapted from Schmiedeke et al. 2021). The beam size and scale-bar are shown in the bottom-left and the top-right corners, respectively. The white-solid contour shows the core boundary calculated using single-dish observations by Pineda et al. (2010).

In the text

Fig. 2 Number of components detected towards different parts of B5, as decided by the Bayesian approach described in Sect. 3.1 (left). These components are then subdivided into 5 separate components according to their velocities (see Sect. 3.2). The black-solid and black-dashed contours show the two filaments and the three condensations in B5, respectively (contours adapted from Schmiedeke et al. 2021). The black-dotted contour shows the boundary of the coherent core in Pineda et al. (2010). The beam size and scale-bar are shown in the bottom-left and the top-right corners, respectively. Examples of NH3 (1,1) spectra from regions with one-, two-, and three-component fits (right, with the positions of the spectra indicated in the left panel). Only a small range in velocity is shown here to clearly highlight the individual components, which are illustrated in blue, red, and purple. In each panel, the fitted model is shown in green, while the residuals are shown in grey. Note: there are multiple, closely-spaced hyperfines in the velocity range shown here. The strongest two of these hyperfines can be seen in the narrow component (blue and purple in panels 2 and 3). However, these hyperfines are blended together in the broad component (green in panel 1, red in panels 2 and 3), and therefore the hyperfine structure is not visible.

In the text

Fig. 3 Velocity map of ‘global component’ (left). Velocity map of the four additional components (right). Each of these four components show a very narrow range of velocities. The median velocities of these components are 9.61 km s−1 (far-blue), 9.99 km s−1 (mid-blue), 10.32 km s−1 (mid-red), and 10.58 km s−1 (far-red). The black-solid and black-dashed contours show the two filaments and the three condensations in B5, respectively. The common colour scale is shown in the left panel. The black-dotted contours show the boundary of the coherent core in Pineda et al. (2010). The scale-bar and beam size are shown in the top-right and bottom-left corners, respectively.

In the text

Fig. 4 Velocity dispersion of the ‘global component’ (left). Velocity dispersions of the four additional components (right). The white-solid and white-dashed contours show the two filaments and the three condensations in B5, respectively. The white/black-dotted contours show the boundary of the coherent core in Pineda et al. (2010). The scale-bar is shown in the top-right corner and the beam size is shown in the bottom-left corner.

In the text

	Fig. 5 Schematic showing the spatial distribution of the different components in the line of sight. The coloured arrows show the direction of the velocity of the corresponding component. The ellipses represent the core and a filament inside the core, as indicated. The grey arrows indicate the north and south directions in the maps presented in this work.
In the text

	Fig. 6 Orthogonal cuts to the major axes of the two filaments. The cuts corresponding to the velocity gradients shown in Fig. 7 are shown with thick lines. The background greyscale shows the integrated intensity of NH3(1,1). The cuts are spaced one beam apart. The beam is shown in the bottom-left corner.
In the text

Fig. 7 Selected examples of the line of sight velocity profiles in the orthogonal cuts in filament-1 (top) and filament-2 (bottom). The red-dotted vertical lines show the position of the filament spine and the black-dotted vertical lines show the boundary of the filament for the respective orthogonal cuts shown. The Δv parameter (Δvh in Eq. (6)) for the cut is shown in the upper-left corner. The Δvh value needed to have Cv = 1 $\left(\text{Δ}{\text{v}}_{C_V=1}\right)$ is shown in the bottom-right corner for the upper panels, and below the Δvh in the lower panels.

In the text

	Fig. 8 Geometry considered to calculate the mass inſall onto filament-1 (top) and filament-2 (bottom).
In the text

	Fig. 9 Geometry considered for calculating the mass infall onto B5.
In the text

	Fig. A.1 Maps of the Bayes factors, $K_0^1$ (top left), $K_1^2$ (top right), and $K_2^3$ (bottom). The black-solid or white-solid contours show the two filaments and black-dashed or white-dashed contours show the three condensations in B5, respectively. The scale-bar and beam size are shown in the top-right and bottom-left corners, respectively.
In the text

Fig. C.1 Line of sight velocity profiles in the orthogonal cuts in filament 1. The red-dotted vertical lines show the position of the filament spine and the black-dotted vertical lines show the boundary of the filament for the respective orthogonal cuts shown. The Δv parameter (Δvh in Eq. 6) for the cut is shown in the upper-left corner. The Δvh value needed to have Cv = 1 $\left(\text{Δ}{\text{v}}_{{\text{C}}_{\text{V}}=1}\right)$ is shown in the bottom-right corner.

In the text

Fig. C.2 Line of sight velocity profiles in the orthogonal cuts in filament-2. The red-dotted vertical lines show the position of the filament spine and the black-dotted vertical lines show the boundary of the filament for the respective orthogonal cuts shown. The Δv parameter (Δvh in Eq. 6) for the cut is shown in the upper-left corner. The Δvh value needed to have Cv = 1 $\left(\text{Δ}{\text{v}}_{{\text{C}}_{\text{V}}=1}\right)$ is also shown.

In the text

	Fig. C.3 Kernel density estimate of the distribution of the Cv values across both filaments.
In the text