© The Authors 2023

1 Introduction

Filamentary structures in molecular clouds (e.g., Ward-Thompson et al. 2010, for the Polaris Flare cloud, Kirk et al. 2013, for the Taurus region, and Alves de Oliveira et al. 2014, for the Chamaeleon cloud complex) have been revealed by Herschel observations (André et al. 2010; Molinari et al. 2010), suggesting that pre-stellar cores may form from the gravitational fragmentation of marginally supercritical and magnetized filaments (Könyves et al. 2015; Benedettini et al. 2018). The role of the magnetic field in the star formation process, and particularly of its orientation with respect to the filament axis, has been widely discussed (Soler et al. 2013; Planck Collaboration Int. XXXIII. 2016) and it represents a crucial point in understanding how the cores embedded in the filament grow in mass and trigger the formation of protostellar objects. Observational data from Zeeman-effect surveys (Crutcher 2012) show that the maximum strength of the interstellar magnetic field is ~10 μG for gas with densities below n_H ~ 300 cm⁻³, increasing significantly at higher densities. This increase is determined by the balance between ambipolar diffusion and the accumulation of the magnetic field due to condensation. The comparison between magnetized and non-magnetized models, either in equilibrium (Toci & Galli 2015a,b) or dynamical (Hennebelle & Audit 2008), shows that magnetic fields have a strong influence on the filamentary structure and the fragmentation process within the filaments (Tilley & Pudritz 2007; Kirk et al. 2015). A relevant driver of these processes is the geometry of the magnetic field.

Polarization measurements often show that the orientation of the magnetic field is nearly perpendicular to the major axis of the star-forming filament, but aligned with low-density filamentary structures (e.g., Soler et al. 2013; Planck Collaboration Int. XXXII. 2016; Jow et al. 2018; Soler 2019), a result also found in numerical simulations (e.g., Nakamura & Li 2008; Soler & Hennebelle 2017; Seifried et al. 2020; Körtgen & Soler 2020). Cox et al. (2016) analyzed the filamentary structure of the Musca cloud, differentiating the main filament and the low-density filamentary structures close to it, showing that they are parallel to the plane-of-the-sky local magnetic field and quasi-perpendicular to the main filament. Additional observations revealed a large-scale network of sub-filaments connected to the filamentary structures (e.g., Schneider et al. 2010), aligned with the direction of the magnetic field, which can represent a mass reservoir for further growth of the filament and the cores embedded into it.

The accretion of filaments is important in the context of star formation, because it could play a key role in preventing dense filaments from collapsing to spindles, and in maintaining constant widths during the evolutionary process (André 2017). Some authors have focused on different ways in which the ambient material could be accreted, while considering the constraints imposed by the magnetic field direction inside molecular clouds. Based on Herschel observations, Shimajiri et al. (2019) showed that the accretion of ambient gas can be driven by the gravitational potential of the filament. This provides a strong support for the scenario of mass accretion along magnetic field lines (oriented nearly perpendicular to the filament) proposed by Palmeirim et al. (2013). Assuming a magnetic field perpendicular to the filament, Hennebelle & André (2013) considered radial accretion in a self-gravitating filament due to the combination of accretion-driven turbulence and dissipation of the turbulence by ion-neutral friction. Gómez et al. (2022) proposed that the width measured by observations may change depending on the tracer used and then suggested a tracer-dependent estimate of the accretion rate onto the filament. In Gómez & Vázquez-Semadeni (2014), a similar behavior was shown, whereby filaments were proposed to be river-like flow structures, with the gas falling onto the filament changing direction as the gas density increases, and accreting mainly in the perpendicular direction.

Magnetic fields also have an impact on the characteristic width of the filament. From an analysis of Herschel observations, Arzoumanian et al. (2011) identified a peaked distribution of filament widths around 0.1 pc, including both low-column density and star-forming filaments. Whether or not this width is a universal filament property is still a matter of debate (e.g., Smith et al. 2014; Panopoulou et al. 2017; Hacar et al. 2022), but several theoretical works have attempted to explain it. For example, in simulations of idealized filaments, Seifried & Walch (2015) found that the 0.1 pc width (Arzoumanian et al. 2011) assumed in their initial conditions could be maintained if the magnetic field is longitudinal.

Ideal magnetohydrodynamics (MHD) may not be sufficient for studying filaments, due to the importance of ambipolar diffusion and other non-ideal processes (Grassi et al. 2019). The coupling between magnetic field and gas depends on the density and the ionization degree (Shu 1983). The resistivity coefficient of ambipolar diffusion is indeed strongly dependent on the abundance of ions (e.g., Grassi et al. 2019), with collisions between charged dust grains and neutral gas particles dominating the momentum transfer. In this context, cosmic rays play a fundamental role as they ionize molecular clouds (Padovani et al. 2009; Padovani & Galli 2011). Chen & Ostriker (2014) simulated core formation in colliding flows including ambipolar diffusion, finding a transient ambipolar diffusion phase during the formation of the shock layer that allows for the formation of cores with higher mass-to-flux ratio. In a subsequent work, Ntormousi et al. (2016) reported wider filaments as a consequence of non-ideal effects. However, numerical simulations including ambipolar diffusion in MHD turbulence have not determined that it plays a role in setting a characteristic scale (Oishi & Mac Low 2006; Burkhart et al. 2015).

In this paper, we explore how different physical parameters, such as the magnetic field strength and cosmic ray ionization rate, affect the evolution and accretion of a self-gravitating filament, including chemistry and non-ideal effects. In Sect. 2, we introduce our initial conditions and the non-ideal magnetohydrodynamic equations, together with details on the chemistry and microphysics employed in this work. In Sect. 3, we discuss the reference model and the parameter study, along with the impact of those parameters on the accretion process. In Sect. 4, we present a discussion of the main results referring to the accretion rates and filament widths. Finally, in Sect. 5, we discuss some limitations of our approach and in Sect. 6, we summarize our main conclusions.

2 Methodology

In the following, we describe the initial conditions and the theoretical framework adopted in this work.

2.1 Initial conditions

Our model consists of a self-gravitating turbulent filament with 0.1 pc width in a 2 pc size periodic box sampled with 1024 cells, which can accrete mass from both sides of the x-axis, namely, the coordinate on which the hydrodynamical equations are solved. For our study, we model the filament as a nonuniform slab with an oblique magnetic field.

To model the filament, we adopted a Plummer-like profile as shown in Arzoumanian et al. (2011), given by: $$\rho \left(x\right)=\frac{{\rho }_{\text{ridge}}}{{\left[1+{\left(x/x_{\text{flat}}\right)}^2\right]}^{\rho /2}},$$(1)

where ρ_ridge is the central density of the filament, x_flat is the characteristic width of the flat inner plateau of the filament, and p is the typical exponent of the profile. For this study, we assumed the same parameters as adopted in Körtgen et al. (2018), namely, p = 2 and x_flat = 0.033 pc, and we aligned the minor axis of the filament with the x-axis. Figure 1 shows a sketch of the filament setup.

The box is initialized with a uniform temperature of 15 K and is filled with a low-density medium with ρ₀ ~ 2 × 10⁻²¹ g cm⁻³ outside the filament region, with the latter modeled according to the Plummer-like profile in Eq. (1), where we set ρ_ridge ~ 3 × 10⁻¹⁹ g cm⁻³. The mass per unit length ratio is calculated as (M/L) = 3(M/L)_crit with ${\left(M\text{/}L\right)}_{\text{crit}}=2c_{\text{s}}^2\text{/}G$ (Inutsuka & Miyama 1992) and $c_{\text{s}}=\sqrt{k_{\text{B}}T\text{/μ}{m}_{\text{p}}}$ is the speed of sound, with G, k_B, μ and m_p being the gravitational constant, the Boltz-mann constant, the mean molecular weight, and the mass of the proton, respectively. Then, (M/L) ~ 86 M_⊙ pc⁻¹. We note that we initialize the filament at the critical line mass to guarantee an axial collapse.

The free-fall time of the background gas, defined as: $$t_{\text{ff}}=\sqrt{\frac{3\pi }{32G{\rho }_0}},$$(2)

corresponds to ~1.3 Myr.

The filament is initialized including turbulence with Mach number ℳ = 2 in the υ_x, υ_y, and υ_z components of the velocity field, following a Burgers-like power spectrum, as described in Bovino et al. (2019; see also Körtgen et al. 2017). Given the randomness of the initial velocity field, to avoid any net momentum that could produce a drift in the filament along the direction of the gravitational force, we subtract the mean initial velocity of the x-component.

Fig. 1 Representative sketch of our set-up. The code evolves a vector variable with three components (x, y, z,) along one coordinate (x); hence, y and z have periodic boundary conditions. The x coordinate traverses the filament (shaded gradient) and the initial density distribution along this component follows a Plummer-like profile with maximum density, ρridge, and background density, ρ0, see Eq. (1). The gas accretes along the same axis, as indicated by the arrows. The magnetic field is inclined by an angle θ with respect to the x-axis, where θ = arctan By/Bx. We note that the color palette employed for the filament is only included for illustration purposes.

2.2 Magnetic field orientation

The magnetic field B is set perpendicular to the z-axis and is tilted relative to the x-axis (the accretion flow direction) by an angle θ = arctan(B_y/B_x), as shown in Fig. 1. We note that when θ = 0 the field is parallel to the flow direction, while for θ = π/2 the field is perpendicular to the flow. Given the dimensionality of the code and the implicit assumption that ∂_y = ∂_z = 0, the solenoidal condition requires ∂_xB_x = 0. Hence, B_x remains constant throughout the box, whereas the initial condition of the y component of the field is assumed to scale as (Crutcher et al. 2010): $$B_y=B_{y,0}{\left(\frac{\rho }{{\rho }_0}\right)}^k,$$(3)

where k = 0.5, and B_y,0 ≡ B₀ sin θ, with B₀ = 10 μG. We note that since B_x is constant by construction, θ also varies with ρ, indicating that θ in our initial conditions refers to B_y = B_y,0.

2.3 Numerical framework

To study the evolution of the filament, we employ the LEMON-GRAB code (Grassi et al. 2019), a 1D code that solves the nonideal MHD equations time-dependently. The code also includes cooling and heating processes as well as grain chemistry (including charged grains). For the purposes of this study, we added self-gravity to LEMONGRAB, as described below.

The coupled non-ideal MHD equations can be written as: $${\partial }_t\rho =-{\partial }_x\left[\rho {\upsilon }_x\right],$$(4) $${\partial }_t\left[\rho {\upsilon }_x\right]=-{\partial }_x\left[\rho {\upsilon }_x^2+P^{\ast }-\frac{B_x^2}{4\pi }\right]-\rho {\partial }_x\text{Φ},$$(5) $${\partial }_t\left[\rho {\upsilon }_y\right]=-{\partial }_x\left[\rho {\upsilon }_x{\upsilon }_y-\frac{B_x{B}_y}{4\pi }\right],$$(6) $${\partial }_t\left[\rho {\upsilon }_z\right]=-{\partial }_x\left[\rho {\upsilon }_x{\upsilon }_z-\frac{B_x{B}_z}{4\pi }\right],$$(7) $${\partial }_t{B}_x=0,$$(8) $${\partial }_t{B}_y=-{\partial }_x\left[{\upsilon }_x{B}_y-{\upsilon }_y{B}_x+\frac{{\eta }_{\text{AD}}}{{\left|B\right|}^2}(F_{B,x}{B}_y-F_{B,y}{B}_x)\right],$$(9) $${\partial }_t{B}_z=-{\partial }_x\left[{\upsilon }_x{B}_z-{\upsilon }_z{B}_x+\frac{{\eta }_{\text{AD}}}{{\left|B\right|}^2}(F_{B,x}{B}_x-F_{B,x}{B}_z)\right],$$(10) $$\begin{array}{c}{\partial }_tℰ={\text{Γ}}_{\text{cr}}-\wedge -\rho {\upsilon }_x{\partial }_x\text{Φ}-{\partial }_x\{\left(ℰ+P^{\ast }\right){\upsilon }_x-\frac{B_x}{4\pi }\left(V·B\right)\\ -\frac{{\eta }_{\text{AD}}}{\pi {\left|B\right|}^2}\left[\left(F_{B,z}{B}_x-F_{B,x}{B}_z\right)B_z-\left(F_{B,x}{B}_y-F_{B,y}{B}_x\right)B_y\right]\},\end{array}$$(11) $${\partial }_t\left(\right(\rho X_i\left)\right)=-{\partial }_x\left(\right(\rho X_i{\upsilon }_x\left)\right)+{𝒫}_i-\rho X_i{\mathcal{L}}_i.$$(12)

In the above equations, Φ is the gravitational potential, Λ is the cooling rate, Γ_cr is the cosmic ray heating rate, and 𝒫_i and ℒ_i are the production and loss terms of the ith chemical species.

The total pressure P* is: $$P^{\star }=P+\frac{{\left|\right|B\left|\right|}^2}{8\pi }.$$(13)

We assume an ideal equation of state for the thermal pressure: $$ℰ=\frac{P}{\gamma -1}+\frac{\rho {\left|\right|v\left|\right|}^2}{2}+\frac{{\left|\right|B\left|\right|}^2}{8\pi },$$(14)

and related to the temperature T, needed by the chemistry, via the ideal gas law $$P=\frac{\rho k_B}{\text{μ}{m}_{\text{p}}}T.$$(15)

The components of the Lorentz force are $$F_{B,x}=-B_y·{\partial }_x{B}_y-B_z·{\partial }_x{B}_z,$$(16) $$F_{B,y}=B_x·{\partial }_x{B}_y,$$(17) $$F_{B,z}=B_x·{\partial }_x{B}_z,$$(18)

The ambipolar diffusion resistivity coefficient is $${\eta }_{\text{AD}}=c^2\left(\right(\frac{{\sigma }_{\text{p}}}{{\sigma }_{\text{p}}^2+{\sigma }_{\text{H}}^2}-\frac{1}{{\sigma }_{\left|\right|}}\left)\right),$$(19)

where σ_Ρ, σ_Η, and σ_‖ are, respectively, the Pedersen, Hall, and parallel conductivities, and c is the speed of light.

Given the spatial discretization of the code, the Poisson equation can be easily discretized over the 1D grid at second-order via central differencing, as $${\partial }_x^2\text{Φ}\approx \frac{{\text{Φ}}_{k+1}-{\text{Φ}}_k+{\text{Φ}}_{k-1}}{△x^2}=4\pi G_{\rho k},$$(20)

where Φ_k and ρ_k are the gravitational potential and density of the kth cell (ranging from 1 to N, with Ν as the number of resolution elements in the simulation), and Δx the spatial resolution of the simulation. For the boundary cells, we have included different boundary conditions (periodic, Dirichlet, or outflowing) among which the user can arbitrarily choose; although in this work we only consider periodic ones, that is, Φ₁ = Φ_Ν. The potential is then calculated via matrix inversion, and the gravitational acceleration is finally obtained via central differencing as: $$g_x=-{\partial }_x\text{Φ}\approx -\frac{{\text{Φ}}_{k+1}-{\text{Φ}}_{k-1}}{2△x}.$$(21)

In this work, we employ a cooling function which depends on both temperature and number density. It corresponds to conditions of collisional ionization equilibrium typical of the interstellar gas, determined using CLOUDY (Ferland et al. 1998) and tabulated by Shen et al. (2013). We imposed a temperature floor of 10 K by including an artificial heating term, defined as Λ(Τ = 10 K), and added to the ∂_tℰ equation.

For details of the calculation of pressure, heating, cosmic ray ionization rate, conductivities, and the chemical evolution of the species, we refer to Grassi et al. (2019). In this work, we employed their reduced chemical network, including eight species: electrons e⁻, X, X⁺, and neutral and charged grains, namely: g, g⁺, g⁺⁺, g⁻, g⁻⁻; The species X and X⁺ are a proxy of all neutrals and cations produced by a chain of fast reactions following H₂ ionization (see Fujii et al. 2011; Grassi et al. 2019, for details). We assumed that the mass of X is the same as the mass of H₂, and the dust grains rates were weighted assuming a MRN distribution (Mathis et al. 1977) with the same characteristics as in Grassi et al. (2019). The cosmic ray ionization rate of X is assumed to be initially uniform at a value ζ_cr = 10⁻¹⁷ s⁻¹. We set the initial electron fraction to $f_i=n_{{\text{e}}^{-}}\text{/}{n}_{{\text{H}}_2}={10}^{-7}$ and the dust-to-gas mass ratio to 𝒟 = 10⁻², as listed in Table 1.

3 Results

In this section, we present the results of our simulations divided in four parts. In Sect. 3.1, we present our reference model based on the initial conditions listed in Table 1. In Sect. 3.2, we show the effects of varying the orientation and strength of the magnetic field. Separately, in Sect. 3.3, we consider different cosmic rays ionization rate and, finally, in Sect. 3.4 we explore the effect of other physical quantities such as the initial turbulence seed, Mach number, density regimes, and k exponent of the magnetic field-density relation, namely, Eq. (3).

3.1 Reference case

In Fig. 2, we show, for the ideal and the non-ideal cases, the evolution of density, temperature, velocity in the x-direction, total velocity, the y-component of the magnetic field, total magnetic field, Mach number, mass flux, ambipolar diffusion resistivity coefficient, and electron fraction as a function of the x-coordinate at three different times, for our reference case (see Table 1). The same panels include a comparison with the ideal MHD case (dashed lines).

At early times, the density shows a peak with an initial maximum of about 3 × 10⁻¹⁹ g cm⁻³ (9 × 10⁴ cm⁻³) in the center, slightly increasing at later times. The peak value of the y-component and the total magnetic field is higher than 100 μG, triggered by the increase in density. In the outer regions, at 0.1t_ff the density is about 2 × 10⁻²¹ g cm⁻³ (6 × 10³ cm⁻³). The initial turbulence affects the first stages of the simulation generating fluctuations in the magnetic field, velocity, and temperature, the latter showing values of 20 K. The density and velocity profiles both indicate that a shock front is generated at around 0.05–0.1 pc from the filament center, after the initial stages.

From the velocity field, we calculated the Mach number as ℳ = υ/c_s, where $\upsilon ={\left({\upsilon }_x^2+{\upsilon }_y^2+{\upsilon }_z^2\right)}^{1/2}$ is the total velocity magnitude, and c_s the current sound speed of each cell. Similarly to the velocity field, the Mach number shows strong fluctuations at the first stages of evolution, going from values around 5 (highly supersonic) in the outer regions to ℳ < 1 (subsonic) in the center. The same behavior is seen at later times. From the density and the velocity field, we estimate the mass flux, reflecting the accretion flow in the x-direction: $$\dot{\text{Σ}}=\rho {\upsilon }_x.$$(22)

At early times, the mass flux shows peaks of ~10³ M_⊙ Myr⁻¹ pc⁻² as a result of the density increase, which remains high also in the subsequent time-steps with a slight decrease.

The gravitational acceleration due to the mass accumulated in the filament leads to peak velocities up to 2 km s⁻¹ moving inward, consistent with the expected free-fall due to the available mass. The y- and z-components of the velocity (not shown), exhibit a similar behavior to the x-component, but reach peak values of 1 km s⁻¹ and 0.5 km s⁻¹, respectively, with the latter close to the edges of the simulation box. This relatively small difference in magnitude is caused by the initial statistics of the turbulence, that has the same dispersion on the three spatial components, but zero average only on the x-component.

The thermal profile, after an initial isothermal state, shows three distinct regions: (i) a cold background gas at around 10 K, (ii) an efficiently heated region with a temperature that is higher than 100 K produced by gas shocking, and (iii) the filament ridge where the high density leads to efficient cooling (i.e., cooling time shorter than the dynamical time). This brings the gas down to 10 K.

The x-component of the magnetic field remains constant in space and time by construction, while the z-component (not shown) is initially zero because of the alignment followed by fluctuations of a few μG due to the interaction with turbulence and velocity fluctuations. The y-component dominates the evolution of the total magnetic field, as seen in panels (ƒ) and (h) of Fig. 2.

We note some differences from the ideal MHD case. In particular, over time, the density peak tends to be broader in the ideal MHD case, due to the increased magnetic pressure and stronger coupling between the magnetic field and the density.

During the evolution, the density profile further steepens. At low densities, namely, ρ ≲ 5 × 10⁻²⁰ g cm⁻³, the electron fraction (panel k) remains above 10⁻⁶, but as the density approaches a maximum value of about 10⁻¹⁸ g cm⁻³, the electron fraction drops to ~10⁻¹⁰ due to efficient recombination. The decrease observed in the ambipolar diffusion at high densities (n_H ~ 10⁵ cm⁻³) is determined by the electrons and X⁺, which are the dominant charge carriers, while grains take over at n_H > 10⁸ cm⁻³ (Marchand et al. 2016). The negatively charged grain abundance (not shown) is strictly linked to the electrons due to recombination; thus, the evolution of the two species is highly correlated. When the electron abundance drops off, the abundance of neutral and positively charged grains increases at high densities. Ambipolar diffusion keeps decreasing until the ion abundance increases enough for the diffusion with respect to the neutrals to become significant.

Fig. 2 Profile along the x-axis of various physical quantities of the model at three different times: t/tff = 0.1, 0.5 and 1.0. (a) density ρ, (b) temperature T, (c) x-component of the velocity υx, (d) total velocity magnitude |v|, (e) y-component of magnetic field By, (f) Mach number ℳ, (g) total magnetic field |B| magnitude, (h) mass flux $\dot{\text{Σ}}$, (i) ambipolar diffusion resistivity coefficient ηAD, and (j) abundance e−/H2. The solid lines correspond to non-ideal MHD, dashed lines to ideal MHD. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments. We note that in panel (j), in the ideal case (dashed lines), ηAD is reported for the sake of comparison, but it is not included in the MHD equations.

3.2 Magnetic field geometry and strength

To study how the initial geometry of the magnetic field influences the evolution of the filament, we tested four different initial inclinations θ with respect to the x-axis, corresponding to θ = 5° (almost parallel to the flow), θ = 45°, θ = 77° (reference case), and θ = 90° (perpendicular to the flow), while the other parameters are kept unchanged. In Figs. 3 and 4, the density, velocity and orientation angle of the magnetic field are shown at different times as a function of position. From now, we present the evolution of the model at two representative times: at half a free-fall time 0.5t_ff, and after one free-fall time 1.0t_ff. The central peak of the density and its width do not change very significantly for the different angles; However, when the magnetic field is nearly parallel to the gas flow, we would expect that the gas flow is strongly impeded by the magnetic field. As a result, the density in the outer parts becomes even lower when θ is small, without significantly changing the density peak in the center, as the extra mass that is added to the center is very small. Similarly, the peak velocity increases as θ becomes aligned to the flow.

In general, the orientation of the magnetic field evolves both as a function of position and as a function of time. In particular, in the outer parts of the filament, the orientation angle decreases as a function of time and the magnetic field tends to align with the flow of the gas. In the filament itself, on the other hand, the velocities are reduced due to increased thermal pressure and the perpendicular component of the magnetic field is compressed and amplified as more gas is accreted onto the filament. As a result, the orientation becomes closer to perpendicular within the dense component of the gas.

It is important to note that because of the assumed symmetry, the simplified setup here may not necessarily hold in a realistic setting. Nevertheless, what we can see clearly from Fig. 4 is the change in the angle of the magnetic field between the outer and the inner parts of the filament. In a realistic setting, such a change of the angle may for instance correspond to a “U”-shaped magnetic field, as suggested by Gómez et al. (2018), or potentially other types of configurations that imply a change of the angle as a function of scale.

To study the effect of the magnetic field strength, we determined how the structure of the filament depends on four different initial field strengths, corresponding to no field (numerically ~10⁻²⁵ μG), 5 μG (weak field), 10 μG (reference case), and 20 μG (strong field), while the other parameters were kept unchanged. We calculated the normalized mass-to-flux ratio (M/Φ_B)/(M/Φ_B)_crit at the initial condition to show the magnetic supercritical nature of our simulations: ~141 for 5 μG, ~62 for 10 μG, and ~31 for 20 μG. In Fig. 5, the density and the x-velocity are shown as a function of position. The magnetic field strength is one of the parameters that mostly affect the evolution of the filament, as highlighted after 0.5t_ff: as the magnetic field strength increases, the density peak of the filament appears smaller, while the filament appears wider. Because of the limited mass available, this also results in a background density that decreases (on average) for stronger magnetic fields (we note that because of the initial turbulent velocity, the local density in some cells does not necessarily reflect this trend). A simple explanation for this behavior is the increasing magnetic pressure, which counteracts gravity, slowing down accretion onto the filament, as can be also seen in the velocity profile, where the peak velocity is higher for weaker magnetic fields. This suppression results in higher densities in the background region, which also correspond to a higher average value of the magnetic field (despite the non-ideal effects) and lower inflow velocities (which arise from the combined effect of a smaller mass concentration producing a weaker gravitational pull and a larger pressure gradient counteracting it). At later times, these differences are almost washed out, since gravity dominates the system evolution, finally leading to an efficient accretion of material onto the central overdensity, as shown in the bottom panels. Nonetheless, mild differences in the density and velocity profiles remain in the underdense region outside the filament.

The density profile shows a similar trend at the earlier stages, whereas the velocity exhibits a peculiar evolution. In particular, the weakest magnetized case (5 μG) now has the lowest velocities. This result is due to the inversion from a magnetically dominated pressure to a thermally-dominated regime, where the previously higher inflow velocities produced a much stronger shock heating of the gas near the filament, which then developed into a strong pressure gradient suppressing the inflow. In any case, these differences remain small, especially considering the corresponding gas densities.

Fig. 3 Density and velocity along the x-axis for different initial angles. The solid line corresponds to the reference case. At 0.5tff in the top panels and at 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.

Fig. 4 Inclination of the magnetic field and y-component for the different initial angles. The solid line corresponds to the reference case. We have t = 0.5tff in the top panels and t = 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.

Fig. 5 Density and velocity along the x-axis for the different initial magnitudes of the magnetic field. The solid line correspond to the reference case. We have t = 0.5tff in the top panels and t = 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.

3.3 Cosmic ray ionization rate

To test the effect of the cosmic ray ionization rate (CRIR) on the evolution of the filament, we performed three simulations with ζ_cr = 10⁻¹⁸ s⁻¹ (low), 10⁻¹⁷ s⁻¹ (reference), and 10⁻¹⁶ s⁻¹ (high), respectively, while the other parameters were kept unchanged. The density, velocity, electron abundance, and ambipolar diffusion resistivity coefficient are shown in Figs. 6 and 7 for the different cases and at two different times. Higher cosmic-ray ionization rates produce an increase in the ionization fraction in the low-density accretion region (see Fig. 7, which corresponds to a stronger coupling between the gas and the magnetic field). As shown in Fig. 7, the lower electron density (left panels) corresponds to a less-ideal MHD, namely, a higher ambipolar diffusion coefficient (right panels), determining a direct relation between CRIR and non-ideal behavior. Analogously, the density profile shows a broader distribution caused by a slow accretion toward the filament center, as shown in the right top panels of Fig. 6. This behavior is also confirmed by the velocity profile at 0.5t_ff, which shows in the high CRIR case a shift of the υ_x peak toward higher radii, determined by the relatively higher magnetic pressure support. This peak is no longer present at 1.0t_ff, where the magnetic pressure halts the gravitational collapse, producing the slowest infall in the high CRIR case.

3.4 Other parameters

To complete our analysis, we studied the effects of additional parameters on the evolution of the filament. We have explored independently different initial Mach numbers, initial central density, the random seeds for the turbulence, and the exponent of the magnetic field-density relation, as in Eq. (3), while the other parameters were kept unchanged.

By changing the initial turbulence seed, we did not observe relevant differences in the final evolution of the filament, suggesting that in our case the global features are dominating over the local turbulence, which dissipates over time. In fact, the turbulence-induced fluctuations (e.g., in the density profile) are more prominent at earlier stages.

For the different density regimes, we reduced ρ_ridge by a factor of 10 and 2, with the free-fall time increasing to ~4.3 Myr and ~1.9 Myr, respectively. As expected, at lower density it becomes harder for the filament to accrete material and there is no real evolution over time, with the density keeping around the ambient values. Similarly, reducing the density by a factor of 2 produces broader density profiles and a slower evolution. By changing the exponent of the magnetic field – density relation from k = 0.5 to k = 0.4, Eq. (3), the filament reaches higher temperatures, velocity, and central density, but almost the same density peak with respect to the reference, since the lower magnetic pressure cannot slow down the accretion efficiently. Conversely, for k = 0.6 (cf. 0.65 in Crutcher et al. 2010) the magnetic field is higher in both the ambient and the filament and it becomes more perpendicular to the flow compared to the previous cases, resulting in a slower accretion. Although these parameters affect the global evolution, their effect is less pronounced than the effect of CRIR and of the initial magnetic field geometry.

Fig. 6 Density and velocity on the x-axis for different cosmic ray ionization rates. The black line correspond to the reference case. We have: t = 0.5tff in the top panels and t = 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.

Fig. 7 Ionization fraction and ambipolar diffusion resistivity coefficient for the different cosmic ray ionization rates. The black line correspond to the reference case. We have t = 0.5tff in the top panels and t=1.0 tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.

4 Discussion

In order to assess the impact of the different parameters during the evolution of the gas, we should define one or more metrics that describe quantitatively the shape and the global morphological characteristics of the filament. The formal definition of its shape and width has been discussed by several authors, as we report in Appendix A. Regarding the aims of our work, we found that our simulation results are better interpretable by defining two global metrics: (i) the average accretion mass flux, $\langle \dot{\text{Σ}}\rangle$, and (ii) the full-width-at-half-maximum, FWHM, obtained by fitting the density profile with a Gaussian function (see e.g., Arzoumanian et al. 2011) including the background up to 0.2 pc for a more suitable fitting of the curve. The former is averaged over time, while the second is calculated at 1.0t_ff, that is, at the end of the simulation. A summary of these two metrics and their relative ratio with respect to the reference case is reported in Table 2.

As both metrics are controlled by the gas pressure support (either thermal or magnetic), the accretion mass flux and the FHWM appear to be anticorrelated: the FWHM increases for higher pressure, while the accretion rate tends to be reduced. This can be explained by observing the ideal MHD case, where the contribution of the magnetic pressure to the total pressure, reduces both the accretion and the capability of the filament to reach narrower configurations, resulting in larger FWHMs. In fact, in this case, the FHWM is 7 times larger than the reference case. For the same reason, the CRIR cases present an increased FWHM in the high-CRIR model, and vice versa for the low-CRIR. The magnetic pressure plays a crucial role also when changing the initial magnetic field; reducing the field strength (no-B and 5 μG) narrows the density profile, while increasing the magnetic strength to 20 μG causes a relatively broader FWHM.

By changing the orientation of the initial magnetic field (θ), we were able to observe a narrower FWHM when B_y tends to be parallel to the flow (i.e., θ = 5°, almost perpendicular to the filament ridge). This decrease is driven by the magnitude of B_y, which is the only component that (by construction) varies spatially¹. When θ is small, B_x is the dominant component of the magnetic field, namely, it is the component that (having no spatial gradient) plays a minor role in the MHD equations. In fact, in the right panels of Fig. 3, υ_x is larger for the smallest θ, being dominated by self-gravity and resulting in higher accretion and a narrower FWHM.

To assess the impact of turbulence on the filament properties, we show in Fig. 8, the time evolution of the FWHM for different values of the Mach number for the reference case. High Mach-number cases tend to behave similarly to ℳ = 2, except at earlier times when turbulence is dominant.

5 Limitations of the model

As in any numerical model, some approximations and assumptions need to be introduced in order to make the problem computationally tractable when including detailed microphysics and when a parameter study is planned. The main limitation in our approach is the 1D approximation employed to model the filament along the x coordinate. Despite the code having three components for the vector variables, such as velocity and magnetic field, the y and z coordinates have periodic boundary conditions that limit the exploration of the spatial variability. Additionally, in a 1D code the solenoidal condition imposes by construction that the x component of the magnetic field is constant in time and space. Another set of periodic boundary conditions is imposed at the boundaries of the simulation domain to avoid infinitely growing accretion onto the filament (outflow boundaries), or nonphysical configurations (Dirichlet boundary conditions, i.e., zero derivatives). Ideally, this could be avoided by using a much larger simulation box, where the boundaries remain unaffected during the simulation due to their distance from the central filament ridge. However, such a set-up requires much more computational resources or an adaptive mesh, which is beyond the aims of the present work.

Although the microphysics is limited by the reduced set of chemical reactions, it is still nevertheless capable of capturing the main features of larger chemical networks, especially the influence on the non-ideal behavior driven by ambipolar diffusion, as discussed in Grassi et al. (2019). In addition, since the chemistry has a noticeable influence on the cooling functions, the current network might reduce the capability of our model to determine time-dependent effects of the thermal evolution. Finally, since our reduced chemical network cannot be used to determine the ratio between the atomic and molecular species accurately, we use a constant molecular weight and adiabatic index that might somewhat affect the results. Additional limitations and method details are discussed in Grassi et al. (2019), where the main features of the code are also presented.

Fig. 8 Time evolution of the FWHM ratio of ℳ = 4 (dot-dashed line) and ℳ = 6 (dashed line) with ℳ = 2 (reference).

6 Summary and conclusions

In this 1D study, we modeled an accreting self-gravitating turbulent filament using LEMONGRAB, a non-ideal MHD code that includes chemistry and microphysical processes. We explored how the main parameters, namely, the configuration and strengths of the magnetic field and the cosmic ray ionization rate, as well as other parameters such as the Mach number, initial turbulence level, initial ridge density, and the exponent of the magnetic field-density relation affect the evolution and accretion of the filament. Our main results can be summarized as follows:

1. Including non-ideal MHD is crucial with regard to the filament evolution. The ideal case produces a wider filament, with an FWHM that is approximately seven times greater than the non-ideal model, due to the increased magnetic pressure support. This suggests that non-ideal MHD is fundamental to understanding the accretion process;

2. Independently of the initial configuration of the magnetic field, the magnetic field lines in the central part of the filament bend to reach a perpendicular configuration with respect to the x-axis (i.e., the accretion direction). This is consistent, for example, with “U”-shaped model of Gómez et al. (2018);

3. Higher cosmic ray ionization rates produce higher ionization degrees that correspond to a more ideal gas. As the coupling with the magnetic field is stronger, the magnetic pressure halts the collapse, resulting in a broader filament;

4. We did not find any strong change in the final FWHM value with the Mach number, pointing to a less pronounced role of the turbulence in regulating the evolution and the final properties of the filament compared to microphysics and the magnetic pressure. However, a one-dimensional model could be inadequate to address this specific question.

We conclude that adding magnetic fields and non-ideal MHD effects significantly affects the evolution of a collapsing filament, its width, and its accretion rate, while other parameters play only a minor role. Special attention should be given to the cosmic ray ionization rate that strongly affect the coupling between the gas and the magnetic field. It is then fundamental for future works to include a proper cosmic rays propagation scheme to accurately study their effect on the accretion rate and on the filament width.

Acknowledgements

We are grateful to the referee for the detailed report. NGV acknowledge support by the National Agency for Research and Development (ANID)/Scholarship Program / Magíster Nacional/2020-22200413. DRGS thanks for funding via Fondecyt Regular (project code 1201280). S.B. and D.S. gratefully acknowledge support by the ANID BASAL projects ACE210002 and FB210003. A.L. acknowledges funding from MIUR under the grant PRIN 2017-MB8AEZ. T.G. acknowledges the financial support of the Max Planck Society.

Appendix A Boundaries of the filament

As there is no unique definition of the boundaries of a filament, we calculated the mass flux for the reference case at different distances from the filament axis to make a comparison with the observational data. These positions are given by: (a) the minimum values of the velocity-divergence (or x-velocity gradient), as in Priestley & Whitworth (2022); (b) x = ±0.05 pc from the axis of the filament (half of the filament width reported by Arzoumanian et al. 2011); (c) the distance where the enclosed mass is 80% of the total mass in the computational box; and (d) where the flow becomes supersonic with a possible shock with Mach number ℳ = 1. Fig. A.1 shows the evolution of the mass fluxes in each time-step for each case, where dots correspond to non-ideal MHD and stars to ideal MHD. By looking at the ideal case (stars) we note that the difference is not very large after 0.3t_ff. The total averaged mass fluxes in time for non-ideal MHD are 206.15 M_⊙ Myr⁻¹ pc⁻² for the criterion (a); 194.02 M_⊙ Myr⁻¹ pc⁻² for the criterion (b); 715.87 M_⊙ Myr⁻¹ pc⁻² for the criterion (c); and 438.38 M_⊙ Myr⁻¹ pc⁻² for ℳ = 1 for the criterion (d).

Observational studies, such as in Palmeirim et al. (2013), provided an estimate of the accretion rate in the B211 filament based on the observed mass per unit length and the ¹²CO (1− 0) inflow velocity, finding a value of ~27–50 M_⊙ Myr⁻¹ pc⁻¹. With this accretion rate, it would take ~1–2 Myr to form the B211 filament, in reasonable agreement with the free-fall time of ~1.3 Myr of our model. Since the accretion rate is estimated in the observations at a distance of 0.4 pc from the filament axis, the method of the width is not appropriate for a comparison, because of the smaller radius used and because the width of the filament changes with time. The methods based on Mach number and 0.8M_tot are not well suited for a comparison, as they assume a position too close to the filament center. On the other hand, the method of the velocity-divergence (which is the chosen method) displays a greater similarity to the observations, where the inflow velocity is ~0.6–1.1 km s⁻¹; whereas in our model, the velocity in the x-axis at the free-fall time t = t_ff is ~1.2 km s⁻¹ at ~0.1 pc from the center of the filament.

Fig. A.1 Time evolution of the mass flux according to the different methods adopted here. Circles correspond to non-ideal MHD, stars to ideal MHD.

References

All Tables

All Figures

Fig. 1 Representative sketch of our set-up. The code evolves a vector variable with three components (x, y, z,) along one coordinate (x); hence, y and z have periodic boundary conditions. The x coordinate traverses the filament (shaded gradient) and the initial density distribution along this component follows a Plummer-like profile with maximum density, ρridge, and background density, ρ0, see Eq. (1). The gas accretes along the same axis, as indicated by the arrows. The magnetic field is inclined by an angle θ with respect to the x-axis, where θ = arctan By/Bx. We note that the color palette employed for the filament is only included for illustration purposes.

In the text

Fig. 2 Profile along the x-axis of various physical quantities of the model at three different times: t/tff = 0.1, 0.5 and 1.0. (a) density ρ, (b) temperature T, (c) x-component of the velocity υx, (d) total velocity magnitude |v|, (e) y-component of magnetic field By, (f) Mach number ℳ, (g) total magnetic field |B| magnitude, (h) mass flux $\dot{\text{Σ}}$, (i) ambipolar diffusion resistivity coefficient ηAD, and (j) abundance e−/H2. The solid lines correspond to non-ideal MHD, dashed lines to ideal MHD. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments. We note that in panel (j), in the ideal case (dashed lines), ηAD is reported for the sake of comparison, but it is not included in the MHD equations.

In the text

	Fig. 3 Density and velocity along the x-axis for different initial angles. The solid line corresponds to the reference case. At 0.5tff in the top panels and at 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.
In the text

	Fig. 4 Inclination of the magnetic field and y-component for the different initial angles. The solid line corresponds to the reference case. We have t = 0.5tff in the top panels and t = 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.
In the text

	Fig. 5 Density and velocity along the x-axis for the different initial magnitudes of the magnetic field. The solid line correspond to the reference case. We have t = 0.5tff in the top panels and t = 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.
In the text

	Fig. 6 Density and velocity on the x-axis for different cosmic ray ionization rates. The black line correspond to the reference case. We have: t = 0.5tff in the top panels and t = 1.0tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.
In the text

	Fig. 7 Ionization fraction and ambipolar diffusion resistivity coefficient for the different cosmic ray ionization rates. The black line correspond to the reference case. We have t = 0.5tff in the top panels and t=1.0 tff in the bottom panels. The gray area in the middle corresponds to 0.1 pc to guide the reader to identify the typical width of the filaments.
In the text

	Fig. 8 Time evolution of the FWHM ratio of ℳ = 4 (dot-dashed line) and ℳ = 6 (dashed line) with ℳ = 2 (reference).
In the text

	Fig. A.1 Time evolution of the mass flux according to the different methods adopted here. Circles correspond to non-ideal MHD, stars to ideal MHD.
In the text