Kelvin-Helmholtz Instability in Self-Gravitating Streams

TL;DR

This paper investigates how Kelvin-Helmholtz instability affects the evolution and fragmentation of self-gravitating gaseous streams in astrophysical contexts, revealing conditions for stability, disruption, and clump formation.

Contribution

It introduces an analytic model and hydrodynamic simulations to determine the critical line-mass for stability and describes the nonlinear evolution of KHI in self-gravitating streams.

Findings

Critical line-mass depends on Mach number and density contrast.

Gravity delays stream disruption by halting shear layer penetration.

Streams above critical line-mass fragment into long-lived clumps.

Abstract

Self-gravitating gaseous filaments exist on many astrophysical scales, from sub-pc filaments in the interstellar medium to Mpc scale streams feeding galaxies from the cosmic web. These filaments are often subject to Kelvin-Helmotz Instability (KHI) due to shearing against a confining background medium. We study the nonlinear evolution of KHI in pressure-confined self-gravitating gas streams initially in hydrostatic equilibrium, using analytic models and hydrodynamic simulations, not including radiative cooling. We derive a critical line-mass as a function of the stream Mach number and density contrast with respect to the background, $\mu_{cr}(M_b,\delta_c)\le 1$, where $\mu=1$ is normalized to the maximal line mass for which initial hydrostatic equilibrium is possible. For $\mu<\mu_{cr}$, KHI dominates the stream evolution. A turbulent shear layer expands into the background and leads to stream deceleration at a similar rate to the non-gravitating case. However, with gravity, penetration of the shear layer into the stream is halted at roughly half the initial stream radius by stabilizing buoyancy forces, significantly delaying total stream disruption. Streams with $\mu_{cr}<\mu\le 1$ fragment and form round, long-lived clumps by gravitational instability (GI), with typical separations roughly 8 times the stream radius, similar to the case without KHI. When KHI is still somewhat effective, these clumps are below the spherical Jeans mass and are partially confined by external pressure, but they approach the Jeans mass as $\mu\rightarrow 1$ and GI dominates. We discuss potential applications of our results to filaments in the ISM and dense streams feeding galaxies at high redshift.

Tables

Table 1. Table 1: Parameters of simulations used for studying the evolution of streams undergoing both GI and KHI. The first three columns list the control parameters, namely the line-mass ratio μ 𝜇 \mu , Mach number M b subscript 𝑀 b M_{\rm b} , and the ratio of central density to background density at the stream boundary δ c subscript 𝛿 c \delta_{\rm c} . The remaining six columns list derived parameters: the ratio of stream to background density on either side of the interface, δ 𝛿 \delta , the GI time scale, t max subscript 𝑡 max t_{\rm max} , the timescale for KHI to destroy the contact discontinuity, t shear subscript 𝑡 shear t_{\rm shear} , the timescale for KHI to destroy the entire stream, t dis subscript 𝑡 dis t_{\rm dis} , the stream sound crossing time, t sc subscript 𝑡 sc t_{\rm sc} , and the shortest unstable wavelength for GI, λ cr subscript 𝜆 cr \lambda_{\rm cr} . All timescales are in units of the stream free-fall time, t ff subscript 𝑡 ff t_{\rm ff} , while λ cr subscript 𝜆 cr \lambda_{\rm cr} is in units of the stream radius, R s subscript 𝑅 s R_{\rm s} . For all cases, the fastest growing wavelength for GI is λ max ∼ < 2 ​ λ cr superscript similar-to subscript 𝜆 max 2 subscript 𝜆 cr \lambda_{\rm max}\lower 2.58333pt\hbox{$\;\buildrel<\over{\sim}\;$}2\lambda_{\rm cr} .

$\mu$	$M_b$	${\delta }_c$	$\delta$	$t_{\max }$	$t_{\mathrm{shear}}$	$t_{\mathrm{dis}}$	$t_{\mathrm{sc}}$	${\lambda }_{cr}$
0.1	1	100	92	2.39	1.34	14.15	0.64	5.95
0.2	1	100	84	2.46	2.01	20.43	0.90	5.95
0.3	1	100	76	2.54	2.65	25.65	1.11	5.95
0.4	1	100	67	2.63	3.32	30.51	1.28	5.95
0.5	1	100	58	2.76	4.09	35.38	1.44	5.98
0.6	1	100	49	2.93	5.06	40.62	1.58	5.98
0.7	1	100	40	3.18	6.39	46.77	1.73	6.02
0.8	1	100	30	3.60	8.52	55.01	1.88	6.09
0.9	1	100	18	4.56	13.20	69.68	2.06	6.23
0.1	1	6.7	6.2	3.21	2.00	6.97	0.64	7.02
0.2	1	6.7	5.6	3.41	3.04	10.26	0.90	7.16
0.3	1	6.7	5.1	3.66	4.07	13.22	1.11	7.36
0.5	2.5	100	58	2.76	2.27	19.57	1.44	5.98
0.6	2.5	100	49	2.93	2.84	22.77	1.58	5.98
0.7	2.5	100	40	3.18	3.65	26.71	1.73	6.02
0.9	2.5	100	18	4.56	8.22	43.40	2.06	6.23
0.7	6	100	40	3.18	1.52	11.13	1.73	6.02
0.8	6	100	30	3.60	2.09	13.47	1.88	6.09
0.9	6	100	18	4.56	3.43	18.09	2.06	6.23

Equations

$$\rho (r)={\rho }_{\mathrm{c}}{\left[1+\frac{1}{8}{\left(\frac{r}{H}\right)}^2\right]}^{-2},H=\frac{c_s}{\sqrt{4\pi G{\rho }_{\mathrm{c}}}},$$

$$\Lambda ={\int }_0^{R_{\mathrm{s}}}2\pi r\rho (r)\mathrm{dr}.$$

$${\Lambda }_{\mathrm{cr},\mathrm{iso}}=2c_{\mathrm{s}}^2/G.$$

$$\frac{\Lambda }{{\Lambda }_{\mathrm{cr},\mathrm{iso}}}={\left[1+8{\left(\frac{H}{R_{\mathrm{s}}}\right)}^2\right]}^{-1}.$$

$$P_{\mathrm{ext}}=P(R_{\mathrm{s}})=c_{\mathrm{s}}^2\rho (R_{\mathrm{s}})=c_{\mathrm{s}}^2{\rho }_{\mathrm{c}}{\left[1+\frac{R_{\mathrm{s}}^2}{8H^2}\right]}^{-2}.$$

$$\frac{\Lambda }{{\Lambda }_{\mathrm{cr},\mathrm{iso}}}=1-\frac{P_{\mathrm{ext}}}{{\rho }_{\mathrm{c}}{c}_{\mathrm{s}}^2}.$$

$$M_{\mathrm{BE}}=1.18\frac{c_{\mathrm{s}}^4}{P_{\mathrm{ext}}^{1/2}{G}^{3/2}}$$

$$t_{\mathrm{ff}}=(4G<\rho >{)}^{-1/2}.$$

$$<\rho >={\rho }_{\mathrm{c}}\left(1-\frac{\Lambda }{{\Lambda }_{\mathrm{cr},\mathrm{iso}}}\right)={\rho }_{\mathrm{c}}{\left(1+\frac{R_{\mathrm{s}}^2}{8H^2}\right)}^{-1}.$$

$$\frac{{\omega }^2}{4\pi G\rho }=-\frac{x{I}_1}{I_0}\left[K_0{I}_0-\frac{1}{2}\right],$$

$$\begin{array}[]{c}\dfrac{\omega^{2}}{4\pi G\bar{\rho}}=-\left[\dfrac{x(\rho_{\rm s}-\rho_{\rm b})^{2}I_{0}K_{0}}{\bar{\rho}^{2}}-\dfrac{x\rho_{\rm s}(\rho_{\rm s}-\rho_{\rm b})}{2\bar{\rho}^{2}}\right]\\ \times\left[\dfrac{\rho_{\rm s}I_{0}}{\bar{\rho}I_{1}}+\dfrac{\rho_{\rm b}K_{0}}{\bar{\rho}K_{1}}\right]^{-1},\end{array}$$

$$I_0{K}_0>\frac{1}{2\left(1-{\delta }^{-1}\right)},$$

$$M_{\mathrm{tot}}=\frac{V_{\mathrm{s}}}{c_{\mathrm{s}}+c_{\mathrm{b}}}.$$

$$h=\alpha V_{\mathrm{s}}t$$

$$h_{\mathrm{s}}=\frac{\alpha V_{\mathrm{s}}t}{1+\sqrt{\delta }},h_{\mathrm{b}}=\frac{\sqrt{\delta }\alpha V_{\mathrm{s}}t}{1+\sqrt{\delta }}.$$

$$t_{\mathrm{dis}}=\frac{\left(1+\sqrt{\delta }\right)R_{\mathrm{s}}}{\alpha V_{\mathrm{s}}}.$$

$$t_{\mathrm{shear}}=\frac{R_{\mathrm{s}}}{\alpha V_{\mathrm{s}}}.$$

$$V_{\mathrm{s}}(t)=\frac{V_{\mathrm{s},0}}{1+t/t_{\mathrm{dec}}},$$

$$t_{\mathrm{dec}}=\frac{\left(1+\sqrt{\delta }\right)\left(\sqrt{1+\delta }-1\right)}{\alpha \sqrt{\delta }}\frac{R_{\mathrm{s}}}{V_{\mathrm{s},0}},$$

$$\alpha \simeq 0.21\times \left[0.8\exp \left(-3M_{\mathrm{tot}}^2\right)+0.2\right].$$

$${\omega }^2=-\frac{{\rho }_{\mathrm{s}}{\rho }_{\mathrm{b}}}{({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{b}}{)}^2}{V}^2{k}^2+\frac{{\rho }_{\mathrm{s}}-{\rho }_{\mathrm{b}}}{{\rho }_{\mathrm{s}}+{\rho }_{\mathrm{b}}}kg,$$

$${\lambda }_{\max }=\frac{2\pi {\rho }_{\mathrm{s}}{\rho }_{\mathrm{b}}{V}^2}{\left({\rho }_{\mathrm{s}}^2-{\rho }_{\mathrm{b}}^2\right)g}.$$

$$\vec{\nabla }P=-\rho \vec{\nabla }\Phi ,$$

$${\nabla }^2\Phi =4\pi G\rho ,$$

$$P=K{\rho }^{\gamma },$$

$$\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{r}{\rho }\frac{\partial \left(K{\rho }^{\gamma }\right)}{\partial r}\right]=-4\pi G\rho ,$$

$$\rho(r=0)=\rho_{\rm c},\quad\frac{\partial\rho}{\partial r}\bigg{|}_{r=0}=0.$$

$$H^2=\frac{c_{\mathrm{s},0}^2}{(\gamma -1)4\pi G{\rho }_{\mathrm{c}}},$$

$$\frac{1}{x}\frac{\partial}{\partial x}\left(x\frac{\partial y^{\gamma-1}}{\partial x}\right)=-y,\quad y(0)=1,\quad\frac{\partial y}{\partial x}\bigg{|}_{0}=0.$$

$${\Lambda }_{\mathrm{cr}}=\frac{c_{\mathrm{s},0}^2}{2(\gamma -1)G}{\int }_0^{R_{\mathrm{equ}}/H}y(x)\times x\mathrm{d}x=a\frac{c_{\mathrm{s},0}^2}{G},$$

Full text

††pagerange: Kelvin-Helmholtz Instability in Self-Gravitating Streams–LABEL:lastpage††pubyear: 0000

1 Introduction

Filamentary structures are present across many astrophysical scales, from $\,{\rm Mpc}$ to sub-$\,{\rm pc}$. On the largest scales, structure formation occurs in the “cosmic web”, a network of sheets and filaments that connect dark matter haloes (Zel’dovich, 1970; Bond et al., 1996; Springel et al., 2005), and is also evident in the distributions of galaxies (e.g. Colless et al., 2003; Tegmark et al., 2004; Huchra et al., 2005). Intergalactic gas cools and condenses towards the centres of the dark matter filaments, forming a network of baryon-dominated intergalactic gas streams (Dekel & Birnboim, 2006; Birnboim et al., 2016). There have been several recent attempts to model such streams self-gravitating Mpc-scale gaseous cylinders, which seems consistent with cosmological simulations (Harford et al., 2008; Harford & Hamilton, 2011; Freundlich et al., 2014; Mandelker et al., 2018).

At the nodes of the cosmic web, the most massive haloes reside at the intersection of several filaments and are penetrated by the gas streams residing at their centres. These streams constitute the main mode of gas accretion onto the central galaxies (Kereš et al., 2005; Dekel et al., 2009a; Danovich et al., 2012; Zinger et al., 2016). At redshifts z\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}2, simulations suggest that streams feeding galactic haloes remain dense and cold, with temperatures of $\sim 10^{4}\,{\rm K}$, as they travel through the hot circumgalactic medium (CGM) towards the central galaxy (Kereš et al., 2005; Dekel & Birnboim, 2006; Ocvirk et al., 2008; Dekel et al., 2009a; Ceverino et al., 2010; Faucher-Giguère et al., 2011; van de Voort et al., 2011, though see also Nelson et al., 2013, 2016). The filamentary structure in such systems can thus be maintained down to scales of tens of $\,{\rm kpc}$ around galaxies (though see below), where it has been suggested that they may fragment due to gravitational instability (hereafter GI; Dekel et al., 2009b; Genel et al., 2012; Mandelker et al., 2018). While these cold circumgalactic streams are difficult to directly detect, recent observations have revealed massive extended cold components in the CGM of high-redshift galaxies, whose spatial and kinematic properties are consistent with predictions for cold streams (Bouché et al., 2013, 2016; Prochaska et al., 2014; Cantalupo et al., 2014; Martin et al., 2014a; Martin et al., 2014b; Borisova et al., 2016; Fumagalli et al., 2017; Leclercq et al., 2017; Arrigoni Battaia et al., 2018).

Within galactic discs, spiral arms have been modeled as one dimensional filaments whose gravitational fragmentation leads to the formation of giant molecular clouds (GMCs) or star-forming clumps (Inoue & Yoshida, 2018). Within individual GMCs, Herschel observations of star forming regions reveal a multi-scale network of filamentary structures and dense cores aligned with them like beads on a string (André et al., 2010; Jackson et al., 2010; Arzoumanian et al., 2011; Kirk et al., 2013; Palmeirim et al., 2013). This has led to the suggestion that turbulence-driven formation of filaments in the interstellar medium (ISM) is the first step towards core and star-formation (Molinari et al., 2010; André et al., 2010; André et al., 2014), a connection which had been speculated for some time (e.g. Schneider & Elmegreen, 1979; Larson, 1985). In this scenario, the densest filaments with widths of order $\sim 0.1\,{\rm pc}$ (Arzoumanian et al., 2011; Hennebelle & André, 2013) collapse due to GI and lead to the formation of dense cores where star-formation occurs. Simulations of molecular clouds in the ISM reveal similar multi-scale filamentary structures, arising from a variety of mechanisms such as turbulence, gravitational collapse of larger structures, thermal instabilities, or colliding flows (e.g. Padoan et al., 2001; Banerjee et al., 2009; Gómez & Vázquez-Semadeni, 2014; Moeckel & Burkert, 2015; Smith et al., 2016).

Studies of the structure and stability of self-gravitating filaments have a long history, mostly in the context of star-formation in ISM filaments. Early analytic work investigated the stability of an infinite incompressible cylinder with and without an axial magnetic field (Chandrasekhar & Fermi, 1953), a compressible yet still homogeneous infinite cylinder (Ostriker, 1964b), a homogeneous stream of finite radius (Mikhaǐlovskiǐ & Fridman, 1972; Fridman & Poliachenko, 1984), and a uniformly rotating isothermal cylinder (Hansen et al., 1976). Hydrostatic equilibrium of a self-gravitating isothermal cylinder is only possible if its mass per unit length (hereafter line-mass) is less than a critical value which depends only on its temperature (Ostriker, 1964a; see eq. (3) below). For non-isothermal filaments, the critical line-mass is similar (§3.1). Filaments with line-mass larger than the critical value must collapse radially. For line-masses smaller than the critical value, a hydrostatic solution exists, but is unstable to long wavelength axisymmetric perturbations. The fastest growing wavelength is roughly eight times radius of the stream, $\lambda\sim 8R_{\rm s}$ (Nagasawa, 1987, hereafter N87), resulting in stream fragmentation as described in more detail below. A collapsing filament with a line mass slightly exceeding the critical value, as may eventually be the case for a filament growing via radial accretion, is also unstable to axisymmetric perturbations and will fragment at a similar wavelength to the hydrostatic case (Inutsuka & Miyama, 1992). Both cases eventually lead to the formation of bound clumps with masses of order the local Jeans mass (Clarke et al., 2016; Clarke et al., 2017). However, if the line-mass greatly exceeds the critical value the filament collapses towards its axis without fragmenting (Inutsuka & Miyama, 1992). On scales smaller than the filament radius, the local stability criterion reduces to the classical Jeans criterion, even in the presence of rotation (Freundlich et al., 2014). This implies that such local collapse is only possible if the filament is larger than its Jeans length.

N87 studied the stability of a self-gravitating isothermal cylinder with line-mass below the critical value, pressure confined by a low density external medium. He found that the system is always unstable to long-wavelength axisymmetric perturbations even at low values of the line-mass. Similarly, Hunter et al. (1998, hereafter H98) found that a self-gravitating cylinder which is pressure confined by an external medium, with a density discontinuity at the boundary, is always unstable to long-wavelength axisymmetric perturbations. These results are contrary to the spherical case, where a hydrostatic sphere with mass below the critical Bonner-Ebert mass (Ebert, 1955; Bonnor, 1956) is stable against gravitational collapse. We elaborate further on these two studies in §2.1.

In addition to GI, cylindrical streams or jets are susceptible to Kelvin-Helmholtz Instability (KHI) whenever there is a shearing motion between the stream and its surroundings. Numerous authors have studied KHI in cylinders, typically focusing on light or equidense jets meant to represent protostellar or AGN jets (e.g. Birkinshaw, 1984; Payne & Cohn, 1985; Hardee et al., 1995; Bassett & Woodward, 1995; Bodo et al., 1998; Bogey et al., 2011). Several authors have also addressed the effects of magnetic fields and/or radiative cooling on KHI in cylindrical jets (Ferrari et al., 1981; Massaglia et al., 1992; Micono et al., 2000; Xu et al., 2000). However, none of the aforementioned studies accounted for the self-gravity of the gas, as this is expected to be negligible for the systems being considered, namely jets from young stars or AGN. It has also been noted that tidally disrupted streams, resulting from stars tidally destroyed by black holes, may also experience KHI (Bonnerot et al., 2016). Recently, in a series of several papers, Mandelker et al. (2016); Padnos et al. (2018); and Mandelker et al. (2019) (hereafter M16, P18 and M19, respectively) presented a detailed study of KHI, without self-gravity or radiative cooling, in a dense supersonic cylinder representing the cold circumgalactic streams feeding high redshift galaxies. These can be up to 100 times denser than their surroundings. They found that KHI can be important in the evolution of such streams, leading to significant deceleration and energy dissipation, and in certain cases to total stream disruption in the CGM. We elaborate further on these studies in §2.2.

Clearly, extensive work has been done studying separately the effects of GI and of KHI in filaments and streams. While the evolution of KHI in a self-gravitating fluid has been studied in planar (Hunter et al., 1997, hereafter H97) and spherical geometry (Murray et al., 1993, hereafter M93), we are unaware of any such work in cylindrical geometry. Since the evolution of KHI in cylindrical geometry is qualitatively different than in planar geometry (M19, and references therein), while GI in cylinders is qualitatively different than in spheres (e.g. N87; H98), it is worth explicitly studying the combined effects of KHI and self-gravity in cylindrical systems, which is the focus of this paper.

This has important astrophysical implications as well, as there are several filamentary systems where both effects are likely to be important. For instance, it has been shown that the cold circumgalactic streams are likely gravitationally unstable in the inner haloes of massive galaxies at high redshift, potentially resulting in star formation and even globular cluster formation along the streams in the CGM (Mandelker et al., 2018). This may explain recent ALMA observations of dense star-forming gas at distances of tens of $\,{\rm kpc}$ away from a massive galaxy at $z\sim 3.5$, which does not appear to be associated with the galaxy or any of its satellites (Ginolfi et al., 2017). Additionally, filaments in GMCs in the ISM occasionally exhibit shearing flows with respect to their background (Hily-Blant & Falgarone, 2009; Federrath et al., 2016; Kruijssen et al., 2019), suggesting that KHI may be important in their evolution.

The rest of this paper is organised as follows. In §2, we review the current theoretical understanding of GI and KHI in pressure-confined cylinders, and present predictions for how the two may behave in unison. In §3, we describe a suite of numerical simulations used to study GI and KHI in cylinders. In §4 we present the results of our numerical analysis and compare these to our analytical predictions. In §5 we discuss our results and their astrophysical applications, present caveats to our analysis and outline future work. Finally, we summarise our main conclusions in §6.

2 Theory of instabilities

In this section we briefly review the existing theory of GI (§2.1) and KHI (§2.2) in pressure confined cylinders. We then make new predictions for how the two effects may be combined in cylindrical systems (§2.3, to be tested using numerical simulations in §4), and compare these to previous results of a combined analysis in spherical systems (§2.4).

2.1 Gravitational instability

We focus here on the results of N87 and H98, as these are the most relevant for our current analysis. These studies both focus on the stability of a self-gravitating cylinder with finite radius and line-mass below the critical value for hydrostatic equilibrium, pressure confined by a uniform external medium.

N87 consider an isothermal cylinder initially in hydrostatic equilibrium, with the density profile

[TABLE]

(Ostriker, 1964a). $\rho_{\rm c}$ is the central density of the cylinder, $H$ is its scale height, $c_{\rm s}$ is the isothermal sound speed, and $G$ is the gravitation constant. The line-mass of such a cylinder out to radius $R_{\rm s}$ is

[TABLE]

For $R_{\rm s}=\infty$, this yields the critical line-mass for hydrostatic equilibrium (Ostriker, 1964a),

[TABLE]

An equilibrium initial condition is only possible for $\Lambda\leq\Lambda_{\rm cr,iso}$. For a cylinder truncated at a finite radius $R_{\rm s}$, the density and line mass profiles at $r<R_{\rm s}$ are still given by eqs. (1) and (2). Thus, the ratio of the cylinder’s line-mass to the critical line-mass is related to the ratio of the cylinder’s radius to its scale height,

[TABLE]

Increasing the central density, $\rho_{\rm c}$, or decreasing the temperature and thus the sound speed, $c_{\rm s}$, reduces the scale height, $H$. For a fixed stream radius, $R_{\rm s}$, this results in an increase of the ratio $\Lambda/\Lambda_{\rm cr,\,iso}$.

In terms of the external pressure confining the truncated cylinder, pressure equilibrium at the boundary dictates that

[TABLE]

Inserting this into eq. (4) yields

[TABLE]

This shows that for a given temperature and external pressure, a cylinder can have any line mass from [math] to $\Lambda_{\rm cr,\,iso}$, by decreasing the central density from $\rho_{\rm c}=P_{\rm ext}/c_{\rm s}^{2}$ to [math]. The critical line-mass therefore does not depend on the external pressure. This is fundamentally different from the spherical case where the maximal mass for which a hydrostatic equilibrium solution exists depends on the external pressure. This is the Bonnor-Ebert mass,

[TABLE]

(Ebert, 1955; Bonnor, 1956). For further comparison of the structure and properties of self-gravitating cylinders and spheres confined by external pressure, see Fischera & Martin (2012). For the remainder of our analysis we will use the scale-height, $H$, and the stream radius, $R_{\rm s}$, rather than the external pressure.

N87 analyzed perturbations about hydrostatic equilibrium in a cylinder with radius $R_{\rm s}$, pressure confined by an external medium with constant pressure and effectively zero density, $\rho_{\rm ext}<<\rho(R_{\rm s})$. The dispersion relation was numerically evaluated for several values of $\Lambda/\Lambda_{\rm cr,\,iso}$. All cases were found to be stable to non-axisymmetric modes. For axisymmetric modes, the system was found to be unstable at long wavelengths, with longitudinal wavenumber $k<k_{\rm cr}$. The system attains a maximal growth rate, $\omega_{\rm max}$, at a finite wavenumber, $k_{\rm max}$, hereafter the fastest growing mode, and then stabilises again at infinite wavelengths, $\omega\rightarrow 0$ as $k\rightarrow 0$. This is unlike the spherical Jeans instability where the growth rate diverges as $k\rightarrow 0$. There is no closed analytic expression for $k_{\rm cr}$, $k_{\rm max}$ or $\omega_{\rm max}$ for the general case, but it is useful to consider two limiting cases.

In the limit $\Lambda\rightarrow\Lambda_{\rm cr,\,iso}$, equivalent to $R_{\rm s}>>H$ (eq. 4), the solution converges to that of an infinite cylinder. In this case, one obtains $k_{cr}\simeq 0.56H^{-1}$, $k_{\rm max}\simeq 0.28H^{-1}$, and $\omega_{\rm max}\simeq 0.60\left(4G\rho_{\rm c}\right)^{1/2}$. For comparison, the free-fall time of a cylinder with average density $<\rho>=\Lambda/(\pi R_{\rm s}^{2})$ is

[TABLE]

For an isothermal cylinder with radius $R_{\rm s}$,

[TABLE]

For $\Lambda=0.90\Lambda_{\rm cr,\,iso}$, we thus have $R_{\rm s}\simeq 8.5H$ and $\omega_{\rm max}/t_{\rm ff}^{-1}\simeq 1.9$. For larger values of $\Lambda$ the ratio $\omega_{\rm max}/t_{\rm ff}^{-1}$ increases.

In the opposite limit, when $\Lambda<<\Lambda_{\rm cr,\,iso}$ or $R_{\rm s}<<H$, the density is roughly constant within $R_{\rm s}$ and the solution converges to that of an incompressible cylinder, first studied by Chandrasekhar & Fermi (1953). The dispersion relation for an incompressible cylinder is given by111Note that there is a minus sign missing from the corresponding equation (4.10) in N87.

[TABLE]

where $I_{\nu}(x)$ and $K_{\nu}(x)$ are modified Bessel functions of the first and second kind of order $\nu$, evaluated at the argument $x=kR_{\rm s}$. This yields $k_{cr}\simeq 1.1R_{\rm s}^{-1}$, $k_{max}\simeq 0.6R_{\rm s}^{-1}$, and $\omega_{\rm max}\simeq 0.4t_{\rm ff}^{-1}$.

To summarise, the shortest unstable wavelength is $\lambda_{\rm cr}=2\pi/k_{\rm cr}\sim 4\pi H$ and $2\pi R_{\rm s}$ in the limits $\Lambda\rightarrow\Lambda_{\rm cr,\,iso}$ and $\Lambda<<\Lambda_{\rm cr,\,iso}$ respectively. In all cases, the most unstable mode occurs at $\lambda_{\rm max}\sim 2\lambda_{\rm cr}$, while $\omega_{\rm max}/t_{\rm ff}^{-1}$ is within a factor $\sim 2$ of unity. Note that since in the latter limit $R_{\rm s}<<H$, we arrive at the somewhat counterintuitive result that for smaller values of the line-mass the shortest and most unstable wavelengths are much shorter. As noted by N87, the instability manifests itself in different ways in these two limits. For large values of the line-mass the system is unstable to body-modes which are maximal near the stream axis and are similar to the classic Jeans instability. On the other hand, for small values of the line-mass the instability is dominated by surface modes, which are maximal near the stream interface and lead to its deformation. In the non-linear regime, these two modes of instability lead to different shapes and orientations of collapsed clumps within the stream (Heigl et al., 2018b).

H98 generalised this analysis by allowing for a finite background density, $\rho_{\rm b}$, confining the stream. However, they assumed a constant stream density, $\rho_{\rm s}$, rather than an isothermal profile. Their scenario is thus analogous to the limit $\Lambda<<\Lambda_{\rm cr,\,iso}$ from N87. H98 derive the following dispersion relation222This is equivalent to equation (68) from H98 using the identity $I_{0}(x)K_{1}(x)+I_{1}(x)K_{0}(x)=1/x$.

[TABLE]

where $\bar{\rho}=0.5(\rho_{\rm s}+\rho_{\rm b})$, and $I_{\nu}(x)$ and $K_{\nu}(x)$ are again modified Bessel functions with $x=kR_{s}$. This converges to eq. (10) in the limit $\rho_{\rm b}\rightarrow 0$.

From eq. (11), the condition for instability is

[TABLE]

where $\delta=\rho_{\rm s}/\rho_{\rm b}$ is the density contrast between the stream and the background. If $\delta<1$, such that the background is denser than the stream, the system is unstable at all wavelengths due to Rayleigh-Taylor instability (RTI). If $\delta>1$, such that the stream is denser than the background, the system is unstable at long wavelengths, i.e. small values of the argument of the Bessel functions on the left-hand side of eq. (12), $x=kR_{\rm s}$. Furthermore, H98 find that the instability always manifests itself as a surface mode, leading to the deformation of the stream-background interface, similar to the conclusion of N87 for the low line-mass case. For $\delta\rightarrow\infty$, corresponding to $\rho_{\rm b}\rightarrow 0$, the system is unstable for $k<k_{\rm cr}\simeq 1.07R_{\rm s}^{-1}$, as for eq. (10). For $\delta=1$ such that there is no density discontinuity at the interface, $k_{\rm cr}=0$ and the system is stable for all finite wavelengths. This highlights the fact that this is an interface instability, caused by a density discontinuity between the stream and the background. For $\delta=4,~{}10,~{}100$, we have $k_{\rm cr}R_{\rm s}\simeq 0.79,~{}0.96,~{}1.06$. The maximal growth rate for these cases is $\omega_{\rm max}/t_{\rm ff}^{-1}\simeq 0.26,~{}0.36,~{}0.43$.

H98 also note that in the case of a dense sphere pressure confined by a lower density background, the analogous surface mode is always stable, in agreement with the known fact that spheres less massive than the Bonner-Ebert mass are stable. However, in planar geometry, such as a dense slab pressure confined by a lower density background, a similar surface instability exists above a critical wavelength (H97).

The GI surface modes can be thought of as RTI analogues, induced by the self-gravity of the fluid rather than by an external gravitational field. An intuitive explanation was offered by H97 for the planar case, and can be adapted to cylindrical geometry as follows. Consider a dense cylinder with constant density $\rho_{\rm s}$ pressure confined by a background medium with constant density $\rho_{\rm b}<\rho_{\rm s}$. Such a system is stable to classical RTI. Now consider an axisymmetric perturbation to the interface of the cylinder with longitudinal wavelength $\lambda$. In some region, say $0<z<\lambda/2$, there is an outward distortion of the interface, $\xi(z)$, which results in a mass excess just outside the original interface, proportional to $(\rho_{\rm s}-\rho_{\rm b})\xi(z)$. Through Poisson’s equation, this leads to a more negative gravitational potential in this region, resulting in a perturbation $\Phi_{1}<0$ to the initial potential. As the fluid is incompressible and at rest, Bernoulli’s equation tells us that $P+\rho\Phi={\rm const}$ along any streamline in either fluid, where $P$ is the pressure. In the incompressible limit, where $\rho={\rm const}$ in each fluid, this implies that $P_{\rm 1,s}=-\rho_{\rm s}\Phi_{\rm 1}$ and $P_{\rm 1,b}=-\rho_{\rm b}\Phi_{\rm 1}$, where $P_{\rm 1,s}$ and $P_{\rm 1,b}$ are the perturbations to the pressure in the stream and the background respectively, on either side of the interface. Since $\rho_{\rm s}>\rho_{\rm b}$ and $\Phi_{1}<0$, we have that $P_{\rm 1,s}>P_{\rm 1,b}$, so the pressure in the stream just inside the interface is larger than the pressure in the background just outside the interface, causing the perturbation to continue growing.

This instability only manifests at long wavelengths, when the mass excess leading to the perturbation of the potential is large enough to overcome the stabilizing effect of RT modes induced by the unperturbed potential. As noted above, the shortest unstable wavelength for cylinders is $\sim 2\pi R_{\rm s}\sim 6.3R_{\rm s}$ (eq. 12). By contrast, the longest available wavelength on the surface of a sphere corresponds to the $l=2$ spherical harmonic, since the $l=0$ mode represents global expansion or contraction of the sphere while the $l=1$ mode represents a rigid displacement. The wavenumber associated with the $l=2$ mode is $k=[l(l+1)]^{1/2}/R_{\rm s}\sim 2.5/R_{\rm s}$, corresponding to a wavelength of $\lambda\sim 2.6R_{\rm s}$. This is too short for GI surface modes to overcome RT stabilization, which is why there are no GI surface modes for spherical systems (H98).

2.2 KH Instability

KHI arises from shearing motion between the interfaces of two fluids, leading to efficient mixing and smoothing out the initial contact discontinuity. We focus here on the recent results of M19, who analysed the non-linear evolution of KHI in a dense 3d cylinder streaming through a static background, expanding on earlier work by M16 and P18. The system is characterised by two dimensionless parameters, the Mach number of the stream velocity with respect to the background sound speed, $M_{\rm b}=V_{\rm s}/c_{\rm b}$, and the density contrast of the stream and the background, $\delta=\rho_{\rm s}/\rho_{\rm b}$. M19 analytically derived timescales for the non-linear mixing of the two fluids and eventual disruption of the stream, as well as for stream deceleration and the loss of bulk kinetic energy, as a function of these two parameters.

We begin by noting that, similar to the dichotomy between surface modes and body modes in GI (N87), there are two modes of KHI. The nature of the instability depends primarily on the ratio of the stream velocity to the sum of the two sound speeds,

[TABLE]

If $M_{\rm tot}<1$, the instability is dominated by surface modes. These are concentrated at the interface between the fluids, and lead to the growth of a shear layer which expands into both fluids. Within the expanding shear layer a highly turbulent medium develops, efficiently mixing the two fluids. Surface modes can have any longitudinal wavenumber333So long as the wavelength, $\lambda=2\pi/k$, is larger than the width of the transition region between the two fluids., $k$, and any azimuthal wavenumber, $m$, representing the number of azimuthal nodes along the stream-background interface. $m=0$ corresponds to axisymmetric perturbations, $m=1$ to helical perturbations, and $m\geq 2$ to more complicated fluting modes. Low order $m$ modes with wavelengths of order $R_{\rm s}$ dominate the early non-linear evolution of the instability, as their eddies reach the largest amplitudes before they break, but the shear layer between the fluids quickly develops into a highly turbulent mixing zone with no discernible symmetry.

The shear layer separating the fluids expands self-similarly through vortex mergers. Independent of the initial perturbation spectrum, the width of the shear layer, $h$, evolves as

[TABLE]

where $\alpha$ is a dimensionless growth rate that depends primarily on $M_{\rm tot}$, and is typically in the range $\alpha\sim 0.05-0.25$ (P18; M19).

The shear layer penetrates asymmetrically into the stream and background due to their different densities. The penetration depth of the shear layer in either medium can be derived from conservation of mass and momentum in the shear layer, and are given by (P18; M19):

[TABLE]

Stream disruption occurs when the shear layer encompasses the entire stream, namely when $h_{\rm s}=R_{\rm s}$. This occurs at time

[TABLE]

The contact discontinuity effectively disappears before the stream is completely disrupted, once the full width of the shear layer is of order the stream radius, namely $h=R_{\rm s}$. This occurs at time

[TABLE]

As the shear layer expands into the background, it entrains background mass. This causes the stream to decelerate as its initial momentum is distributed over more mass. As shown by M19, the stream velocity as a function of time is well fit by

[TABLE]

where $V_{\rm s,0}$ is the initial velocity of the stream, and

[TABLE]

is the time when the background mass entrained in the shear layer equals the initial stream mass, such that momentum conservation implies the velocity is half its initial value.

An empirical expression for the dimensionless shear layer growth rate, $\alpha$, was proposed by Dimotakis (1991),

[TABLE]

M19 found eq. (20) to be a good fit to shear layer growth in simulations of 2d slabs, regardless of whether one measures $h$, $h_{\rm s}$, or $h_{\rm b}$. However, they found that $h_{\rm s}$ expanded more rapidly in 3d cylinders due to an enhanced eddy interaction rate near the stream axis. This yielded $\alpha$ values $\sim 50\%$ larger than eq. (20) when measuring $h_{\rm s}$ and using eq. (15). On the other hand, $h_{\rm b}$ was found to expand at a similar rate in 2d and 3d so long as h_{\rm b}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}2R_{\rm s}. Since the shear layer width is dominated by $h_{\rm b}$ for $\delta>1$, we use eq. (20) together with eq. (17) to evaluate the time when the contact discontinuity is destroyed.

Once h_{\rm b}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}2R_{\rm s}, its growth rate is reduced by roughly half, due to a turbulent cascade to small scales which removes energy from the largest eddies driving the expansion. For $\delta>8$, this occurs before the stream reaches half its initial velocity (eqs. 18-20). M19 found that in these cases, a good fit to the velocity evolution of streams can be obtained simply by using $0.5\alpha$ in eq. (19) with $\alpha$ taken from eq. (20).

When $M_{\rm tot}>1$, surface modes of low azimuthal order (low values of $m$) stabilise444The formal condition for stabilization of $m=0,\,1$ surface modes is $M_{\rm b}>(1+\delta^{-1/3})^{3/2}$, similar to $M_{\rm tot}>1$.. The nature of the instability then depends on the width of the initial transition region between the fluids (which is likely set by transport processes such as viscosity and thermal conduction). If this is relatively narrow, the instability becomes dominated by high-$m$ surface modes, and the above description, summarised in eqs. (14)-(19), remains valid, with $\alpha\sim 0.05$ according to eq. (20). However, if the initial transition region is wide, of order \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}0.25R_{\rm s} or larger, high-$m$ surface modes are also stable and the instability becomes dominated by body modes. These do not result in shear layer growth but rather in the global deformation of the stream into a helical, $m=1$, shape with a characteristic wavelength of $\sim 10R_{\rm s}$ and an amplitude of \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}R_{\rm s}. The timescale for this to occur depends on the initial perturbation amplitude and spectrum, though it is almost always longer than the timescale for stream disruption by surface modes when these are unstable. Following the formation of the sinusoid, small scale turbulence develops near its peaks and leads to stream disruption within roughly one stream sound crossing time. Interestingly M19 find that eqs. (18)-(19) are a good description of stream deceleration due to body modes as well, despite the different processes involved.

We will hereafter ignore KHI body modes, and assume that KHI is dominated by surface modes of some order $m$ for all Mach numbers. If KHI surface modes are suppressed by a large initial transition region, then GI surface modes will also likely be suppressed, based on the analysis of H97 and H98.

2.3 Combined treatment

We now wish to combine the above two processes, and discuss the evolution of a pressure-confined self-gravitating cylinder undergoing KHI. In addition to $M_{\rm b}$ and $\delta$, a third parameter is required to describe such a system, namely the line-mass of the cylinder in units of the critical line-mass for hydrostatic equilibrium, $\mu\equiv\Lambda/\Lambda_{\rm cr}$. We begin by making the assumption, to be justified below, that any coupling between GI and KHI in the linear regime is relatively small, such that the region of parameter space where each process results in instability is unchanged, and the linear growth rates are only mildly altered. Under this assumption, it is clear from §2.1 and §2.2 that for all values of $(M_{\rm b},\delta,\mu)$, the system is unstable over some wavelength range. We assume that the initial perturbation spectrum spans this range.

GI enhances density contrasts and leads to the formation of long-lived collapsed clumps, while KHI smooths the interface between the fluids and dilutes the mean density of the stream. The question is which process will win. The timescale for GI is the inverse growth rate of the fastest growing mode discussed in §2.1, $t_{\rm max}\equiv\omega_{\rm max}^{-1}$. At low values of $\mu$, GI is dominated by surface modes (N87), which require the presence of a contact discontinuity (H98). Thus, the timescale for KHI to prevent gravitational collapse is $t_{\rm shear}$ (eq. 17), the timescale for nonlinear KHI to destroy the contact discontinuity. On the other hand, for high values of $\mu$, GI is dominated by body modes which are unrelated to the contact discontinuity (N87). In this case, the relevant timescale for KHI to prevent collapse is $t_{\rm dis}$ (eq. 16), the timescale for nonlinear KHI to disrupt the stream itself.

Since $t_{\rm shear}<t_{\rm dis}$ for all $\delta>1$, we distinguish between three regimes. If $t_{\rm max}<t_{\rm shear}<t_{\rm dis}$, we expect GI to win and the stream to fragment into long-lived clumps. If $t_{\rm shear}<t_{\rm dis}<t_{\rm max}$, we expect KHI to win and disrupt the stream by mixing it into the background. We hereafter refer to this process as “shredding the stream”. In the intermediate case where $t_{\rm shear}<t_{\rm max}<t_{\rm dis}$, the outcome may depend on the value of $\mu$. If $\mu$ is small, such that GI is dominated by surface modes, then we expect KHI to win and shred the stream since $t_{\rm shear}<t_{\rm max}$. On the other hand, if $\mu$ is large such that GI is dominated by body modes, GI may still win and lead to stream fragmentation and the formation of bound clumps, since $t_{\rm max}<t_{\rm dis}$. However, this is uncertain, since the shear layer will penetrate somewhat into the stream within $t_{\rm max}$, reducing the effective line-mass of the unperturbed (non-turbulent) region. If this is reduced below the threshold for GI body modes to be effective, KHI may still win and suppress clump formation.

Since $t_{\rm max}\propto\rho_{\rm c}^{-1/2}\propto\mu^{-1/2}$, as $\mu$ is increased at fixed $(M_{\rm b},\delta)$, $t_{\rm max}$ decreases while $t_{\rm shear}$ (and $t_{\rm dis}$) remain constant. Thus, for each $(M_{\rm b},\delta)$ there exists a critical value of $\mu$, $\mu_{\rm cr}\equiv\mu_{\rm cr}(M_{\rm b},\delta)$, such that $t_{\rm max}<t_{\rm shear}$ for $\mu>\mu_{\rm cr}$ (see Fig. 5 in §4.1 below). Therefore, GI will win and lead to stream fragmentation and clump formation whenever $\mu>\mu_{\rm cr}$. If $\mu_{\rm cr}$ is small enough to be in the regime where GI is dominated by surface modes, then KHI will win and shred the stream for $\mu<\mu_{\rm cr}$. On the other hand, if $\mu_{\rm cr}$ is in the regime where GI is dominated by body modes, the fate of the stream at $\mu<\mu_{\rm cr}$ depends on the ratio of $t_{\rm max}$ to $t_{\rm dis}$.

At first glance, it may seem inconsistent to compare a linear timescale for GI, $t_{\rm max}=\omega_{\rm max}^{-1}$, to a nonlinear timescale for KHI, $t_{\rm shear}$ or $t_{\rm dis}$. While $t_{\rm max}$ is formally the timescale for the growth of linear perturbations, once density perturbations grow the free-fall times become ever shorter and the collapse accelerates. Full collapse is thus dominated by the linear growth time. On the other hand, KHI tends to saturate following the linear phase, because it is driven by the presence of a contact discontinuity which is destroyed by the instability. Continued growth in the nonlinear regime is dominated by the merger of eddies within the shear layer on timescales of $t_{\rm shear}$ and $t_{\rm dis}$, as described in §2.2.

The above discussion notwithstanding, one may ask whether density fluctuations within the stream induced by KHI can trigger local gravitational collapse when $\mu<\mu_{\rm cr}$. Note that this is different than the global fragmentation of the stream induced by GI. Such local collapse can occur in filaments on scales larger than the spherical Jeans length, $\lambda_{\rm J}=[\pi c_{\rm s}^{2}/(G\rho)]^{1/2}$, but smaller than the stream radius, $R_{\rm s}$ (Freundlich et al., 2014). This implies that this is only possible if $\lambda_{\rm J}<R_{\rm s}$. A lower limit to the Jeans length is obtained by inserting $\rho=\rho_{\rm c}$, the density along the stream axis. This yields $\lambda_{\rm J}=2\pi H$, with $H$ given by eq. (1). The condition that $\lambda_{\rm J}<R_{\rm s}$ thus implies that $R_{\rm s}>>H$, so GI is dominated by body modes (N87). We conclude that KHI induced density fluctuations can only trigger local gravitational collapse if $\mu<\mu_{\rm cr}$ but GI is still dominated by body modes.

We must now justify our initial ansatz that the linear coupling between GI and KHI does not fundamentally alter the instability region of parameter space. We rely here on the analysis of H97, who derived the dispersion relation of a self-gravitating system undergoing KHI in the vortex sheet limit, i.e. two semi-infinite fluids separated by a single, planar interface. In their derivation they made the simplifying assumption that the gravitational field in the unperturbed system was weak compared to the perturbed forces induced by both pressure and potential perturbations. This is equivalent to assuming that the wavelengths are much shorter than the gravitational scale-height of the unperturbed system, which itself is equivalent to assuming constant density and pressure in both fluids. The resulting dispersion relation contains terms associated with KHI, RTI, and surface mode GI. We refer the reader to H97 for the expression and its derivation. Relevant to our discussion is the fact that the coupling between self-gravity and shearing motions does not modify the stability region of the system, only mildly affects the linear growth rates of KH modes at short to intermediate wavelengths, and does not suppress GI surface modes at long wavelengths. Deriving an analogous dispersion relation for cylinders is beyond the scope of this paper. Rather, we assume that the same conclusions hold for cylindrical systems, in particular because KHI in cylinders is even more unstable than for planar vortex sheets (M16; M19). The validity of this assumption and our subsequent analysis will be tested with numerical simulations in §4.

2.4 Comparison to the Spherical Case

It is worth comparing our analysis to that of M93, who addressed the question of when self-gravity would prevent KHI from disrupting a cold, dense spherical cloud moving through a hot, dilute background. They assumed that the cloud was pressure confined by the background fluid, and that its mass was less than the Bonnor-Ebert mass, making it gravitationally stable and in hydrostatic equilibrium. In this case, unlike for self-gravitating cylinders, there is no GI, and the only effect of the self-gravity is to induce RT modes at the cloud surface. Since the cloud is denser than the background, these RT modes can counteract the KHI and stabilise the system, due to the restoring buoyancy force. They showed this by considering the combined dispersion relation of KHI and RTI in the incompressible limit,

[TABLE]

where $V$ is the velocity of the cloud in the static background and $g$ is the gravitational acceleration at its surface. This implies that KHI is stable for all wavelengths greater than

[TABLE]

M93 then assumed that KHI would only disrupt the cloud if $\lambda_{\rm max}>R_{\rm cl}$, the cloud radius. This was based on the assumption that KHI surface modes saturate at an amplitude comparable to their wavelength, thus neglecting the subsequent shear layer growth. This assumption together with $g=GM_{\rm cl}/R_{\rm cl}^{2}$ and $M_{\rm cl}=(4\pi/3)\rho_{\rm cl}R_{\rm cl}^{3}$ results in a minimum mass for self-gravity to stabilise the sphere against KHI. For velocities of order the background sound speed, the critical mass is of order the Bonnor-Ebert mass, $M_{\rm BE}$. Such a system is thus always unstable, either to KHI at $M_{\rm cl}<M_{\rm BE}$ or to global gravitational collapse at $M_{\rm cl}>M_{\rm BE}$.

Our main prediction for the cylindrical case is qualitatively similar. We predict that a self-gravitating stream will always be unstable either to KHI at $\mu<\mu_{\rm cr}$ or to GI at $\mu>\mu_{\rm cr}$, depending on whether the timescale for GI, $t_{\rm max}$, is longer or shorter than the timescale for KHI to destroy the contact discontinuity, $t_{\rm shear}$, and/or the stream itself, $t_{\rm dis}$. However, unlike M93, we do not rely on a similar criterion of gravity stabilizing wavelengths longer than $R_{\rm s}$. First of all, unlike in spherical systems, self-gravity actually destabilises cylinders at long wavelengths (N87; H98; §2.1). Furthermore, even if KHI is stable for wavelengths longer than $R_{\rm s}$ in the linear regime, it can still lead to stream disruption in the nonlinear regime by shear layer growth caused by initially shorter wavelength perturbations.

3 Numerical Methods

In this section we describe the details of our simulation code and setup, as well as our analysis method. We use the Eulerian AMR code RAMSES (Teyssier, 2002), with a piecewise-linear reconstruction using the MonCen slope limiter (van Leer, 1977), an HLLC approximate Riemann solver (Toro et al., 1994), and a multi-grid Poisson solver.

3.1 Hydrostatic Cylinders

Unlike the isothermal cylinder described in §2.1, there is no closed analytic expression for the density profile of an isentropic cylinder in hydrostatic equilibrium, so this must be evaluated numerically. We briefly review here how this is done, beginning with the equilibrium solution of an isolated cylinder following Ostriker (1964a). The equation of hydrostatic equilibrium,

[TABLE]

is solved together with Poisson’s equation

[TABLE]

and an isentropic equation of state (EoS),

[TABLE]

where we assumed $K$ to be constant and the adiabatic index of ideal monoatomic gas, $\gamma=5/3$, throughout. These equations can be combined to yield

[TABLE]

with the boundary conditions

[TABLE]

Eqs. (26)-(27) can be cast into unitless form by defining $y=\rho/\rho_{\rm c}$ and $x=r/H$, with

[TABLE]

the scale radius of the cylinder, where $c_{\rm s,0}^{2}=\gamma P_{\rm c}/\rho_{\rm c}=\gamma K\rho_{\rm c}^{\gamma-1}$ is the sound speed along the filament axis, with $P_{\rm c}=P(r=0)$ the pressure along the filament axis. The resulting equation is

[TABLE]

Analytic solutions exist only for $\gamma=1$ (isothermal cylinder), $\gamma=2$, and $\gamma=\infty$ (incompressible cylinder) (Ostriker, 1964a). For other values of $\gamma$ eq. (29) must be solved numerically.

While the isothermal cylinder discussed in §2.1 extends to $r=\infty$, all cases with $\gamma>1$ have a finite radius, $R_{\rm equ}$, defined as the radius where the density profile first reaches $\rho=0$ (Ostriker, 1964a). We can thus generalise the notion introduced in §2.1 of a critical line-mass above which hydrostatic equilibrium is not possible

[TABLE]

where $y(x)$ is the solution to eq. (29). The factor $a$ on the right-hand-side of eq. (30) depends on the EoS. For $\gamma=5/3$, $R_{\rm equ}\simeq 2.648H$, the half-mass radius is $R_{\rm 1/2}\simeq 1.168H$, and $a\simeq 0.796$. For comparison, an isothermal cylinder has $R_{\rm 1/2}\simeq 2.828H$, with $H$ defined in eq. (1), and $a=2$ (eq. 3). In Fig. 1 we show the normalised equilibrium density and line-mass profiles of an isolated, isentropic, $\gamma=5/3$ cylinder.

Equilibrium profiles with $\Lambda<\Lambda_{\rm cr}$ can be constructed for cylinders pressure confined by an external medium and truncated at some radius $R_{\rm s}<R_{\rm equ}$. In Fig. 2 we show the stream radius, $R_{\rm s}/H$, as a function of $\mu=\Lambda/\Lambda_{\rm cr}$. For $\mu=0,\,1$ we have $R_{\rm s}=0,\,R_{\rm equ}$ respectively. For $\mu=0.5$ we have $R_{\rm s}=R_{1/2}\simeq 1.17H$. We adopt model units where $G=\rho_{\rm c}=1$ and $R_{\rm s}=1/32$. For a given value of $\mu$, we can obtain $H$ in model units from Fig. 2 and then eq. (28) can be used to obtain $c_{\rm s,0}=(8\pi/3)^{1/2}H$ and $P_{\rm c}=K_{\rm s}=3c_{\rm s,0}^{2}/5$. Note that the stream and the background fluid have different entropy, and hence different values of $K$.

In addition to $\mu$, the system is defined by

[TABLE]

the ratio of the density along the stream axis to the background density just outside the stream. For a given $\mu$ and $\delta_{\rm c}$ we may evaluate the density contrast between the stream and background on either side of the interface,

[TABLE]

We show the ratio $\rho(R_{\rm s}^{-})/\rho_{\rm c}=\delta/\delta_{\rm c}$ as a function of $\mu$ in Fig. 2. For $\mu=0.1,\,0.5,\,0.9$ we have $\delta/\delta_{\rm c}\simeq 0.92,\,0.58,\,0.18$ respectively.

To construct equilibrium profiles for pressure confined cylinders with given values of $\mu$ and $\delta_{\rm c}$, we first evaluate $R_{\rm s}/H$ and $\delta$ from Fig. 2. We then solve eq. (29) separately for $r<R_{\rm s}$ and $r>R_{\rm s}$. For $r<R_{\rm s}$, the boundary conditions are $y(0)=1$, and ${\rm d}y/{\rm d}x|_{0}=0$, and the profile is unchanged from the isolated cylinder. For $r>R_{\rm s}$, the boundary conditions are given in terms of the pressure, rather than the density. Specifically, the pressure is continuous at the interface, $P(R_{\rm s}^{-})=P(R_{\rm s}^{+})$, while the pressure gradient is discontinuous, with

[TABLE]

which follows from eq. (23).

Fig. 3 shows the resulting density and pressure profiles for $\mu=0.1$ and $0.9$. For $\mu=0.1$ the density and pressure are nearly constant in either medium, while for $\mu=0.9$ there are strong gradients within the stream.

3.2 Initial Conditions

3.2.1 Simulation Domain & Boundary Conditions

The simulation domain is a cube of side $L=1$, extending from $-0.5$ to $0.5$ in all directions. We hereafter adopt the standard cylindrical coordinates, $(r,\varphi,z)$. The axis of our cylindrical stream is placed along the $z$ axis, at $r=0$, and we adopt a stream radius of $R_{\rm s}=1/32$. The stream fluid occupies the region $r<R_{\rm s}$ while the background fluid occupies the rest of the domain. The equation of state (EoS) of both fluids is that of an ideal monoatomic gas with adiabatic index $\gamma=5/3$.

We use periodic boundary conditions at $z=\pm 0.5$, and zero force boundary conditions, often called outflow boundary conditions, at $x=\pm 0.5$ and $y=\pm 0.5$, such that gas crossing the boundary is lost from the simulation domain. At these boundaries, the gradients of density and velocity are set to 0, while the pressure gradient is taken from the hydrostatic profile computed following §3.1. The potential at the boundary is set to be that at the outer edge of an isolated and infinitely long cylinder with total mass $M$, equal to the total mass in the simulation domain, $\Phi(r)=2G(M/L){\rm ln}(r)$ with $r=(x^{2}+y^{2})^{1/2}$ on the boundary. We note that this does not produce perfect equilibrium due to fitting a cylindrical profile in a cubic box. However, we find that our configuration is extremely stable in simulations with no initial perturbations and no shear flow, exhibiting \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}3\% change in the density and pressure profiles after 4 stream free-fall times.

3.2.2 Smoothing the Discontinuity

As noted by many previous studies of KHI, the presence of a sharp discontinuity at the interface of two fluids leads to numerical perturbations on the grid scale. These grow faster than the intended perturbations in the linear regime, and may dominate the instability at late times depending on their amplitude. Furthermore, since smaller scales grow more rapidly in the linear regime, these numerical perturbations become more severe as the resolution is increased, preventing convergence of the solution. To alleviate this issue, we smooth the velocity and density around the interfaces using the ramp function proposed by Robertson et al. (2010), also used by M16, P18, and M19. Specifically, we normalise each quantity in the stream and the background by its value at $R_{\rm s}$, denoted $f_{\rm s}$ and $f_{\rm b}$ respectively. We then smooth between these values using

[TABLE]

[TABLE]

and multiply the normalised profiles in either medium by $f(r)$. The parameter $\sigma$ determines the width of the transition zone. The function $\mathcal{R}(r)$ transitions from $0.05$ to $0.95$ over a full width of $\sim 3\sigma$ in $(r-R_{\rm s})$. We adopt $\sigma=R_{\rm s}/32$ for all of our simulations, which is sufficient to suppress artificial perturbations with small longitudinal wavelength, while still allowing azimuthal modes with m\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}12 to grow (M19).

3.2.3 Perturbations

The stream is initialised with velocity $\vec{v_{s}}=M_{\rm b}c_{\rm b}\hat{z}$, where $c_{b}=[\gamma P(R_{\rm s}^{+})/\rho(R_{\rm s}^{+})]^{1/2}$ is the sound speed at the outer boundary of the stream. The background gas is initialised at rest, with velocity to $\vec{v_{b}}=0$.

We then perturb our simulations with a random realization of periodic perturbations in the radial component of the velocity, $v_{\rm r}=v_{\rm x}{\rm cos}(\varphi)+v_{\rm y}{\rm sin}(\varphi)$, as in M19. In practice, we perturb the Cartesian components of the velocity,

[TABLE]

[TABLE]

The velocity perturbations are localised on the stream-background interface, with a penetration depth set by the parameter $\sigma_{\rm pert}$. We set $\sigma_{\rm pert}=R_{\rm s}/16$ in all of our simulations, as in M19. To comply with periodic boundary conditions, all wavelengths were harmonics of the box length, $k_{j}=2\pi n_{j}$, where $n_{j}$ is an integer, corresponding to a wavelength $\lambda_{j}=1/n_{j}$. In each simulation, we include all wavenumbers in the range $n_{j}=2-64$, corresponding to all available wavelengths in the range $R_{\rm s}/2-16R_{\rm s}$. Each perturbation mode is also assigned a symmetry mode, represented by the index $m_{j}$ in eqs. (36) and (37), and discussed in §2.2. As in M19, we only consider $m=0,1$. For each wavenumber $k_{j}$ we include both an $m=0$ mode and an $m=1$ mode. This results in a total of $N_{\rm pert}=2\times 63=126$ modes per simulation. Each mode is then given a random phase $\phi_{j}\in[0,2\pi)$. The stochastic variability from changing the random phases was extremely small, as shown in P18 and M19. The amplitude of each mode, $v_{0,j}$, was identical, with the rms amplitude normalised to $0.01c_{\rm s}$.

3.2.4 Resolution and Refinement Scheme

We used a statically refined grid with resolution decreasing away from the stream axis. The highest resolution region is ${\rm max}(|x|,|y|)<3R_{\rm s}$, with cell size $\Delta=2^{-10}$. For $R_{\rm s}=1/32$ this corresponds to 64 cells per stream diameter. The cell size increases by a factor of 2 every $3R_{\rm s}$ in the $x$ and $y$-directions, up to a maximal cell size of $2^{-6}$. The resolution is uniform along the $z$ direction, parallel to the stream axis. For uniform density cylinders, KHI surface modes are converged at this resolution (M19). We also ran two cases with a factor 2 higher resolution (128 cells per stream diameter) in order to test convergence of our results for self-gravitating streams. As described in §4.2 and §4.3, we find that the majority of our results are indeed converged.

3.2.5 Simulations Without Self-Gravity

In addition to the simulations of self-gravitating cylinders described above, we performed several simulations without self-gravity for comparison, hereafter our “no-gravity” simulations. In the no-gravity simulations, the boundary conditions at $x=\pm 0.5$ and $y=\pm 0.5$ are simply zero gradients in all fluid variables, including the pressure. These were then initialised with the same density profiles as the corresponding self-gravitating streams, but with constant pressure throughout the simulation domain, since there is no gravitational field. We set the pressure to be the same as the pressure at the stream boundary in the corresponding self-gravity simulations, $P_{\rm no-gravity}(r)=P_{\rm self-gravity}(R_{\rm s})$. This allows us to separate the effects of the density profile from those of self-gravity on the evolution of KHI.

3.3 Tracing the Two Fluids

Following P18 and M19, we use a passive scalar field, $\psi(r,\varphi,z,t)$, to track the growth of the shear layer and the mixing of the two fluids. Initially, $\psi=1$ and [math] at $r<R_{\rm s}$ and $r>R_{\rm s}$ respectively, and is then smoothed using eqs. (34)-(35). During the simulation, $\psi$ is advected with the flow such that the density of stream-fluid in each cell is $\rho_{\rm s}=\psi\rho$, where $\rho$ is the total density in the cell.

The volume-weighted average radial profile of the passive scalar is given by

[TABLE]

The resulting profile is monotonic (neglecting small fluctuations on the grid scale) and can be used to define the edges of the shear layer around the stream interface, $r(\overline{\psi}=\epsilon)$ on the background side and $r(\overline{\psi}=1-\epsilon)$ on the stream side, where $\epsilon$ is an arbitrary threshold. We set $\epsilon=0.04$ following M19, though our results are not strongly dependent on this choice. The background-side thickness of the shear layer is then defined as

[TABLE]

while the stream-side thickness is defined as

[TABLE]

While $h_{\rm b}$ as defined in eq. (39) is always well defined, at late times the perturbed region encompasses the entire stream and $\overline{\psi}(r=0)<1-\epsilon$. In this case, we define $h_{\rm s}=R_{\rm s}$. The total width of the perturbed region is given by $h\equiv h_{\rm b}+h_{\rm s}$.

4 Results

In this section we present the results of our numerical simulations. In §4.1, we examine when the combined evolution of GI and KHI leads to the formation of long-lived clumps or to stream shredding, and compare to our theoretical predictions. In §4.2 and §4.3, we discuss the late time evolution of the system in the cases when KHI and GI dominate, respectively.

4.1 KHI vs GI

As detailed in §2.3, we predict that a dense, self-gravitating filament shearing against a dilute background will either fragment into long-lived, bound clumps due to GI, or disrupt and mix into the background due to KHI, depending on the ratio of their respective timescales. The timescale for GI, $t_{\rm max}(\mu,\delta_{\rm c})$, is well approximated by eq. (11) (see the Appendix §A). The timescales for KHI are $t_{\rm shear}(M_{\rm b},\delta_{\rm c})$ (eq. 17) or $t_{\rm dis}=(1+\delta^{1/2})t_{\rm shear}$ (eq. 16). For given values of $(M_{\rm b},\delta_{\rm c})$ there is a critical line mass ratio, $\mu_{\rm cr}$, such that for $\mu>\mu_{\rm cr}$, $t_{\rm max}<t_{\rm shear}$ and GI will dominate. If $\mu_{\rm cr}$ is small enough to be in the regime where GI is dominated by surface modes, then KHI will dominate for $\mu<\mu_{\rm cr}$. However, if $\mu_{\rm cr}$ is in the regime where GI is dominated by body modes, then the fate of the stream when $\mu<\mu_{\rm cr}$ depends also on the ratio of $t_{\rm max}$ to $t_{\rm dis}$.

Solid curves in Fig. 4 show the ratio $t_{\rm max}/t_{\rm shear}$ as a function of $\mu$ for $(M_{\rm b},\delta_{\rm c})=(1.0,100)$, $(1.0,6.7)$, $(2.5,100)$, and $(6.0,100)$. The corresponding values of $\mu_{\rm cr}$ are $\sim 0.36$, $0.28$, $0.62$, and $0.96$. Note the very weak dependence of $t_{\rm max}/t_{\rm shear}$ on $\delta_{\rm c}$ for $M_{\rm b}=1$, since $t_{\rm max}$ depends weakly on $\delta_{\rm c}$ for \delta_{c}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}4, while $t_{\rm shear}$ depends weakly on $\delta_{\rm c}$ only through $\alpha(M_{\rm tot})$ (eq. 20). The dependence of $t_{\rm max}/t_{\rm shear}$ on $M_{\rm b}$ is much stronger, since $t_{\rm shear}$ decreases roughly linearly with $M_{\rm b}$.

In Fig. 5 we show $\mu_{\rm cr}$ as a function of $M_{\rm b}$ and $\delta_{\rm c}$. The general trend is the same as inferred from Fig. 4, namely $\mu_{\rm cr}$ increases strongly with $M_{\rm b}$ and has only a slight tendency to increase with $\delta_{\rm c}$. The exception is a narrow strip near $M_{\rm b}\sim(1-2)$ where $\mu_{\rm cr}$ decreases with $M_{\rm b}$. In this region, the increase of $t_{\rm shear}$ due to decreasing $\alpha$ is stronger than the decrease in $t_{\rm shear}$ due to increasing $V$, leading to a net increase in $t_{\rm shear}$ with $M_{\rm b}$ and thus a net decrease in $\mu_{\rm cr}$. For density contrasts \delta_{\rm c}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}100, $\mu_{\rm cr}>0.5$ only for supersonic flows with M_{\rm b}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}2.5. This implies that for massive streams, KHI can only overcome GI for highly supersonic flows (recall that the Mach number of the flow with respect to the sound speed in the stream is $\sim\delta^{1/2}M_{\rm b}$). In this regime, KHI is dominated by high-order azimuthal surface modes (see §2.2), which have a short eddy turnover time leading to rapid shear layer growth.

Consider, for example, $\delta_{\rm c}\sim 30$. $\mu_{\rm cr}$ increases from $\mu_{\rm cr}<<1$ at $M_{\rm b}<<1$ towards \mu_{\rm cr}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}0.3 at $M_{\rm b}\sim 0.6$, then decreases to \mu_{\rm cr}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}0.2 at $M_{\rm b}\sim 1.2$, before strongly increasing at $M_{\rm b}>>1$. Thus, as $M_{\rm b}$ is increased from \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}0.2 to \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}2 for $\delta_{\rm c}\sim 30$ and $\mu\sim 0.25$, the stream fluctuates from being dominated by GI, to KHI, to GI, to KHI. The high $M_{\rm b}$ KHI regime is dominated by surface modes with high azimuthal wavenumber. While these modes are always unstable, at lower Mach numbers they tend to be sub-dominant compared to axisymmetric or helical modes, with $m=0,1$ (M19).

To test our predictions, we performed a series of simulations with the same combinations of $(M_{\rm b},\delta_{\rm c})$ as shown in Fig. 4, and different values of $\mu$. For $(M_{\rm b},\delta_{\rm c})=(1.0,100)$, we performed nine simulations spanning the line-mass range $\mu=0.1,\,0.2,\,...,\,0.9$. For the other combinations of $(M_{\rm b},\delta_{\rm c})$, we performed three to four simulations each, with $\mu$ spanning a small region around the predicted $\mu_{\rm cr}$. The full list of simulations is presented in Table 1, along with several relevant parameters. The stream sound crossing time555The sound crossing times listed in Table 1 refer to the self-gravity simulations only. In the no-gravity runs at $r<R_{\rm s}$, $\rho_{\rm no-gravity}(r)=\rho_{\rm gravity}(r)$ while $P_{\rm no-gravity}(r)=P_{\rm gravity}(R_{\rm s})<P_{\rm gravity}(r)$. This results in a lower sound speed at each $r<R_{\rm s}$, and hence a longer sound crossing time., $t_{\rm sc}$, is defined as

[TABLE]

where $c_{\rm s}(r)=(\gamma P(r)/\rho(r))^{1/2}$ is the sound speed at radius $r$.

The markers in Fig. 4 indicate for each of our simulations whether or not the stream has fragmented into long-lived collapsed clumps. To identify such clumps, we examine the PDF of stream fluid density, $\rho_{\rm s}=\psi\rho$. If the density distribution is a result of pure turbulence, we expect it to have a roughly lognormal shape. If, however, the highest density regions have collapsed due to gravity, we expect a break in the PDF at high densities (e.g. Vázquez-Semadeni et al., 2008; Elmegreen, 2011; Hopkins et al., 2012; Federrath & Banerjee, 2015). An example is shown in Fig. 6, where we show the density PDFs for the gravity and no-gravity simulations with $(M_{\rm b},\delta_{\rm c},\mu)=(1.0,100,0.9)$ at $t=8t_{\rm sc}$. While the no-gravity simulation has a unimodal PDF which is roughly lognormal except at the lowest densities, the gravity simulation produces a bi-modal PDF, and we associate all cells with densities larger than the break density, $\rho_{\rm s,th}$, as being in clumps. As discussed in §4.3 below, these clumps are indeed long-lived. If a simulation never exhibits a similar break in the density PDF we determine that this simulation has not formed any clumps. In particular, isolated high density regions produced in no-gravity simulations at late times (see Figs. 7 and 11 below) are not clumps, but rather transient features associated with the high-density part of a turbulent PDF.

All of our simulations with $\mu>\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$ form gravitating clumps, as predicted by our model. Furthermore, for $M_{\rm b}=1,2.5$, when \mu_{\rm cr}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}0.63, streams in simulations with $\mu<\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$ are disrupted by KHI and mixed into the background before forming bound clumps, as predicted by our model. In these cases, GI is dominated by surface modes, so the comparison of $t_{\rm max}$ and $t_{\rm shear}$ is justified. On the other hand, for $M_{\rm b}=6.0$, $\mu_{\rm cr}=0.96$ is in the body mode regime for GI, and our simulation with $\mu=0.9$ fragments into bound clumps, as discussed in §2.3. However, in this same regime we find that streams with $\mu=0.8$ and $0.7$ are disrupted by KHI and do not form bound clumps. So the effect of GI body modes is to lower $\mu_{\rm cr}$ from $\sim 0.96$ to $\sim 0.85$.

Overall, we conclude that our model adequately predicts the fate of streams under the combined effects of KHI and GI when GI surface modes dominate. When GI body modes dominate, the actual value of $\mu_{\rm cr}$ is $\sim 10\%$ lower than our prediction, since the relevant timescale for KHI to prevent clump formation is no longer $t_{\rm shear}$. In the next two sections, we turn to studying the evolution of streams and clumps in the regime where each instability dominates.

4.2 Stream Disruption due to KHI

We here examine the evolution of streams with $\mu<\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$, where KHI dominates over GI and prevents the formation of long-lived collapsed clumps. Specifically, we examine whether the self-gravity of the gas, while unable to completely overcome the KHI, affects its evolution in any way.

Figure 7 shows the evolution of streams with $(M_{\rm b},\delta_{\rm c},\mu)=(1.0,100,0.1)$, with and without self-gravity. At early times, t\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}4t_{\rm sc}, the evolution in the two cases is extremely similar, and the shear layer seems to expand at roughly the same rate, mixing the two fluids and diluting the stream density. At later times, t\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}6t_{\rm sc}, while the expansion of the shear layer into the background continues similarly in both simulations, the penetration into the stream has stalled in the gravity run. The self-gravity of the stream thus seems to partly shield its inner core from mixing and disruption. As we will show below, this is due to restoring buoyancy forces caused by the stream’s gravitational field.

We examine this more quantitatively in Fig. 8, where we compare the evolution of $h_{\rm b}$ and $h_{\rm s}$, the penetration of the shear layer into the background and stream respectively (eq. 15), in gravity and no-gravity simulations with $\mu<\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$. Focusing on the top panel, we see that $h_{\rm b}$ evolves similarly with and without self-gravity, and is consistent with the results of M19. During the first sound crossing time, $h_{\rm b}$ remains roughly constant as the initial velocity perturbations trigger perturbations in the stream-background interface associated with growing eigenmodes of the system. Following this phase, $h_{\rm b}$ grows approximately linearly following eq. (15) until it reaches $h_{\rm b}\sim 2R_{\rm s}$. Up until this point, the gravity and no-gravity runs are nearly indistinguishable. Following this, the growth rate of $h_{\rm b}$ is reduced by roughly half in both cases, as the developing turbulent cascade transfers power from the largest scales driving the expansion to smaller scales (M19). During this phase, the growth rate of $h_{\rm b}$ is reduced in the gravity simulations, by $\sim 25\%$ for $\mu=(0.5-0.6)$ and $\sim 12\%$ for $\mu=(0.1-0.2)$. Overall, we find that the self-gravity of the stream has a relatively minor effect on the growth of $h_{\rm b}$.

On the other hand, as inferred from Fig. 7, there is a qualitative difference in the evolution of $h_{\rm s}$, shown in the bottom panel of Fig. 8. For the first $\sim 2t_{\rm sc}$, until $h_{\rm s}\sim 0.3R_{\rm s}$, the gravity and no-gravity runs evolve similarly. After this, the growth rate in the gravity runs is a factor $\sim 3$ smaller than in the no-gravity runs. In the latter, the shear layer reaches $h_{\rm s}/R_{\rm s}=1$ and consumes the entire stream at $t\sim t_{\rm dis}$ (eq. 16), and the evolution is similar to that seen in M19 (see their figure B1). However, in the runs with self-gravity, $h_{\rm s}\sim(0.3-0.4)R_{\rm s}$ at this time, and does not exceed $\sim 0.5R_{\rm s}$ at $t=10t_{\rm sc}$. This is consistent with the visual impression in Fig. 7, where the density remains high in the interior of the self-gravitating stream even after the non-gravitating stream has been completely diluted. Although the growth rate of $h_{\rm s}$ does not depend on $\mu$, there appears to be a tendency for more penetration for larger $\delta_{\rm c}$.

We propose that the stalling of $h_{\rm s}$ is due to restoring buoyancy forces in the stream interior. This can be seen by considering the Richardson number, ${\rm Ri}=[N_{\rm BV}/({\rm d}u/{\rm d}r)]^{2}$, where ${\rm d}u/{\rm d}r$ is the gradient of longitudinal velocity inside the shear layer, and $N_{\rm BV}$ is the Brunt-Vais$\ddot{{\rm a}}$l$\ddot{{\rm a}}$ frequency,

[TABLE]

with $g(r)$ the magnitude of the gravitational field, $\gamma$ the adiabatic index, $K(r)=P(r)\rho^{-\gamma}(r)$ the entropy profile of the gas. Note that $K$ is piecewise constant in our initial conditions, with a non-zero gradient only at the stream-background interface. However, as the shear layer expands, mixing between the fluids creates a non-zero entropy gradient throughout the shear layer. Had our initial conditions been such that the initial stream was not isentropic, this may have increased $N_{\rm BV}$ and ${\rm Ri}$ in the stream interior.

For a 2d plane-parallel system in a constant external gravitational field, it can be shown that a sufficient (but not necessary) criterion for buoyancy to stabilize the system against shearing induced mixing is that ${\rm Ri}>0.25$ (Miles, 1961; Howard, 1961). While our situation is more complex in that the geometry is cylindrical and the gravitational field is due to self-gravity rather than an external field666To our knowledge, no analogous criterion exists for the stability of self-gravitating flows or for cylindrical flows. Deriving such a criterion is beyond the scope of this paper., we may use this as a benchmark to asses the role of buoyancy in stabilizing the inner stream. In Fig. 9 we show the mass weighed average of ${\rm Ri}$ within the shear layer, $[R_{\rm s}-h_{\rm s}(t)]<r<R_{\rm s}$, as a function of time. In all simulations, ${\rm Ri}<<1$ at early times, and crosses ${\rm Ri}=0.25$ at $t\sim(2-3)t_{\rm sc}$, corresponding to the sharp decline in the growth rate of $h_{\rm s}$. Further growth of ${\rm Ri}$ is rather slow and it does not exceed ${\rm Ri}\sim(0.3-0.4)$. We find very similar behaviour when evaluating ${\rm Ri}$ locally at the inner boundary of the shear layer, $r=[R_{\rm s}-h_{\rm s}(t)]$. This supports our assertion that buoyancy stabilizes the inner stream and slows the growth of $h_{\rm s}$, significantly delaying stream disruption.

In Fig. 10 we show the deceleration of streams in simulations with and without gravity. We show the centre of mass velocity of the stream fluid, i.e. weighted by the passive scalar $\psi$, normalised by its initial value, $V_{\rm i}$, as a function of time normalised by the predicted decelration timescale, $t_{\rm dec}$ (eq. 19). The time axis has been shifted to begin at $t_{0}$, the time when the stream velocity reaches $98\%$ of its initial value. In all cases, $t_{0}\sim t_{\rm sc}$. The gravity and no-gravity simulations behave similarly, and both are well fit by the theoretical prediction (eq. 18). This was expected given the similarity between the evolution of $h_{\rm b}$ in the gravity and no-gravity simulations (Fig. 8), since KHI-induced deceleration is primarily driven by entrainment of background material in the shear layer (P18,M19), not affected by buoyancy in the stream.

We ran the $(M_{\rm b},\delta_{\rm c},\mu)=(1.0,100,0.1)$ simulation with a factor two higher spatial resolution throughout the simulation domain. The results of this simulation are shown in Figs. 8 and 10. The evolution of $h_{\rm b}$ and stream velocity, $V$, are nearly indistinguishable from our fiducial resolution. The penetration of the shear layer into the stream is slightly enhanced, with $h_{\rm s}\sim 10\%$ larger in the high-resolution run at $t\sim 8t_{\rm sc}$. This is still significantly less than the no-gravity simulation, supporting our conclusion that self-gravity significantly suppresses shear layer growth inside the stream.

4.3 Stream Fragmentation due to GI

We here examine the evolution of streams undergoing GI in our simulations, and in particular the properties of clumps formed within them. Regardless of whether GI is dominated by surface or body modes in the linear regime, the end result is always expected to be the collapse of dense, long-lived clumps along the stream axis (N87, H98, Heigl et al., 2016, 2018b).

Figure 11 shows the evolution of three simulations, each with $(M_{\rm b},\delta_{\rm c})=(1.0,100)$. The left-hand column shows the no-gravity simulation with $\mu=0.9$, while the centre and right-hand columns show the gravity simulations with $\mu=0.9$ and $0.4$, respectively. By $t=2t_{\rm sc}$, the non-gravitating stream has developed a well defined shear layer which has penetrated into both the background and the stream, inducing a turbulent mixing zone and diluting the stream density. Meanwhile, the interior of the stream shows numerous density fluctuations caused by turbulence and shocks, with overdensities of up to \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}1.5 times the unperturbed density. At later times the shear layer continues to grow, reaching $h_{\rm s}\sim 0.4R_{\rm s}$ at $t\sim 4t_{\rm sc}$, when the fraction of unmixed fluid in the stream, with $\psi>0.96$, is $\sim 50\%$. This is very similar to the no-gravity simulation shown in the left-hand column of Fig. 7, and is consistent with the evolution of KHI in a constant density stream with $\delta=100$ (M19, figure B1), showing that the steep density profile associated with $\mu=0.9$ does not qualitatively alter the evolution.

Comparing to the corresponding self-gravitating stream, we see that the initial KHI has been suppressed by the introduction of gravity. At $t=4t_{\rm sc}$, the stream appears relatively unperturbed, with no shear layer and only minor density perturbations. The fraction of unmixed fluid in the stream is $77\%$. By $t\sim 6t_{\rm sc}$, small density perturbations can be seen along the stream axis, with a wavelength of $\sim 6.5R_{\rm s}$, slightly larger than the shortest unstable wavelength for GI predicted by H98777While the fastest growing mode in this case is $\lambda_{\rm max}\sim 11R_{\rm s}$, corresponding to 3 clumps, the growth rate at $\sim 6.5R_{\rm s}$ is only $0.85$ times the growth rate at $\lambda_{\rm max}$, and the resulting power spectrum is roughly flat in the range $\lambda\sim(6.5-12)R_{\rm s}$. (Table 1). These density peaks are associated with an axisymmetric distortion of the stream-background interface, despite the fact the the initial perturbations had an equal mix of axisymmetric ($m=0$) and helical ($m=1$) modes. As described in §2.1 and §2.2, GI is unstable only to $m=0$ modes, while the dominant KHI mode at late times has either a long-wavelength and $m=1$ or a short wavelength and $m>1$. This supports the fact that these density perturbations were not amplified by nonlinear KHI, but rather by GI. By $t\sim 8t_{\rm sc}$, these density perturbations have evolved into five dense clumps along the box length of $32R_{\rm s}$, two of which merge by $t\sim 10t_{\rm sc}$.

The evolution of the lower line-mass stream, with $\mu=0.4$, is different. Despite being in the regime where GI dominates over KHI (Fig. 4), at early times the evolution appears dominated by KHI. By $t\sim 4t_{\rm sc}$, a shear layer has developed around the stream, turbulent density fluctuations are visible, and the fraction of unmixed fluid in the stream is $65\%$. This is because the ratio $t_{\rm max}/t_{\rm shear}$ is larger and closer to 1, allowing KHI to develop further before GI takes over. However, by $t\sim 6t_{\rm sc}$, GI has begun to dominate, developing an axisymmetric pattern in the stream-background interface associated with density perturbations along the stream axis, characteristic of GI but not of nonlinear KHI. By $t\sim 8t_{\rm sc}$, five proto-clumps are visible along the stream axis, consistent with the predicted $\lambda_{\rm cr}$. Two of these clumps merge by $t\sim 10t_{\rm sc}$, leaving four large clumps. Assymptotically, for both $\mu=0.4$ and $0.9$, the spacing between clumps is predicted to be $\lambda_{\rm max}\sim 11R_{\rm s}$, the fastest growing GI mode, corresponding to $3$ clumps across $32R_{\rm s}$.

To study the properties of clumps in the simulations, we first select all cells with stream density greater than the break in the PDF of the corresponding snapshot, $\rho_{\rm s,th}$ (Fig. 6). We then group together neighbouring cells above this threshold, removing groups containing fewer than 30 cells to avoid spurious density fluctuations. Varying $\rho_{\rm s,th}$ by 0.1 dex, or using $\rho$ rather than $\rho_{\rm s}$, does not change the number of identified clumps, changes the clump masses by \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}20\%, and the other clump properties discussed below by \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}10\%.

Figure 12 shows several properties of clumps identified in our simulations as a function of time, where $t=0$ is set to the first timestep where clumps have been identified. We show the clump mass, $M_{\rm c}$, the turbulent Mach number within the clumps, $\mathcal{M}_{\rm turb}=\sigma_{\rm turb}/c_{\rm s}$, and the clump virial parameter, defined as

[TABLE]

where the factor $5/3$ comes from assuming a constant density profile inside the clump. If $\alpha_{\rm vir}\sim 1$, the clump is in virial equilibrium, while $\alpha_{\rm vir}<1$ implies the clump is collapsing and $\alpha_{\rm vir}>1$ implies it is unbound. For each property we display the average over all clumps identified in a given snapshot, typically four to five clumps.

Following the initial collapse when the clump mass grows significantly, it tends to saturate at a well defined value despite some oscillations. These oscillations, on the order of $\sim 10-20\%$, are due in part to our density threshold for clump cells, which is recalibrated at each snapshot. We have normalised the mass in Fig. 12 by $M_{i}=M_{\rm stream}/N_{\rm clump}$, where $N_{\rm clump}$ is the number of clumps in the stream and $M_{\rm stream}=\pi R_{\rm s}^{2}L{\overline{\rho_{\rm s}}}$ is the initial stream mass with ${\overline{\rho_{\rm s}}}$ the mean density in the stream. $M_{\rm i}$ is thus the typical clump mass one would expect if the entire initial stream fragments into clumps. We find that $M_{\rm c}/M_{\rm i}$ increases with $\mu$, rising from $\sim(0.4-0.9)$ for $\mu=(0.3-0.9)$, independent of $M_{\rm b}$ or $\delta_{\rm c}$.

The spherical Jeans mass obtained using the average properties in the initial stream is $M_{\rm J}=(\pi^{5/2}/6){\overline{c_{\rm s}}}^{3}G^{-3/2}{\overline{\rho_{\rm s}}}^{-1/2}$. For $N_{\rm clumps}=4$ and $L=32R_{\rm s}$ we obtain $M_{\rm i}/M_{\rm J}\sim 0.14(t_{\rm sc}/t_{\rm ff})^{3}$, with $t_{\rm sc}\simeq 2R_{\rm s}/{\overline{c_{\rm s}}}$ and $t_{\rm ff}=(4G{\overline{\rho_{\rm s}}})^{-1/2}$. This corresponds to $M_{\rm i}/M_{\rm J}\sim(0.2-1.2)$ for $\mu=(0.3-0.9)$ (Table 1), yielding clump masses $M_{\rm c}\sim(0.1-1)M_{\rm J}$. For small $\mu$, when the density profile in the initial stream is roughly constant, the Bonnor-Ebert mass (eq. 7) is $M_{\rm BE}\sim 0.5M_{\rm J}$. In general, $M_{\rm BE}>M_{\rm c}$ for $\mu<1$.

The turbulent Mach number increases by a factor of $\sim 3$ as $\mu$ is decreased from 0.9 to 0.3. However, in all cases $\mathcal{M}_{\rm turb}$ is \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}0.3 asymptotically, and does not exceed $\sim 0.6$ during the initial collapse of the clump. Turbulent support is thus negligible compared to thermal pressure. The clump virial parameter increases from $\alpha_{\rm vir}\sim 1$ for $\mu=0.9$, consistent with $M_{\rm c}/M_{\rm J}\sim 1$ in this case, to $\alpha_{\rm vir}\sim 2.3$ for $\mu=0.3$.

The additional support for clumps in simulations with lower values of $\mu$ comes from the external pressure, which also played a larger role in confining the initial stream. This can be seen by considering the full virial parameter including the surface pressure term (e.g. Krumholz, 2015, Chapter 6). We approximate this as

[TABLE]

This is shown by dot-dashed lines in the rightmost panel of Fig. 12. For $\mu=0.9$, the external pressure is negligible and the two virial parameters are nearly identical. However, for $\mu=0.3$, ${\tilde{\alpha}}_{\rm vir}\sim 1.4$, indicating that the clumps in this case are primarily confined by external pressure. While this is still larger than 1, eq. (44) is only an approximation, assuming a spherical clump with constant density and uniform external pressure. Properly accounting for the density profile within the clump tends to reduce the virial parameter compared to eqs. (43)-(44) (e.g. Mandelker et al., 2017). Given this, a value of ${\tilde{\alpha}}_{\rm vir}\sim 1.4$ is indicative of the clumps being in approximate virial equilibrium due to a combination of gravitational and pressure confinement.

Contrary to the strong dependence of clump properties on $\mu$, their dependence on $(M_{\rm b},\delta_{\rm c})$ at fixed $\mu$ is extremely weak. $M_{\rm c}$ and $\alpha_{\rm vir}$ vary by only a few percent as $\delta_{\rm c}$ varies from $6.7-100$ or $M_{\rm b}$ from $1-2.5$. Furthermore, clumps formed in simulations of pure GI, with $(M_{\rm b},\delta_{\rm c})=(0,100)$ (see the Appendix §A), have masses only $\sim 10\%$ larger than those in simulations with $M_{\rm b}=1$ for both $\mu=0.9$ and $0.4$. We conclude that once GI dominates over KHI and leads to clump formation, KHI has little effect on the resulting clump properties even if $\mu$ only slightly exceeds $\mu_{\rm cr}$.

To check convergence, we repeated the $(M_{\rm b},\delta_{\rm c},\mu)=(1.0,100,0.9)$ simulation with a factor two higher spatial resolution, and show the results in Fig. 12. No significant change was found in the number of clumps, their formation time, or their properties. The clump mass increases by \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}4\%, while $\mathcal{M}_{\rm turb}$ and $\alpha_{\rm vir}$ are unchanged. We conclude that our fiducial resolution is sufficient to resolve the stream fragmentation and resulting clumps.

In summary, clumps forming in higher line-mass streams are more massive, have lower turbulent Mach numbers and lower virial parameters. This is primarily due to the larger degree of external pressure support for low line-mass streams present in the initial conditions, with a small contribution from enhanced mixing and dilution in lower line-mass streams caused by more efficient KHI. At fixed $\mu$, the variation of clump properties with $M_{\rm b}$ and $\delta_{\rm c}$ is very small. For $\mu=0.9$, roughly $90\%$ of the initial stream mass winds up in clumps, which following collapse are in approximate virial equilibrium at the thermal Jeans scale. For $\mu=0.3$, only \lower 2.58333pt\hbox{;\buildrel<\over{\sim};}50\% of the initial stream mass is in clumps, the rest having mixed into the background due to KHI. The collapsed clumps have $M_{\rm c}\sim(0.1-0.2)M_{\rm J}$, and are confined by external pressure. In all cases, the turbulent pressure in the collapsed clumps is negligible, with turbulent Mach numbers $\sim(0.1-0.3)$ for $\mu=(0.9-0.3)$.

5 discussion

5.1 Astrophysical Applications

Our results on the combined evolution of KHI and GI in self-gravitating filaments have several astrophysical implications. In this section, we highlight potential applications for studies of star-forming filaments in the ISM and for cold streams feeding massive galaxies at high redshift.

5.1.1 High-z Intergalactic Streams

Massive galaxies with baryonic masses \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}10^{11}{\rm M}_{\odot} at $z\sim(1-4)$ reside in halos with virial masses M_{\rm vir}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}10^{12}{\rm M}_{\odot}. The CGM of these galaxies is thought to contain hot gas with T\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}10^{6}\,{\rm K} in approximate hydrostatic equilibrium. However, the star-formation rates measured in these galaxies of \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}100\,M_{\odot}\,{\rm yr}^{-1} is significantly larger than expected from the cooling of the hot CGM, and their prevalence exceeds that expected from mergers (Dekel et al., 2009a). As outlined in §1, such galaxies are fed by cold, $T\sim 10^{4}\,{\rm K}$ gas streams from the cosmic web, which efficiently penetrate the hot halo all the way to the central galaxy (Kereš et al., 2005; Dekel & Birnboim, 2006; Ocvirk et al., 2008; Dekel et al., 2009a; Ceverino et al., 2010; Faucher-Giguère et al., 2011; van de Voort et al., 2011). The shearing against the hot CGM makes these streams susceptible to KHI. This has motivated several detailed studies of KHI in such systems, with $\delta\sim(30-100)$ and $M_{\rm b}\sim(0.5-2)$ (M16; P18; M19). As cosmological simulations lack the spatial resolution to properly resolve KHI in the streams, these studies have been idealized, accounting thus far only for non-radiative hydrodynamics without gravity.

These studies find that sufficiently narrow streams, with R_{\rm s}/R_{\rm v}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}(0.005-0.05) where $R_{\rm v}$ is the halo virial radius, will disrupt in the CGM before reaching the central galaxy. The threshold value of $R_{\rm s}$ depends on $(M_{\rm b},\delta)$. However, our results suggest that in a certain regime of parameter space, self-gravity may stabilize streams and halt their disruption. Even if the line mass is very low compared to the critical value, $\mu\sim 0.1$, we find that buoyancy can prevent the shear layer from penetrating the inner stream (Figs. 7 and 8). For $\delta=100$, we find that the penetration rate of the shear layer into the stream is reduced by a factor of \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}3 when self-gravity is included (Fig. 8). This implies that the previous estimates of the upper limit on the radius of streams that can disrupt in the CGM should be reduced by a similar factor, namely R_{\rm s}/R_{\rm v}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}(0.0015-0.015). Very narrow streams may thus survive the journey to the central galaxy, though they are likely to reach it somewhat wider and more diluted than they began.

M19 also found that typical streams can significantly decelerate in the CGM, dissipating $\sim(10-50)\%$ of their bulk kinetic energy before the central galaxy. If this energy is subsequently radiated away, it can significantly contribute to the Ly$\alpha$ emission observed in the CGM of massive high-$z$ galaxies. Our results show that the self-gravity of the gas is unlikely to alter this conclusion, because the deceleration rates and the entrainment of background mass are unaffected (Figs. 8 and 10).

Other studies have suggested that at higher redshift, z\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}5, the streams feeding massive galaxies may be gravitationally unstable, with $\mu\sim 1$ (Mandelker et al., 2018). These authors speculated that such streams could gravitationally fragment while still in the halo, and that this could lead to the formation of metal-poor globular clusters and stars directly in the halos of high-$z$ galaxies. While this study did not account for KHI in the streams, our results suggest that this is unlikely to affect their conclusions, since for $\mu\sim 1$ and $M_{\rm b}\sim 1$, GI is unaffected by KHI (Figs. 11 and 12). We note that for cosmic web filaments far from haloes, only GI operates as no shear is expected.

5.1.2 ISM Filaments

As outlined in §1, numerous filametary structures are observed in the ISM, in particular in star-forming regions such as giant molecular clouds. While much attention has been payed to the gravitational stability and fragmentation of such filaments, these studies do not consider KHI induced by shearing motions between the filament and its surroundings. This is despite the fact that strong shearing motions and even signatures of KHI have been detected in molecular clouds and around filaments (e.g. Rodriguez-Franco et al., 1992; Berné et al., 2010; Berné & Matsumoto, 2012). Numerical simulations of molecular clouds in the central molecular zone have also revealed strong shearing motions which generate turbulence and reduce the SFR by a factor of $\sim 7$ compared to nearby clouds (Federrath et al., 2016). It is thus important to consider how KHI might affect the fragmentation of ISM filaments.

We note that in this case, the shearing motion is thought to be due to a background “wind” flowing across a roughly static filament, rather than a stream flowing through a static background. However, due to Galilean invariance, these two scenarios should behave identically.

The regions surrounding ISM filaments are often extremely turbulent, with turbulent Mach numbers of order 10 or higher, and the filaments themselves are often supervirial, with \mu\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}1. This is obviously very different from our initial conditions of a smooth filament in hydrostatic equilibrium (see §5.2). However, subvirial filaments with $\mu<1$ have been observed (Henshaw et al., 2016; Hacar et al., 2018, and references therein; Orkisz et al., 2019). In some cases, these low line-mass filaments host pre-stellar cores which are at least partly supported by external pressure (Kirk et al., 2017; Seo et al., 2018), consistent with GI in filaments with $\mu_{\rm cr}<\mu<1$ (Fig. 12). If $\mu$ is known, this will constrain $\mu_{\rm cr}$, which in turn can be used to place constraints on the properties of the confining medium, and in particular on $M_{\rm b}$, the velocity of the shear flow between the filament and the background (Fig. 5).

5.1.3 Tidal Disruption Events

Stars that wander too close to a supermassive black hole, such as found in the centres of most massive galaxies, can be disrupted by the strong tidal forces exerted by the black hole (Rees, 1988). Following the disruption, the stellar debris often evolves into a gas stream which partly accretes onto the black hole producing a luminous flare. Following their formation, the tidal shear of the black hole renders these streams gravitationally stable, so bound clumps are unlikely to form along the stream. Furthermore, the streams can be treated as approximately in hydrostatic equilibrium in the cylindrically-radial direction (Coughlin & Nixon, 2015; Coughlin et al., 2016b, a). Recently, it has been argued that interactions between the debris stream and the ambient tenuous gas near the galactic centre can render such streams unstable to KHI, with nominal disruption times shorter than the infall timescale of the stream onto the black hole (Bonnerot et al., 2016). If true, this would significantly reduce the expected luminosity of the accretion flare. However, as we have shown, even in weakly self-gravitating streams, total stream disruption is significantly delayed due to buoyancy within the stream. This would mean that KHI in the streams below $\mu_{\rm cr}$ will be stopped by buoyancy, and the decrease in the flare-luminosity predicted by Bonnerot et al. (2016) may be overestimated. Such a scenario can be tested with dedicated simulations.

5.2 Caveats and Additional Physical Effects

While our analysis has focused on elucidating the interplay between KHI and GI in filaments, applications of our results to astrophysical scenarios require careful consideration of additional physical processes that have not yet been taken into account. These include the assumed isentropic initial conditions and lack of radiative cooling, the assumption of line mass ratios $\mu<1$ and hydrostatic equilibrium in the initial conditions, the lack of magnetic fields, and (in the case of cold streams feeding massive galaxies at high redshift) the lack of a dark matter component to the gravitational potential. In this section, we speculate as to the possible effects of these processes, all of which will be explored in future work.

Radiative cooling is clearly very important for both ISM filaments and intergalactic gas streams. Both of these are expected to have cooling times much shorter than their sound crossing times, which is why they are often modeled as isothermal. Radiative cooling can either enhance or suppress KHI in the linear regime, depending on the slope of the cooling function and on the ratio of the cooling time in each fluid to the sound crossing time (Massaglia et al., 1992; Bodo et al., 1993; Vietri et al., 1997; Hardee & Stone, 1997; Xu et al., 2000). However, when these ratios are either much larger or much smaller than unity, the linear growth rates are similar to the adiabatic case at longitudinal wavelengths \lambda\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}R_{\rm s} (Mandelker et al., in prep.). Even in this case, cooling can substantially alter the nonlinear evolution of KHI (Vietri et al., 1997; Stone et al., 1997; Xu et al., 2000; Micono et al., 2000), though the net effect again depends on details of the cooling function and the stream parameters. Some authors have found that cooling leads to more violent disruption of the stream (Stone et al., 1997; Xu et al., 2000), while others have found that it prevents stream disruption by limiting the penetration of the shear layer into the stream (Vietri et al., 1997; Micono et al., 2000). If shear layer growth is suppressed and the contact discontinuity maintained, then $t_{\rm shear}$ will increase and $\mu_{\rm cr}$ will decrease (Fig. 4). Thus, the regime where GI dominates over KHI will expand. Furthermore, it is also found that KHI in a cooling medium leads to much larger density fluctuations, and to the formation of dense knots and filaments inside the stream. These are likely to further enhance GI and filament fragmentation. Cooling is also likely to allow the clumps to collapse to higher densities and reach lower temperatures, thus decreasing their Jeans mass and leading to further fragmentation and collapse.

Magnetic fields are likely to be dynamically important in ISM filaments. This can have a stabilizing effect on GI, especially when $\mu<1$ (e.g. N87, H98), and also on KHI, where magnetic fields parallel to the flow have been found to stabilise high-$m$ modes and suppress shear layer growth (Ferrari et al., 1981; Birkinshaw, 1990). It is therefore unclear what the net effect will be in terms of the competition between these two processes, and this will likely depend sensitively on the properties of the field. For intergalactic gas streams at high redshift, magnetic fields are likely dynamically unimportant (e.g. Bagchi et al., 2002). Nevertheless, they may significantly weaken thermal conductivity and viscosity, which will influence the width of the shear layer (M19) and thus affect the instability. All these effects should be accounted for simultaneously in future work.

When considering intergalactic gas streams, we must also account for the contribution of the host dark matter filament to the gravitational potential. To our knowledge, the gravitational stability of a gas stream embedded in a dark matter filament has not been studied. The dark matter may stabilise the stream by making it more buoyant, or it may destabilise the stream by increasing the inward radial gravitational force, thus requiring non-thermal turbulent motions to support the stream against radial collapse. This may also suppress KHI by further limiting shear layer growth and stream disruption (see Figs. 7-8). The central dark matter halo into which the streams are flowing will also affect their evolution. The central potential focuses the stream into a conical shape with its radius decreasing towards the halo centre, $R_{\rm s}\propto r$. (Dekel et al., 2009a; van de Voort & Schaye, 2012). This decreases the KHI timescales, which are proportional to $R_{\rm s}$ (eqs. 16-17). However, this focusing also increases the stream density, with $\rho\propto R_{\rm s}^{-2}\propto r^{-2}$, resulting in a decrease of the free-fall time, $t_{\rm ff}\propto\rho^{-1/2}\propto r$. Since $t_{\rm max}\propto t_{\rm ff}$, the ratio $t_{\rm max}/t_{\rm shear}$ is unlikely to vary significantly throughout the halo, as is the critical line-mass ratio, $\mu_{\rm cr}$. However, this must be studied in more detail, as must the effect of gravitational acceleration towards the halo centre on the evolution of KHI and GI in intergalactic cold streams.

Throughout our analysis, we assumed that filaments began in hydrostatic equilibrium, and without any internal non-thermal support such as turbulence or vorticity. This is unlikely to be the case for either ISM filaments or intergalactic streams. Theoretical studies of GI in ISM filaments growing self-consistently via radial accretion have shown that turbulence builds up inside the stream with Mach numbers of order unity and contributes to its support (Heitsch, 2013; Clarke et al., 2016; Clarke et al., 2017; Heigl et al., 2018a). Despite this, the filament was found to fragment when its line mass reached the critical value for hydrostatic equilibrium, namely at \mu\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}1, in a similar manner to the $\mu<1$ filaments considered here, leading to the formation of Jeans-scale clumps (Clarke et al., 2016; Clarke et al., 2017). It is unclear how these results will change in the presence of KHI. Likewise, it has been suggested that accretion onto cosmic gas streams from the intergalactic medium creates specific profiles (Fillmore & Goldreich, 1984; Birnboim et al., 2016), induces roughly sonic turbulence (Mandelker et al., 2018) and vorticity (Codis et al., 2012; Codis et al., 2015; Laigle et al., 2015), and grows streams to $\mu>1$ (Mandelker et al., 2018). Such non-equilibrium effects must be considered in order to describe stream evolution.

6 Summary and Conclusions

Self-gravitating gaseous filaments are ubiquitous in astrophysics, from sub-pc filaments within the interstellar medium, to Mpc scale streams feeding galaxies along the cosmic web. As such, they may be subject to gravitational instability (GI), which leads to stream fragmentation and to the formation of long-lived, collapsed clumps along the stream axis. In many cases, such filaments are also susceptible to Kelvin-Helmholtz Instability (KHI) due to a shear flow against a confining background medium, which acts to mix the filament with the background fluid via a turbulent shear layer. Motivated by this, we have performed the first ever study of the evolution of a self-gravitating filament or stream undergoing KHI, using simple analytic models and hydrodynamic simulations. Such a system is characterised by three dimensionless parameters: the Mach number of the stream with respect to the sound speed in the (static) background, $M_{\rm b}$, the ratio of the central density in the stream to the background density outside the stream, $\delta_{\rm c}$, and the ratio of the mass-per-unit-length (line-mass) of the stream to the maximal line-mass for which initial hydrostatic equilibrium is possible, $\mu$. The current analysis is restricted to filaments with $\mu<1$ initially in hydrostatic equilibrium. Our main results can be summarised as follows:

The competition between GI and KHI is governed by the ratio of the timescale for linear growth of the fastest growing GI mode, $t_{\rm max}$, and the relevant nonlinear KHI timescale. When GI is dominated by surface modes, this is the time for the KHI-induced shear layer to expand to a size comparable to the stream radius and destroy the initial contact discontinuity, $t_{\rm shear}$. If $t_{\rm max}/t_{\rm shear}<1$, GI causes the stream to fragment into long-lived clumps and suppresses mixing with the background medium. Likewise, if $t_{\rm max}/t_{\rm shear}>1$, KHI mixes the stream with the background medium, dilutes its density and suppresses clump formation (Fig. 4). Regardless, the stream is always unstable. When GI is dominated by body modes, clumps may form even when $t_{\rm max}$ is slightly longer than $t_{\rm shear}$, since the contact discontinuity no longer plays a role in GI. 2. 2.

The timescale criterion can be rephrased as a criterion on the line-mass ratio $\mu$. If this is smaller than a critical value which depends on the Mach number and density contrast, $\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$, then KHI will win and mix the stream and background. However, if $\mu>\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$, the stream will fragment into long-lived, bound clumps. $\mu_{\rm cr}$ increases strongly with $M_{\rm b}$, and has a weak tendency to increase with $\delta_{\rm c}$ (Fig. 5). For M_{\rm b}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}2.5 we have \mu_{\rm cr}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}0.5. At larger Mach numbers, when KHI is dominated by high-order azimuthal surface modes, \mu_{\rm cr}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}0.9. In practice, values of $\mu_{\rm cr}>0.9$ are not relevant, as GI body modes will cause clump formation even if \mu\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}\mu_{\rm cr}. 3. 3.

When $\mu<\mu_{\rm cr}$, the evolution of KHI outside the stream boundary is similar to the case of a non-gravitating uniform density stream, studied in detail by Mandelker et al. (2019). Self-gravity slows the expansion of the shear layer into the background by less than $20\%$ for large $\mu$, and significantly less than that for smaller $\mu$ (Fig. 8). Consequently, the stream deceleration due to entrainment of background mass in the shear layer is also unaffected, and follows the analytical prediction (Fig. 10). 4. 4.

However, gravity does qualitatively affect the penetration of the shear layer into the stream. At t\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}(2-3)t_{\rm sc}, the penetration rate of the shear layer into the stream is slowed by a factor of \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}3 compared to the no-gravity case (Fig. 8). This is due to restoring buoyancy forces in the stream interior, corresponding to values of the Richardson number, ${\rm Ri}>0.25$ (Fig. 9). This significantly slows the final disruption of the stream by KHI, as a dense central core remains partly shielded against mixing (Fig. 7). 5. 5.

The clumps that form by GI when $\mu>\mu_{\rm cr}$ are largely unaffected by KHI. They are typically less massive than the Jeans mass, and supported partially by external pressure. However, as $\mu\rightarrow 1$ the clumps approach the Jeans mass and the external pressure support becomes negligible. In all cases, the internal turbulent motions are subsonic and turbulent pressure support is negligible, though the turbulent Mach number increases towards lower $\mu$ (Fig. 12). KHI seems to have a minor effect on the clump properties, which are largely insensitive to the Mach number of the flow, even in the static limit, $M_{\rm b}=0$. 6. 6.

Our finding that self-gravity may shield the inner core of filaments from disruption by KHI, implies that recent studies of KHI in gas streams feeding massive galaxies at high-$z$ may have overestimated the disruption of these streams in the CGM. However, the dissipation and deceleration rates should not be affected. Additionally, our finding that GI induced fragmentation only occurs when $\mu>\mu_{\rm cr}(M_{\rm b},\delta_{\rm c})$ can be used to place constraints on the properties and kinematics of the confining medium surrounding low mass filaments in the ISM. However, in order to properly address these phenomena, additional physics such as radiative cooling, magnetic fields, external gravitational potential, and non-thermal turbulent motions, will have to be added to our models.

Acknowledgments

We thank Romain Teyssier for many helpful suggestions while running the simulations. We thank Frank van den Bosch, Frederic Bournaud, Andreas Burkert, Drummond Fielding, Shuo Kong, Diederik Kruijssen, and Xun Shi for insightful discussions. NM acknowledges support from the Klaus Tschira Foundation through the HITS Yale Program in Astrophysics (HYPA). The simulations were performed on the Omega and Grace HPC clusters at Yale. This work is supported in part by the facilities and staff of the Yale Center for Research Computing. AD was partly supported by the grants BSF 2014-273, GIF I-1341-303.7/2016 and NSF AST-1405962.

Appendix A GI Growth Rates and Surface Vs Body Modes

For an incompressible pressure-confined cylinder, the linear growth-time of the fastest growing GI mode is given by eq. (11), based on H98. In our case, however, the stream is highly compressible and its density can be far from constant (Fig. 3), so it is unclear whether eq. (11) remains valid. On the other hand, for isothermal cylinders confined by a zero density background, the ratio $t_{\rm max}/t_{\rm ff}$ does not vary much with line-mass (N87). If the same is true for a non-isothermal cylinder confined by an arbitrary density background, then eq. (11) may also apply to our case. To test this, we performed simulations with self-gravity but without shear flow, $M_{b}=0$, and with a single perturbation wavelength, $\lambda=1/4=8R_{\rm s}$. We examined two different values of the line-mass, $\mu=0.4$ and $0.9$, with $\delta_{\rm c}=100$. Based on the analysis of N87, the former is expected to be unstable to surface modes, while the latter to body modes.

Fig. 13 shows the perturbation amplitude as a function of time, measured in each simulation at both the surface of the stream and along its axis. The former is defined by the deformation of the stream-background interface,

[TABLE]

where $r$ specifies the radius at which the tracer variable $\psi=0.5$. For the perturbation along the stream axis we use the density

[TABLE]

As can be seen in the top panels of Fig. 13, these two measurements of the perturbation amplitude yield similar growth rates. At t\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}7t_{\rm ff}, the perturbation amplitude is roughly constant. This is approximately four times the stream sound crossing time, $t_{\rm sc}\propto 2R_{\rm s}/c_{\rm s}$ (see eq. 41 and Table 1), which is the sound crossing time of the perturbation wavelength, $\lambda=8R_{\rm s}$. This is the coherence time of the perturbation, during which the initial velocity perturbation is converted into a growing eigenmode of the system (see M16 and M19 for a discussion of a similar phenomenon in KHI). At later times the perturbation amplitude is well fit by $\Delta\propto\exp{(\omega t)}$. The best fit growth rates are $\omega_{\rm fit}\approx 0.96\omega_{\rm Hunter}$ for $\mu=0.4$, and $\omega_{\rm fit}\approx 0.90\omega_{\rm Hunter}$ for $\mu=0.9$, where the profile deviates further from the constant density assumption in H98. We conclude that the growth time for GI can be approximated by the H98 dispersion relation even for large values of the line-mass, $\mu$.

The bottom panels compare the 1-D power spectra of density perturbations along the stream axis, at $r=0$, and near its edge, at r\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}R_{\rm s}. We show here the power measured at the wavelength $8R_{s}$, but note that the results are nearly identical when showing the total power, as power at all other scales is small. For $\mu=0.4$, the power near the stream edge is larger, while the opposite is true for $\mu=0.9$, consistent with these two simulations corresponding to surface and body modes, respectively.