Growth of Massive Molecular Cloud Filament by Accretion Flows. I. Slow-shock Instability versus Ambipolar Diffusion

The Herschel Gould Belt Survey showed that stars form in dense filaments in nearby molecular clouds. Recent studies suggest that massive filaments are bound by the slow shocks caused by accretion flows onto the filaments. The slow shocks are known to be unstable to corrugation deformation of the shock front. Corrugation instability could convert the accretion flow's ram pressure into turbulent pressure that influences the width of the filament, which, according to theory, determines the self-gravitational fragmentation scale and core mass. In spite of its importance, the effect of slow-shock instability on star-forming filaments has not been investigated. In addition, the linear dispersion relation obtained from ideal magnetohydrodynamics (MHD) analysis shows that the most unstable wavelength of shock corrugation is infinitesimally small. In the scale of dense filaments, the effect of ambipolar diffusion can suppress the instability at small scales. This study investigates the influence of ambipolar diffusion on the instability of the slow shock. We perform two-dimensional MHD simulations to examine the linear growth of the slow-shock instability, considering the effect of ambipolar diffusion. The results demonstrate that the most unstable scale of slow-shock instability is approximately 5 times the length scale of ambipolar diffusion ℓ AD calculated using post-shock variables, where ℓ AD corresponds to the scale where the magnetic Reynolds number for ambipolar diffusivity is unity.

In general, mechanisms G, C, and O result in supercritical filaments, which are of interest for star formation.The common feature of these mechanisms is that filaments are formed by gas flow along the local magnetic field in a shock-compressed sheet.
Several recent works highlight the importance of accretion in the context of a filament's evolution.For example, molecular emission-line observations provide evidence of perpendicular accretion onto filaments (Palmeirim et al. 2013;Shimajiri et al. 2019;Chen et al. 2020).In particular, Shimajiri et al. (2019) reported the occurrence of accretion onto filaments in a shocked sheet.Clarke et al. (2016) demonstrated that the most unstable length scale for self-gravitational fragmentation along a filament depends on the accretion rate onto the filament.Hennebelle & André (2013) developed an analytical model that can be applied to self-gravitating and accreting filaments.They considered turbulence driven by accretion onto the filament and its dissipation by ion-neutral friction.
A shock wave with an Alfvén mach number 1 A  < and a sonic Mach number 1 s  > is called a "slow(-mode) shock."Given that massive filaments are formed in the post-shock layer threaded by a strong magnetic field with energy exceeding the kinetic energy of accretion flows (the Type O mechanism; Inoue et al. 2018), the filament surface is naturally bound by slow shocks.Lessen & Deshpande (1967) found through linear stability analysis that the slow-shock front is corrugationally unstable (slow-shock instability, hereafter SSI).The mechanism of SSI is explained as follows.In contrast to the case of fast shock, the component of the magnetic field tangential to the shock surface decreases across the front.Thus, when the shock corrugates, the magnetic field lines kink as denoted by the red lines in Figure 1.Because gas flows along the magnetic field, the gas converges behind the peak of the shock, while it diverges behind the valley.Such flow patterns increase (decrease) the pressure behind the peaks (valleys), which further pushes up (pulls down) the shock front.Édel'Man (1989) showed that the approximated dispersion relation of SSI for  1 A  can be written as where ω, v sh , and k denote the frequency, shock velocity, and wavenumber of the shock corrugation.As a more accurate solution, Édel'Man (1989) derived the approximate dispersion relation for γ = 5/3: Equations (1) and (2) indicate that the most unstable scale is infinitesimally small.This unphysical feature stems from the ideal approximation and the resulting discontinuous treatment of the shock.To know the physical scale length of the SSI, we consider a nonideal effect.Since the corrugation of the shock front generally produces turbulent flows behind the shock (e.g., Inoue & Inutsuka 2012;Inoue et al. 2012), we can expect that the SSI will deposit additional energy to the filament.
In molecular clouds, the ambipolar diffusion is effective and potentially modifies the SSI dynamics.The magnetic Reynolds number of the flow with ambipolar diffusivity is given as where γ in ≡ 〈σ in v in 〉/(m + m i ) = 3.5 × 10 13 cm 3 g −1 s −1 and ρ i denotes the ion mass density; σ in , v in , m, and m i represent the ion-neutral cross section (the Langevin cross section), the relative velocity between a neutral molecule and an ion, the mean molecule mass, and the mean ion mass, respectively.Assuming a balance between ionization by cosmic rays and recombination, ρ i can be expressed as Cρ 1/2 .In this study we apply C = 3 × 10 −16 cm −3/2 g 1/2 (Shu 1992).The characteristic length scale below which the effect of ambipolar diffusion becomes nonnegligible can be obtained by solving Comparing the actual observation, for densities around 5000 cm −3 , the mean magnetic field from the Crutcher (2012) plot (Figure 6 in that review) is about 5 μG.Then Equation (4) gives ℓ AD ∼ 2.2 × 10 −4 pc.The observed maximum magnetic field strength for the same density range is about 50 μG, for which we get ℓ AD ∼ 2.2 × 10 −2 pc.Therefore, ℓ AD can take a wide range of values.This suggests that the characteristic scale of the ambipolar diffusion can be comparable to the filament width and hence the ambipolar diffusion can affect the filament dynamics.In the context of the solar chromosphere, Snow & Hillier (2021) performed twodimensional two-fluid simulations of SSI for partially ionized gas.They demonstrated that the neutral fluid stabilizes the SSI on a small scale and found new features such as gas accumulation at valleys.However, the situation in their simulations is different from the one in molecular clouds (e.g., the ionization degree, the ion-neutral collision cross section, etc.), and they did not study the dispersion relation and the dependence on density and the magnetic field.The linear analysis of SSI including ambipolar diffusion is fundamentally challenging although there is a recent study that has challenged it (Béthune 2023).Our strategy is to directly simulate the SSI including the ambipolar diffusion and the measurement of the growth rate.
In this paper, as a first step to understanding the effect of SSI on filaments, we study the effect of ambipolar diffusion on the SSI and derive the most unstable scale.As a result of this study, we can determine the typical length scale of SSI in filaments that will meet the resolution requirement in future simulations.More realistic simulations of filament evolution with slow shocks caused by the converging accretion including self- Type C perpendicular Gas coagulation along the magnetic field by local turbulent velocity perturbations within shock-compressed layers.
Type S parallel Shear flows associated with turbulence stretch in existing clumps.gravity will be the object of our future studies.The paper is organized as follows: In Section 2, we provide the setup of our simulations, and we show and interpret the results in Section 3.
In Section 4, we discuss the stabilizing scale of SSI versus ambipolar diffusion.Finally, we summarize the results in Section 5.

Setup for Simulations
We perform two-dimensional and three-dimensional ideal/ nonideal magnetohydrodynamic (MHD) simulations using the Athena++ code (Stone et al. 2020).To determine the physical scale of the SSI, two-dimensional simulation seems to be sufficient, because linear stability analysis does not show any difference.In addition, by ideal MHD simulations, Stone & Edelman (1995) demonstrated that the growth rates of the SSI in three-dimensional cases are not different from those in twodimensional cases.We confirm this expectation even with the effect of ambipolar diffusion in Section 3.2.1.We use the second-order-accurate van Leer predictor-corrector scheme and piecewise linear method applied to primitive variables to integrate the equations.The constrained transport method (Stone & Gardiner 2009) ensures the divergence-free condition, ∇ • B = 0.In this paper, we do not solve the Poisson equation for self-gravity because we concentrate on the physics of SSI under the influence of ambipolar diffusion as the first step of this sort of study.The effect of self-gravity will be considered in our future studies.We solve the following equations: where P * = p + B 2 /(8π) and E = e + ρv 2 /2 + B 2 /(8π) denote the total pressure and total energy density; ρ, p, v, and B represent the density, pressure, velocity, and magnetic field; and J = ∇ × B represents the current.We introduce the viscous stress tensor to prevent the carbuncle phenomenon (Quirk 1994;Liou 2000; Kim et al. 2003) and the growth of a grid-scale SSI seeded by the carbuncle instability; ν denotes the coefficient of physical kinematic viscosity, which is adjusted to stabilize a grid-scale (eight cells) fluctuation.The box size L box and ν are chosen so that the stabilizing scale by ambipolar diffusion is sufficiently larger than this grid scale.The variable η AD denotes the ambipolar diffusion coefficient, which is given by where ρ n ; ρ represents the neutral gas mass density, and the ion mass density is denoted by ρ i = Cρ 1/2 .In this study we apply C = 3 × 10 −16 cm −3/2 g 1/2 (Shu 1992).
We numerically solve Equations ( 5)-( 8) on a two-dimensional domain of size [−4L box , 4L box ] × [0 pc, L box ] in the shock rest frame.We select L box = 0.2, 0.25, or 0.5 pc.The specific heat ratio γ = 1.01 is used.The initial density, velocity, and pressure fields are set as respectively, where ρ 0 , v x0 , and p 0 denote the initial density, xcomponent of velocity, and pressure in the pre-shock region, respectively.We set the upstream gas sound speed c s to 0.2 km s −1 so that p 0 is given by p c . The compression ratio r and pressure jump r pres can be written as As a seed of instability, the density perturbation is introduced as follows: For a three-dimensional simulation, the density perturbation is where λ p denotes the wavelength of perturbation.These initial conditions lead to a perturbed stationary shock at x = 0. Since the star-forming filaments are perpendicular to the magnetic field, the initial uniform magnetic field is set along the x-axis B x 0 ˆ.The numerical domain is a two-dimensional box with a uniform grid of 4096 × 512 cells, which leads to a spatial resolution of Δx = L box /512.We apply zero-gradient boundary conditions (with continuous gas flow) at the boundaries x = −4L box and x = 4L box .For y = 0, L box boundaries, we use the periodic boundary conditions.

SSI in an Ideal MHD Case
Since isothermal treatment is justified in dense regions of molecular clouds, we adopt γ = 1.01.We use the Roe solver because of its more numerically stable features in nonlinear regimes (see Appendix B).The method to measure the growth rate of the SSI is the same as the method developed by Stone & Edelman (1995).It is convenient to use compression-weighted averages because we need to evaluate physical quantities in the vicinity of a shock wave.The compression-weighted transverse magnetic energy is written as follows: In Figure 2, we show the evolution of B y 2 á ñ in model n1000b30v1.We can confirm linear growth for λ p = 0.005-0.2pc modes for t ∼ 0.1-0.6Myr.For λ p = 0.002 pc, B B y 2 0 2 á ñ decreases until ∼0.2 Myr, then a larger-scale (>0.002 pc) grid noise grows after t ∼ 0.2 Myr, which is different from the growth of the λ p = 0.002 pc mode of the SSI.Also for λ p = 0.2, 0.1, 0.02, and 0.005 pc, we can see the slope increments of the perturbed magnetic field after t = 0.7 Myr caused by the grid noise.(In these cases, the scale of noise is smaller than λ p .) B B Since the initial perturbation is not given as the eigenstate of the SSI, the SSI does not start growing at t = 0. Thus we define the measuring range as [t 0 = t start + ft growth , t 0 + t range ] to observe the linear growth of the SSI, where t growth ≡ 1/ω ana is the growth timescale.In this section, we select t start = 0.05 Myr, t range = 0.5 Myr, and f = 0.6.We show the dispersion relation for the isothermal ideal MHD case including the physical shear viscosity (model n1000b30v1, 19.5 shear  = ) as gray cross marks in Figure 3.The vertical dotted line represents the scale of 8Δx.We find that if we use 19.5 shear  = , a physical dispersion relation is successfully obtained by suppressing the carbuncle phenomenon.

SSI versus Ambipolar Diffusion
We perform a similar analysis to that in Section 3.1 for the simulation results including ambipolar diffusion.In Figure 4, we show the evolution of the mean value of the perturbed magnetic field for model n1000b30v1AD.Because the ambipolar diffusion diminishes the phase speed of the Alfvén wave, the eigenstate requires more time than the ideal MHD case to develop from the given initial perturbation.Thus, to measure the linear phase growth rate, we take a longer t start of 1.2 Myr, and f = 0.6 and t range = 0.3 Myr.We show the dispersion relation for the isothermal MHD case including ambipolar diffusion (model n1000b30v1AD) as black cross marks in Figure 3.We can see the reduction of the SSI growth by the ambipolar diffusion.We find that the most unstable scale  ℓ 0.2 max pc and the damping scale ℓ damp ; 0.02 pc.To investigate the parameter dependence on the most unstable scale, we perform a parameter survey for the unperturbed magnetic field strength, density, and shock velocity.In Figure 5, we show the dispersion relation for models n1000b24v1AD, n1000b40v1AD, n800b30v1AD, n1600b24v1AD, and n1000b30v0.8AD.We can confirm that the most unstable scale varies with the magnetic field density and velocity, such that a higher density/velocity corresponds to a smaller most unstable scale.Conversely, a larger magnetic field corresponds to a larger most unstable scale.These trends can be understood based on the scale of ambipolar diffusion 4)).

Three-dimensional Simulations
We also perform three-dimensional simulations to investigate more realistic cases in molecular clouds.Because of the   higher computational cost for three-dimensional simulations, we use the super-time-stepping method for the diffusion term.In Appendix C, we show the results of tests for the super-timestepping method (Meyer et al. 2014) and confirm that the results do not change even with considerable acceleration of the calculations.In Figure 6, we show the dispersion relations for models n1000b24v1AD3D (black solid line) and n1000b24v1AD (gray dashed line).We measure the growth rate in the same way as in Section 3.2.We can confirm that those two are not different despite the different dimensions.

Discussion
In Section 3.2, we have stated that the most unstable scale depends on the magnetic field, density, and shock velocity.In this section, we discuss how the most unstable wavelength of the SSI is scaled.In Figure 7, we show the damping length scale, which scales with the most unstable length, as a function of the initial density (panel (a)), the initial magnetic field (panel (b)), and the shock velocity (panel (c)).The gray dashed line is the scale of ambipolar diffusion evaluated in post-shock quantities, which can be written as  - denote the postshock magnetic field strength, density, and velocity, respectively.The gray solid line represents 1.7 × ℓ AD,post .The top and bottom edges of the vertical lines denote the minimum lengths Figure 5. Top panel: effect of varying the initial density.Blue, black, and green crosses represent the growth rate for models n800b30v1AD, n1000b30v1AD, and n1600b30v1AD, respectively.The vertical dashed line for each color represents the 1/(1.7ℓAD ) for the model corresponding to that color.Middle panel: same as top panel, but for cases with varying initial magnetic field (models n1000b24v1AD (blue), n1000b30v1AD (black), and n1000b40v1AD (green)).Bottom panel: same as top panel, but for cases with varying shock velocity (models n1000b30v1AD (black), n1000b30v0.9AD(green), and n1000b30v0.8AD(blue)).with a positive growth rate and the maximum lengths with a negative growth rate, i.e., the vertical lines show the range in which the damping scale exists.For panels (a) and (b), the damping length scale well follows 1.7ℓ AD,post .The dependence on shock velocity deviates more strongly from the predictions of Equation ( 21), but the difference from 1.7ℓ AD,post is within a factor of 2. While our estimation for the length scale of ambipolar diffusion utilizes physical quantities from the downstream region, a more appropriate approach should be based on the shock transition layer, where ambipolar diffusion actually works.The dependency on the magnetic field and upstream density in Equation (21) remains consistent because the change of the magnetic field or upstream density does not shift the compression ratio (Figures 7(a) and (b)).However, in contrast to the density and magnetic field strength, the shock velocity affects the compression ratio in the isothermal case ), which brings some error to the ambipolar diffusion scale estimation (Equation ( 21)).The density in the transition layer is lower than the downstream density, which implies that the actual length scale of ambipolar diffusion is expected to be larger than Equation (21).This may account for the observed deviation between the simulation results and Equation (21).We conclude that the characteristic damping scale of the SSI can be approximately estimated as ℓ damp ∼ 1.7ℓ AD,post , and the most unstable scale as ℓ SSI ∼ 5ℓ damp ∼ 9ℓ AD,post .
The filament width is one of the significant quantities that determine the initial conditions for star formation.The critical line mass including the support of the magnetic field depends on the filament width (Tomisaka 2014).According to linear theory, the self-gravitational fragmentation length scale of the filament depends on the filament width (Stodólkiewicz 1963;Inutsuka & Miyama 1992).Arzoumanian et al. (2011Arzoumanian et al. ( , 2019) ) found that the characteristic width of the Herschel Gould Belt filaments is 0.1 pc (see also Juvela et al. 2012;Koch & Rosolowsky 2015).Remarkably, filaments maintain their width even when their line mass exceeds 100 M e pc −1 .It should be noted that some studies have questioned the universality of 0.1 pc filament widths.For example, Ossenkopf-Okada & Stepanov (2019) used wavelet decomposition to analyze Herschel survey data and did not find a characteristic length scale.Panopoulou et al. (2022aPanopoulou et al. ( , 2022b) ) suggested that the estimated filament width appears to depend weakly on the distance to the observed cloud.On the other hand, André et al. (2022) show that the finding by Panopoulou (2022aPanopoulou ( , 2022b) ) of "an apparent width depending on distance" is entirely consistent with an intrinsic filament width of ∼ 0.1 pc when finite resolution/beam convolution effects are properly taken into account.Considering only thermal support against gravity, such a high-line-mass structure cannot maintain a length scale of 0.1 pc.Several authors have studied the effects of turbulence and/or the magnetic field and have shown that subcritical and mildly supercritical filaments have a width of 0.1 pc (Fischera & Martin 2012;Auddy et al. 2016;Federrath 2016;Priestley & Whitworth 2022); however, the reason for the constant width of filaments, especially of massive filaments, remains a mystery.Seifried & Walch simulated the evolution of a massive filament with a line mass of 75 M e pc −1 and found that filaments perpendicular to the magnetic field are thinner.They claimed that their very narrow filaments could be interpreted as the fibers reported by Hacar et al. (2013).Smith et al. (2014) performed simulations of filament formation, which indicated that the width of a massive filament is approximately 0.3 pc, which is inconsistent with the constant 0.1 pc filament width proposed by Arzoumanian et al. (2011Arzoumanian et al. ( , 2019)).We need to understand the origin of the universal width, especially for massive filaments, and can do so by examining the detailed process of gas accretion flows onto filaments.

Summary
We performed SSI simulations of molecular clouds using two-dimensional isothermal (γ = 1.01) or adiabatic (γ = 5/3) models using Athena++.Furthermore, we investigated the most unstable length scale of SSI in molecular clouds caused by ambipolar diffusion, which provides the resolution requirement in future simulations.The major findings of this study are described as follows.
1. Ambipolar diffusion suppresses SSI on a small scale.We find that the most unstable scale is of the order of ∼0.1 pc and the damping scale is of the order of ∼0.01 pc in molecular clouds.2. The most unstable and damping length scales depend on the density, shock velocity, and magnetic field strength.The scaling is roughly described by ℓ ℓ 5 In a follow-up paper (Paper II), we will show that a natural extension of this work leads to a better understanding of the filament evolutionary process.In realistic situations, filaments are bound by two shocks.Because the separation of these two shocks is narrow and the two shock surfaces share the magnetic field lines that thread them, we can expect that the two shocks dynamically influence each other.In Paper II, we will discuss the effect of the two interacting shocks and demonstrate that the SSI creates inhomogeneous post-shock flows and can provide additional dynamical pressure to the filament.).The simulation using HLLD suffers from the carbuncle phenomenon.The LHLLD scheme is designed to alleviate the carbuncle phenomenon, but after a long-time integration, growth of grid noise occurs.Both schemes provide similar results, but we fail to measure for λ p > 1.0 pc modes if we do not involve physical shear viscosity by the effects of the carbuncle phenomenon due to the slower growth of the SSI than of the carbuncle phenomenon (see filled blue circles and blue plus marks).It should be mentioned that the results using LHLLD are closer to the approximated analytical solution than those using HLLD at λ p 0.5 pc.The simulations using either HLLD or LHLLD successfully reproduce the growth of SSI for 9.8 shear  = (see pink circles and plus marks).In conclusion, we find that the HLLD, LHLLD, and Roe solvers with adjusted physical shear viscosity can correctly calculate the growth rate of SSI over a wide scale range.We can use any of the HLLD, LHLLD, and Roe solvers to measure the linear growth rate, but in the following sections, we use the Roe solver because of its more numerically stable features in nonlinear regimes (see Appendix B).

Appendix B Unphysical Numerical Explosion in HLLD/LHLLD
In Figures 10 and 11, we show snapshots of the density and pressure maps of models n1000b30v1.3aD ¥ and n1000b30v1.3aLD ¥ , respectively.For long-term simulations using HLLD/LHLLD without the physical shear viscosity, the numerical errors around the shock front cause unphysical numerical explosions, which do not appear in simulations with the Roe scheme.Miyoshi & Kusano (2005) is close to zero.The latest version of Athena++ has been designed to prevent this issue to some extent, but it cannot prevent the unphysical explosion under the initial conditions dealt with in this study.Although such a numerical effect can be quenched by physical viscosity, we select the Roe solver with a physical shear viscosity to ensure safe long-term integration.(red).We can confirm that the results do not change significantly even when dt dt max 1000 parabolic ( )= .
to −4 after t = 0.75 Myr for λ p = 0.01 due to the saturation of SSI (Stone & Edelman 1995).The slope of each line in Figure2reflects the growth rate.The growth rate can be measured from the slope of B y 2

Figure 2 .
Figure 2. Evolution of mean value of the perturbed magnetic field for the case with isothermal ideal MHD including physical shear viscosity (model n1000b30v1).

Figure 4 .
Figure 4. Evolution of mean value of the perturbed magnetic field for cases with isothermal ideal MHD including ambipolar diffusion and physical shear viscosity (model n1000b30v1AD).

Figure 6 .
Figure 6.The dependence on the number of spatial dimensions (models n1000b3024v1AD3D (black) and n1000b24v1AD (gray)).The k is defined as k k 2 2 y z 2 2 p p l + = for model n1000b3024v1AD3D.