The impact of magnetic fields on cold streams feeding galaxies

TL;DR

This study investigates how magnetic fields influence the stability and longevity of cold gas streams feeding high-redshift galaxies, revealing that magnetic fields can significantly extend stream survival despite enhancing instabilities.

Contribution

The paper provides the first spatially resolved idealized model showing magnetic fields can prolong cold stream survival by suppressing mixing, challenging previous assumptions about their destabilizing effects.

Findings

Magnetic fields do not inhibit but enhance Kelvin-Helmholtz instability growth.

Magnetic fields can increase stream survival times by 2-8 times.

Thinner streams with magnetic fields can reach galaxies, aiding sustained star formation.

Abstract

High redshift, massive halos are observed to have sustained, high star formation rates, which require that the amount of cold gas in the halo is continuously replenished. The cooling time scale for the hot virialized halo gas is too long to provide the source of cold gas. Supersonic, cold streams have been invoked as a mechanism for feeding massive halos at high redshift and deliver the cold gas required for continued star formation at the rates observed. This mechanism for replenishing the cold gas reservoir is motivated by some cosmological simulations. However, the cold streams are likely to be subject to the supersonic version of the Kelvin-Helmholtz instability (KHI), which eventually leads to stream disruption. Cosmological simulations have yet to obtain the spatial resolution required for understanding the detailed stability properties of cold streams. In this paper, we consider instead an idealized model of magnetized cold streams that we spatially resolve. Using linear theory we show how magnetic fields with dynamically important field strengths do not inhibit the KHI but rather enhance its growth rate. We perform nonlinear simulations of magnetized stream disruption and find that magnetic fields can nevertheless increase stream survival times by suppressing the mixing rate of cold gas with the circumgalactic medium. We find that magnetic fields can allow streams to survive $\sim 2-8$ times longer and, consequently, that streams $\sim 2-8$ times thinner can reach the central galaxy if the magnetic field strength is $\sim 0.3-0.8 \mu$G.

Tables

Table 1. Table 1: Fastest growing eigenmodes for the KHI in a cylindrical cold stream with diameter L 𝐿 L , density ρ s / ρ 0 = 50 subscript 𝜌 s subscript 𝜌 0 50 \rho_{\mathrm{s}}/\rho_{0}=50 , Mach number ℳ 0 = 2 subscript ℳ 0 2 \mathcal{M}_{0}=2 , smoothing length a / L = 0.05 𝑎 𝐿 0.05 a/L=0.05 and β = 1 𝛽 1 \beta=1 .

$m$	$k_{\max }L$	$\sigma L/c_0$
0	3.43372	0.60051
1	4.01473	0.61800
2	2.59792	0.81842
3	3.11557	0.96195
4	3.48196	0.99990
5	3.73777	0.97093

Equations

$${\text{varv}}_{\mathrm{a}}=\frac{B}{\sqrt{{\mu }_0\rho }},$$

$${\mathcal{M}}_0=\frac{V}{c_{\mathrm{ad},0}}.$$

$$\beta =\frac{p}{p_{\mathrm{mag}}}=2\frac{c^2}{{\text{varv}}_{\mathrm{a}}^2},$$

$$\text{varv}(z)$$

$$\frac{\rho (z)}{{\rho }_0}$$

$$\text{varv}(r)$$

$$\rho (r)$$

$$E_{\mathrm{kin},\parallel }=L_x{\int }_{-\infty }^{\infty }\frac{\rho (z)}{2}\text{varv}(z{)}^{2}dz$$

$$E_{\mathrm{kin},\parallel }=2\pi L_z{\int }_0^{\infty }\frac{\rho (r)}{2}\text{varv}(r{)}^{2}rdr$$

$$V{t}_{1/2}=\xi M\sqrt{\gamma }L,$$

$$\frac{L}{R_{\mathrm{v}}}\lesssim\frac{1}{\xi M\sqrt{\gamma}}=\left\{\begin{array}[]{lll}2\times 10^{-2}&\mathrm{for}&\mathrm{hydro}\\ 9\times 10^{-3}&\mathrm{for}&\beta=10\\ 2\times 10^{-4}&\mathrm{for}&\beta=1\\ \end{array}\right.\ ,$$

$$B=\sqrt{\frac{2p{\mu }_0}{\beta }}=0.8{\left(\frac{n_{\mathrm{H}}}{10^{-4}{\mathrm{cm}}^{-3}}\right)}^{1/2}{\left(\frac{T}{10^6\mathrm{K}}\right)}^{1/2}{\beta }^{-1/2}\mu \mathrm{G}.$$

$$\frac{\partial \rho }{\partial t}$$

$$\rho \frac{\mathrm{d}\text{varv}}{\mathrm{d}t}$$

$$\frac{\partial \mathbf{B}}{\partial t}$$

$$\frac{p}{\gamma -1}\frac{\mathrm{d}\ln \left(p{\rho }^{-\gamma }\right)}{\mathrm{d}t}$$

$$-\mathrm{i}\left(\omega -k\text{varv}\right)\frac{\delta \rho }{\rho }$$

$$-\mathrm{i}\left(\omega -k\text{varv}\right)\delta {\text{varv}}_x$$

$$-\mathrm{i}\left(\omega -k\text{varv}\right)\delta {\text{varv}}_z$$

$$-\mathrm{i}\left(\omega -k\text{varv}\right)\frac{\delta A}{B}$$

$$-\mathrm{i}\left(\omega -k\text{varv}\right)\frac{\delta T}{T}$$

$$-\mathrm{i}(\omega -k\text{varv})\delta {\text{varv}}_r$$

$$-\mathrm{i}(\omega -k\text{varv})\delta {\text{varv}}_{\varphi }$$

$$-\mathrm{i}(\omega -k\text{varv})\delta v_z$$

$$-\mathrm{i}(\omega -k\text{varv})\frac{\delta p}{p}$$

$$-\mathrm{i}(\omega -k\text{varv})\delta b_r$$

$$-\mathrm{i}(\omega -k\text{varv})\delta b_{\varphi }$$

$$-\mathrm{i}(\omega -k\text{varv})\delta b_z$$

$${\mathcal{D}}_r\equiv \frac{\partial }{\partial r}+\frac{1}{r},$$

$$\nabla \cdot \delta \mathbf{b}={\mathcal{D}}_r\delta b_r+\frac{\mathrm{i}m}{r}\delta b_{\varphi }+\mathrm{i}k\delta b_z=0,$$

Full text

††pubyear: 2019††pagerange: The impact of magnetic fields on cold streams feeding galaxies–A.4

1 Introduction

High redshift ($z\sim 2-3$), massive galaxies residing in $\sim 10^{12}M_{\odot}$ dark matter halos are observed to have high star formation rates of $\sim 100M_{\odot}\mathrm{yr}^{-1}$ (Genzel et al., 2006; Förster Schreiber et al., 2006). This peak epoque of cosmic star formation requires a continuous source of new cold gas. The galactic cold gas reservoir, which is needed to fuel intense star formation, would otherwise be depleted on a time scale shorter than observed. Determining whether the star formation is sustained by recycling of gas, in situ cooling or whether it is provided by an external source is an important topic in galaxy formation (see e.g. Somerville & Davé 2015, Naab & Ostriker 2017 and Dayal & Ferrara 2018 for reviews of galaxy formation).

Theories of structure formation predict that the halo gas is heated to the virial temperature, $T_{\mathrm{v}}$, by an accretion shock when the halo mass is $\gtrsim 10^{12}M_{\odot}$ (Birnboim & Dekel, 2003). Such hot ($T_{\mathrm{v}}\gtrsim 10^{6}$ K) virialized gas has a prohibitively long cooling time for making in situ cooling a plausible source of cold gas (Rees & Ostriker, 1977; White & Rees, 1978; Silk, 1977). An alternative source of cold gas could be mergers but observations have indicated that star formation takes place in a large rotating disk, disfavoring this scenario as well (Genzel et al., 2006; Stark et al., 2008).

In the cold stream model of galaxy formation in massive halos (Dekel & Birnboim, 2006; Dekel et al., 2009), the central galaxy in a massive hot halo is instead fed with cold gas that streams along filaments in the cosmic web and penetrates down into the central galaxy. Such cold streams have been observed in some cosmological simulations (e.g. Kereš et al. 2005; Ocvirk et al. 2008; Van De Voort et al. 2012; Goerdt & Ceverino 2015), where they are able to provide the halo with $\sim 100M_{\odot}\mathrm{yr}^{-1}$ cold gas, in rough agreement with the observed star formation rate.

State-of-the-art zoom-in cosmological simulations are however currently in disagreement on whether the cold stream model outlined above works. While some simulations do find cold streams that feed the central galaxy, others find that the streams are disrupted and heated inside $\lesssim 0.25-0.5R_{\mathrm{v}}$ (where $R_{\mathrm{v}}$ is the virial radius, Nelson et al. 2013). The discrepancy between different cosmological simulations has been attributed to numerical differences in simulations employing smoothed particle hydrodynamics (SPH), Eulerian adaptive mesh refinement (AMR), and a quasi-Lagrangian finite volume approach, i.e., a moving mesh (Nelson et al., 2013).

A key question that needs to be addressed for the cold stream model is whether streams can in fact penetrate sufficiently far into the halo or whether they will be disrupted by hydrodynamical instabilities during their propagation. The question cannot be directly answered with current-day zoom-in cosmological simulations, because such simulations have not attained the required spatial resolution to accurately capture potentially disrupting hydrodynamical instabilities. Cold streams have been observed to have widths of 1-10 per cent of the virial radius while the true width could be even smaller since those values are close to the numerical resolution limit. For a halo with $\sim 50$ kpc virial radius, this means that the thinnest streams are expected to have width of $\sim 1$ kpc (Mandelker et al., 2019) or smaller. They are therefore not well resolved in cosmological simulations, despite recent progress where the spatial resolution in the circumgalactic medium (CGM) is enhanced to $\lesssim 1$ kpc (Suresh et al., 2019; Hummels et al., 2018; van de Voort et al., 2019; Peeples et al., 2018). We note that the resolution requirement for studying the multi-phase gas of the CGM might be a general hindrance for obtaining high fidelity cosmological simulations of the CGM in the near future, as some studies show that it is necessary to resolve the gas all the way down to sub-pc scales (McCourt et al., 2018; Sparre et al., 2019).

As an alternative to studying cold flows in full cosmological simulations, Mandelker et al. (2016); Padnos et al. (2018); Mandelker et al. (2019) have pursued the problem of stream disruption using analytical theory and idealized simulations. They found that the main disruption mechanism for a cold, dense stream propagating through a less dense but hotter background, is the Kelvin-Helmholtz instability (KHI). For supersonic propagation, the KHI can manifest itself as sound waves reflecting and growing in amplitude inside of the stream (called reflective or body modes, Payne & Cohn 1985; Hardee & Norman 1988). This type of KHI differs from the familiar textbook version (e.g. Chandrasekhar 1961; Drazin & Reid 2004) and the resulting breakup of the stream does not contain the characteristic cat’s eye vortex (see Figs. 5 and 15 in Padnos et al. 2018 for a comparison between the two types of instability).

Mandelker et al. (2019) found that mixing of the stream with the background CGM (and deceleration of the stream) occurs much faster in three dimensions in comparison to two dimensions. It was concluded that the KHI can lead to disintegration of cold streams before they reach the galaxy residing at the center of a massive halo. For typical parameters, estimates showed that streams with width less than $\sim 1-10$ per cent of the virial radius will be disintegrated (Mandelker et al., 2019). The deceleration and consequent conversion of kinetic to thermal energy by the KHI could be a viable power for observed Ly$\alpha$ blobs (Dijkstra & Loeb, 2009; Goerdt et al., 2012; Wisotzki et al., 2018).

The analyses in Mandelker et al. (2016), Padnos et al. (2018) and Mandelker et al. (2019), while extremely comprehensive, considered the KHI using ideal hydrodynamics, and thus neglected potentially important influences on stream behavior from physical effects such as radiative cooling, thermal conduction, gravity and magnetic fields. In this paper, we present an idealized model of magnetized cold streams in order to understand the impact of magnetic fields on cold streams feeding galaxies. This inclusion of magnetic fields in the analysis is motivated by widespread observations of dynamically important magnetic fields in a large variety of astrophysical systems (see Durrer & Neronov, 2013; Han, 2017; Subramanian, 2019, for reviews on generation, evolution and observation of magnetic fields).

Numerical simulations of galaxy formation have shown that a small-scale dynamo can amplify magnetic seed fields in idealised setups to observed field strengths (Wang & Abel, 2009; Pakmor & Springel, 2013; Rieder & Teyssier, 2016, 2017a; Steinwandel et al., 2019) potentially aided by a cosmic-ray driven dynamo (Hanasz et al., 2004; Pakmor et al., 2016). Cosmological simulations show that tiny magnetic seed fields grow exponentially by a small-scale dynamo driven by the gravitational collapse and by supernova feedback until it saturates at around redshift $z\simeq 4$ with a magnetic energy of about 10 per cent of the kinetic energy in the interstellar medium of the protogalaxy with typical strengths of 10 to 50 $\mu$G (Beck et al., 2012; Pakmor et al., 2014; Marinacci et al., 2015; Pakmor et al., 2017; Rieder & Teyssier, 2017b; Pakmor et al., 2018).

Amplification takes place in two phases: in the kinematic phase the magnetic field grows exponentially, primarily at the smallest nonresistive scale. This is followed by the nonlinear phase, which transfers the magnetic energy toward larger scales until the dynamo saturates on the turbulent forcing scale (Schober et al., 2013). Exponential growth occurs irrespective of the seeding mechanism, which can range from primordial magnetogenesis to battery processes during the proto-galaxy formation (modelled by injecting dipole-shaped magnetic fields in Beck et al., 2013). After redshift $z\simeq 6$, strong galactic outflows driven by supernova and active galactic nucleus feedback transport metals and magnetic fields into the CGM and cosmological filaments so that the magnetization reaches levels of 0.1–1 $\mu$G close to galactic halos in simulations (Vazza et al., 2015; Marinacci et al., 2018; Nelson et al., 2019).

Observationally, there is little known about magnetic field strengths in cosmic filaments, cold streams and the CGM of massive halos, in particular at the peak of cosmic star formation around $z\simeq 2$. Using strong gravitational lensing of polarized background quasars by galaxies enables probing magnetic fields in cosmologically distant galaxies. Using differential polarization properties (Faraday rotation and fractional polarization) of such a lensing system at $z=0.44$, Mao et al. (2017) detected coherent $\mu$G magnetic fields in the lensing disk galaxy. The observation of radio halos in nearby edge-on galaxies provides direct evidence that galactic outflows transport magnetic fields several kpcs into the halo (Tüllmann et al., 2000; Heesen et al., 2018; Stein et al., 2019), and motivates our study of how magnetic fields impact cold streams feeding galaxies.

The paper is divided as follows: we present idealized models of cold streams in two- and three dimensions in Section 2. We analyze these idealized models using linear theory in Section 3, and find that fields modify the KHI in the supersonic regime. Remarkably, we find that magnetic fields can increase the growth rate of the supersonic KHI with respect to the hydrodynamic case. The non-linear dynamics of cold streams is then studied in Section 4 by performing numerical simulations. We present two-dimensional (2D) simulations of cold streams in Section 4.1 and proceed to study the more realistic, three-dimensional (3D) models in Section 4.2. Here we find that magnetic fields can suppress mixing of cold streams with the CGM during the nonlinear evolution of the instability. We assess the differences found between 2D and 3D in Section 4.3 and discuss the astrophysical implications of our study in Section 5, in particular how magnetic fields could allow streams that are $\sim 2-8$ thinner to reach the central galaxy. We conclude in Section 6 by summarizing and pointing to future work.

2 Idealized, Magnetized, Cold Stream Model

We study cold streams using both a 2D and a 3D model. The stream is a slab of width $L$ in the 2D model while it is a cylinder of diameter $L$ (radius $R=L/2$) in the 3D model, see Fig. 1. Both models consider a cold, dense stream moving at supersonic speed through the ambient hot, dilute CGM. They thus act as idealized models of filamentary accretion of cold gas in the CGM, a physical picture which is illustrated in fig. 1 in Tumlinson et al. (2017).

We take the stream density to be $\rho_{\mathrm{s}}=(1+\delta)\rho_{0}$ where $\rho_{0}$ is the density of the CGM with $\delta>0$ for a dense stream. It is assumed that the stream initially is in pressure equilibrium with the CGM which has pressure $p_{0}$. This implies that the stream temperature is $T_{\mathrm{s}}=T_{0}/(1+\delta)$ where $T_{0}$ is the CGM temperature. The stream is thus colder and denser than the CGM. As a consequence, the speed of adiabatic sound waves in the stream, $c_{\mathrm{ad},\mathrm{s}}=\sqrt{\gamma p_{0}/\rho_{\mathrm{s}}}$ where $\gamma=5/3$ is the adiabatic index, is lower in the stream than in the CGM by a factor $1/\sqrt{1+\delta}$. The abrupt change in sound speed at the surface of the stream means that it can act as waveguide for sound waves. Such sound waves are called body modes because they permeate the entire body of the stream and have an amplitude which decays exponentially with distance outside the stream. Body modes are distinguished from surface waves which are localized at the surface of the stream and do not penetrate far into it. This distinction between surface and body modes is important for understanding the physics of the supersonic KHI (Payne & Cohn, 1985; Hardee & Norman, 1988; Mandelker et al., 2016).

The key difference between our models of cold streams and the ones studied by Mandelker et al. (2016), Padnos et al. (2018) and Mandelker et al. (2019) is that we include a magnetic field and treat the system with ideal MHD (see e.g. Freidberg 2014). We consider a homogeneous magnetic field with magnitude $B$ which is oriented along the stream direction. The introduction of a magnetic field changes the wave properties of the fluid, i.e., sound waves are replaced by the compressive slow and fast magnetosonic waves and the shear Alfvén wave. The properties of these MHD waves depend on whether the propagation direction is parallel, perpendicular, or inclined with respect to the direction of the magnetic field (see e.g. Spruit 2013 for an introduction to MHD waves). The new characteristic wave speed introduced by the magnetic field is the Alfvén speed

[TABLE]

where $\mu_{0}$ is the vacuum permeability. As the Alfvén speed depends on the gas density, it is lower in the stream than in the CGM in our model (similarly to sound waves, by a factor $1/\sqrt{1+\delta}$). The stream can therefore also act as a waveguide in the MHD limit. MHD waveguide models have a long history in modeling space plasmas in the solar system (see e.g. Mazur et al. 2010 and references therein). A key property of waveguides is that they have resonant modes, and this leads to resonances in the growth rate of the KHI.

The waves that arise in an MHD model of the stream can still be divided into body modes, that have a large amplitude inside the stream, and surface modes, which have their amplitude localized at the surface of the stream. The difference with respect to hydrodynamics is that the waves are further classified as e.g. fast or slow body modes and fast and slow surface modes. These classifications are carefully outlined in Edwin & Roberts (1982) for 2D slabs and in Edwin & Roberts (1983) for 3D cylinders. Their studies considered stationary systems but can still help us gain some intuition for the dynamics of a moving stream.

The stream moves through the CGM with supersonic speed, $V$, and we define the sonic Mach number of the stream velocity with respect to the CGM111Because of the difference in sound speeds in the CGM and in the stream, $\mathcal{M}_{0}$ differs from the Mach number defined with respect to the stream itself, $\mathcal{M}_{\mathrm{s}}=V/c_{\mathrm{ad,s}}$, by $\sqrt{\delta+1}$., as

[TABLE]

In our simplified model, a stream moving at speed $V$ through a stationary medium, is physically equivalent to a stationary stream embedded in a medium moving at speed $V$. This is a consequence of the system being invariant under a Galilean transformation. Moving to this frame of reference, the stream can be regarded as a waveguide with moving boundaries. Body modes, also known as reflective modes, are primarily localized inside of the stream and gain energy by tapping into the free energy of the moving boundaries. This can be explicitly shown by calculating the reflectivity and transmission coefficients of waves as they impinge on the stream boundary. Such calculations show that the amplitude of both transmitted and reflected waves can be higher than the amplitude of the incident wave (Payne & Cohn 1985 for a 3D cylinder, Hardee & Norman 1988 for a slab). This energy gain at each reflection gives rise to exponential growth of wave amplitude, i.e., an instability. Ingoing and outgoing waves can constructively interfere at specific wavelengths, called resonant wavelengths, and the growth rate peaks for these waves. Transmitted and reflected waves have nearly equal propagation angles when such resonances occur and the resonant propagation angle can be related to the Mach number of the stream (Payne & Cohn 1985, see also Mandelker et al. 2016).

The waveguide analogy is not perfect and only describes the linear dynamics of the stream. In particular, the boundaries of the waveguide are not immovable and not perfectly reflecting, and the surfaces of the stream will eventually be disrupted by the large shear velocities present there. This could occur in the form of the classical KHI which is localized at the surfaces. For supersonic flow, however, a stabilization of surface modes occurs (Landau, 1944; Mandelker et al., 2016). The instability that arises due to body modes can therefore be the dominant one when the flow is supersonic.

The key parameter, which determines the importance of the magnetic field, is the plasma-$\beta$

[TABLE]

where $p_{\mathrm{mag}}=B^{2}/(2\mu_{0})$ is the magnetic pressure and $c=\sqrt{p/\rho}$ is the isothermal sound speed. Since we assume an initially constant pressure and magnetic field strength ($B=\mathrm{const.}$ and $p=\mathrm{const.}$ in both stream and background), and $c$ and $\varv_{\mathrm{a}}$ have the same density dependence, $\beta$ is also initially constant in space.

The idealized model outlined above of cold streams feeding galaxies at high redshift, contains three dimensionless parameters, $\rho_{\mathrm{s}}/\rho_{0}=1+\delta$, $\mathcal{M}_{0}$, and $\beta$. For a $10^{12}M_{\odot}$ dark matter halo at redshift $z=2$, estimated values for the density and Mach number of the stream are $\rho_{\mathrm{s}}/\rho_{0}\approx 10-100$ and $\mathcal{M}_{0}\approx 0.75-2.25$ (Mandelker et al., 2019). We take $\rho_{\mathrm{s}}/\rho_{0}=50$ and $\mathcal{M}_{0}=2$ throughout the main body of the paper and vary only222A parameter study of growth rates for the 2D slab, where we vary $\rho_{\mathrm{s}}/\rho_{0}$ and $\mathcal{M}_{0}$, is included in Appendix A.4. the magnetic field strength, i.e., $\beta$. We consider two different regimes for the magnetic field strength, i.e., $\beta=1$, where the magnetic pressure is as large as the thermal pressure, and $\beta=10$ where the magnetic pressure is 10 per cent of the thermal pressure. These models are compared to a purely hydrodynamic reference model ($\beta^{-1}=0$, no magnetic field).

The transition in values of density and velocity between stream and CGM have so far been described as discontinuous. A discontinuous step function however leads to growth rates that diverge with wavenumber which prohibits convergence of numerical simulations (Robertson et al., 2010; McNally et al., 2012; Lecoanet et al., 2016; Berlok & Pfrommer, 2019). We therefore use a smooth function, i.e., hyperbolic tangent profiles which have a smooth transition on a length scale $\sim 2a$, where $a$ is the smoothing parameter. This introduces a fourth dimensionless parameters, $a/L$. Note that this is not just a numerical necessity, as a smooth profile is also physically motivated due to turbulent mixing and diffusion. A high value of $a$ leads to a significant reduction of the growth rate of the KHI while a low value makes it difficult for simulations to converge (Berlok & Pfrommer, 2019; Mandelker et al., 2019). We choose a smoothing length $a/L=0.05$ which is a compromise between these two extremes.

3 Linear theory

We analyze the linear stability properties of the models by linearizing the ideal MHD equations in Cartesian and cylindrical geometry. We use the linear theory as an aid in understanding the results of the Athena++ simulations that we present in Sections 4.1 and 4.2. The details of the linear stability analysis are presented in Appendix A, here we provide a brief outline of the steps involved.

The key assumption for the linear theory is that the deviations from equilibrium have the form $f(z)\mathrm{e}^{-\mathrm{i}\omega t+\mathrm{i}kx}$ for the 2D model and $f(r)\mathrm{e}^{-\mathrm{i}\omega t+\mathrm{i}kz+\mathrm{i}m\phi}$ for the 3D model. Here, $\omega$ is a complex frequency and $\sigma=-\mathrm{Im}(\omega)>0$ is the growth rate of the KHI. Linearizing the ideal MHD equations with this ansatz, one finds that the resulting equations constitute one-dimensional eigenvalue problems. We solve these eigenvalue problems by discretizing the equations with two different types of rational Chebyshev polynomials (Boyd, 1987a, b, 2000) in the direction perpendicular to the stream. The procedure is performed with psecas333psecas, Pseudo-Spectral Eigenvalue Calculator with an Automated Solver, is freely available online. which automates some of the steps involved (Berlok & Pfrommer, 2019).

Given $\rho_{\mathrm{s}}/\rho_{0}=1+\delta$, $\mathcal{M}_{0}$, $a/L$, $\beta$ and a wavenumber $k=2\pi/\lambda$ parallel to the stream (and for the cylinder, an azimuthal mode number, $m$), psecas returns the eigenvalues and eigenvectors of the eigenvalue problem. Here $\lambda$ is the wavelength and the eigenvalues tells us how fast the KHI will grow. The eigenvectors can be used to construct linear solutions for all the perturbed variables, i.e., the density, temperature, magnetic field components etc.

3.1 2D stream

We consider a 2D, dense and cold stream of gas moving at supersonic speed through a hot and dilute background. Due to the abrupt change in velocity at the surfaces of the stream, the stream is likely to be unstable to the KHI. We are interested in understanding whether magnetic fields change this picture. Do they stabilize the KHI? Do they change the non-linear evolution and prevent the cold stream of gas from mixing with the CGM? Before introducing the extra complexity of a magnetic field, we briefly outline how the KHI in a supersonic, hydrodynamic stream already differs from the textbook picture of the KHI that many readers are likely to have in mind.

For very high velocities, i.e., supersonic motions, the KHI changes it behavior. In the classic textbook picture where a fluid moving to the right is separated from a fluid moving to the left (see e.g. Chandrasekhar 1961), the KHI arises at the interface between the two fluids. If the motion is sufficiently supersonic, disturbances on the interface are removed by the fast flow before they can grow. The precise velocity at which this stabilization occurs depends on the density contrast, $\delta$, and is given by equation 22 in Mandelker et al. (2016) who generalized the result of Landau (1944).

In a supersonic stream, however, a different version of the KHI can arise in which waves that are primarily localized inside the body of the stream grow in amplitude as they are reflected back and forth between the inside surfaces. In 2D hydrodynamics, the body modes are the dominant type of instability in the supersonic regime (Mandelker et al., 2016).

Let us now introduce a magnetic field, which is aligned with the stream. In the classic textbook picture (with two separated fluids), this system has been analyzed in Chandrasekhar (1961) at subsonic velocities by considering the incompressible limit of the MHD equations. Chandrasekhar (1961) found that the KHI occurs at the interface of the two fluids (i.e., it is a surface mode) and that the magnetic field stabilizes the KHI if the magnetic field is sufficiently strong, that is if the Alfvén speed is larger than the gas velocity.444The same criterion also applies to the case of a thin magnetic layer that separates the two counter-flowing fluids where the magnetic layer can suppress instabilities on scales significantly larger than its thickness (Dursi, 2007).

The magnetized KHI for two fluids moving at supersonic speeds was studied in Pu & Kivelson (1983) who treated the problem using the full compressible equations. As in Chandrasekhar (1961) it was found that the flow velocity needs to exceed the Alfvén speed in order for the KHI to be unstable. Additionally, it was found that highly supersonic flows stabilize the KHI and that the solution found in this limit is given by stable, magnetosonic body modes.

In the hydrodynamic case, the supersonic stabilization only takes place when there is a single interface between two fluids. When a 2D stream is considered instead, sound waves are able to grow as body modes inside of it. We find that this picture carries over to the MHD version of the KHI. In particular, we find that surface modes are stabilized by supersonic motion and magnetic tension but that the KHI can instead take place as magnetosonic, body modes that are primarily localized inside the stream. This is in agreement with Hardee et al. (1992) who studied the magnetized slab in the context of an astrophysical, under-dense jet.

We have calculated the growth rate of the KHI in the 2D stream as a function of wavenumber for the unmagnetized, $\beta=10$ and $\beta=1$ stream (see Fig. 9 in Appendix A.1). In agreement with Chandrasekhar (1961), we find that the magnetic field inhibits the growth rate of the KHI at low wavenumbers where the modes are surface modes. At higher wavenumbers, where the modes are body modes, the growth rate however increases with increasing magnetic field strength. We have the following heuristic explanation for this surprising behavior. A strong magnetic field can increase the restoring force of compressive waves, i.e., as is well known for fast magnetosonic waves in a uniform medium. We believe that the supersonic, magnetized KHI grows faster than the unmagnetized version, because the magnetic field in a similar way increases the phase velocity of reflective modes inside the stream.

We show the fastest growing eigenmode for $\beta=1$ in Fig. 2. The magnetic field lines, shown in the first panel with grey solid lines, are tied to the fluid in ideal MHD and thus trace the fluid displacement. Opposite sides of the surface of the stream are displaced in the same direction. This type of mode is called an anti-symmetric sinusoidal mode. The mode has non-zero perturbations in pressure, density and parallel magnetic field, which shows that the instability is compressive.

3.2 3D stream

We analyze how the 3D cylindrical stream differs from the 2D, rectangular stream. The primary difference is, of course, the geometry and that the extra dimension allows non-axisymmetric disturbances. The linear theory for our 3D cold stream model is closely related to the stability theory for astrophysical jets, which can also be unstable to the KHI (e.g. Ferrari et al. 1981; Bodo et al. 1989; Appl & Camenzind 1992; Bodo et al. 1996). These studies analyze the vortex sheet approximation, i.e., they consider a discontinuous transition between a jet and the background, and derive analytical dispersion relations. Our linear theory considers instead a smooth transition between stream and background which requires a numerical treatment.

The eigenmodes are obtained with the aid of psecas and have the general form $f(r)\mathrm{e}^{-\mathrm{i}\omega t+\mathrm{i}kz+\mathrm{i}m\phi}$ where $f(r)$ is a radial dependence, $k$ is the wavenumber along the stream and $m$ is the azimuthal wavenumber. The azimuthal dependence is illustrated in Fig. 3 where we show the cross-section of the stream (an $xy$ slice at $z=0$) for the fastest growing modes with $m=0-5$ and $\beta=1$. These eigenmodes have been obtained with psecas using linearized equations and a procedure which is detailed in Appendix A.2. In Fig. 3, the dotted circle shows the initial cylinder ($r=R$) and the solid black lines how the stream surface is displaced by the KHI (exaggerated here for illustrative purposes). The underlying false color images show the perturbation to the $z$-component of the magnetic field.

The $m=0$ mode is a body mode with disturbances throughout the interior of the stream. In the MHD literature this mode is known as a sausage555The deformed cylinder looks like a symmetrically squeezed sausage (see e.g. fig. 7 in Fujimura & Tsuneta, 2009). mode. The $m=1$ mode is the 3D equivalent of the sinusoidal mode that we found to be the fastest growing mode for the 2D stream, see Fig. 2. This mode is also known as a kink mode in the MHD literature (e.g. Fujimura & Tsuneta, 2009).

For higher $m$ (called fluting modes), the instability appears as a mixture of surface and body modes, with a pronounced disturbance at the surface of the cylinder and disturbances penetrating gradually less into the cylinder as $m$ increases. That is, the instability gradually changes from being dominantly of the body mode type to being dominantly a surface mode as $m$ increases. These surface modes survive stabilization by the supersonic flow because the variation in $z$ has long wavelength, making the disturbance effectively subsonic666Surface modes are suppressed when a supersonic flow removes a perturbation faster than it grows (Mandelker et al., 2016). This process is faster on short spatial scales than on large spatial scales. Mandelker et al. (2016) shows explicitly how so-called fundamental modes are unstable at long wavelength for any value of $\mathcal{M}_{0}$ (and $\delta$, see their Sections 2.3.2 and 2.4.1 for a 2D slab and 3D cylinder, respectively)., and because the azimuthal variation is perpendicular to the flow and thus unaffected by magnetic tension.

The growth rate of the KHI in 3D depends on the additional parameter, the azimuthal wavenumber, $m$. We have calculated the growth rates as a function of $kL$ and $m$ in Appendix A.2, see Fig. 10. This allowed us to find the value of $kL$ which gives the fastest growth for each azimuthal wavenumber, $m$. The fastest growing mode is an $m=4$ surface mode in both the unmagnetized and the magnetized limit.

4 MHD Simulations

The computer simulations are performed with the MHD code Athena++ (White et al., 2016) on a Cartesian domain. We perform three simulations in each geometry (2D and 3D), i.e., a pure hydrodynamic simulation and MHD simulations with $\beta=10$ and $\beta=1$, for a total of 6 simulations.

The 2D simulations have size $L_{x}\times L_{z}=(10\times 200)L$ and the 3D simulations have size $L_{x}\times L_{y}\times L_{z}=(200\times 200\times 10)L$ where $L$ is the width of the stream. The stream is moving in the $x$-direction in the 2D simulations and in the $z$-direction in the 3D simulations. The long domain size along the stream is chosen in order to resolve high wavelength modes which can have significant growth rates (see Appendix A). The highly elongated domain in the perpendicular direction(s) is chosen in order to mimic an infinite domain. We employ periodic boundary conditions in the parallel direction and outflow boundary conditions in the perpendicular direction(s).

The interesting stream dynamics and disruption take place near the initial position of the stream, and covering the entire domain with uniform resolution would be wasteful and computationally unfeasible. We therefore enable the static mesh refinement capabilities of Athena++ with 10 levels of refinement in the 2D simulations and 8 levels of refinement in 3D simulations. This increases the resolution by a factor of two at each level of refinement. The 2D simulations have $\Delta{}x=\Delta z=10\times 2^{-12}L\approx 2.4\times 10^{-3}L$ inside the central $2L$ while the 3D simulations are four times coarser. This corresponds to $\sim 400$ ($\sim 100$) cells per stream width in 2D (3D).

The simulations are initialized with subsonic, random Gaussian velocity noise, with a standard deviation of $10^{-3}c_{0}$ in each (Cartesian) component. We seed every cell in the computational domain and do not restrict the disturbances to the surface of the stream. All simulations employed the second order accurate Van Leer scheme with a Courant number of 0.4 (0.3) for the 2D (3D) simulations and the HLLE (HLLC) Riemann solver for the MHD (hydrodynamic) simulations.777HLLE is more diffusive than the HLLD Riemann solver but the 2D simulations with $\beta=1$ crash when we use HLLD. We believe this happens due to the carbuncle instability(Quirk, 1994), a numerical instability that can arise with grid-aligned shocks (see discussion in Appendix C in Stone et al. 2008).

The simulations are performed in a frame in which the stream is moving at speed $V/2$ and the background (the CGM) is moving in the opposite direction at the same speed. This is of course physically equivalent to a stream moving at speed $V$ through a stationary medium. Numerically, however, the change of reference frame reduces the speed with respect to the grid, which reduces advection errors and increases the allowed time step. All analysis of simulations is performed in the natural frame, i.e., a moving stream and stationary CGM.888Note that the dense stream in our MHD simulations moves to the left, unlike our linear theory analysis, where it moves to the right. The physics results are, of course, invariant under this change of frame.

4.1 2D stream

The initial condition for the stream velocity and density in the simulations is

[TABLE]

and the magnetic field is aligned with the stream velocity, $\bm{B}=B\bm{e}_{x}$.

Fig. 5 compares the density evolution of the stream in a hydrodynamic simulation (left-hand column) with MHD simulations with $\beta=10$ (middle column) and $\beta=1$ (right-hand column). Density disturbances due to the KHI arising as body modes are seen in the top row of panels. The criss-cross pattern inside the streams at $t=20L/c_{0}$ in the hydrodynamic and $\beta=10$ simulations are signatures of the spectrum of body modes that grow inside of them (see insets). At this point in time, the simulations are about to enter the nonlinear stage in which weak oblique shocks form (Padnos et al., 2018). The criss-cross pattern arises already earlier in the linear regime, which distinguishes the pattern as a superposition of body modes that are constructively interfering. The case of $\beta=1$, which has a larger growth rate, is already in the nonlinear regime at $t=20L/c_{0}$, and the interior of the stream is being disrupted by a series of overlapping weak oblique shock and rarefaction waves. The shock waves propagate from the stream into the background medium at an angle and transfer parallel kinetic energy and momentum from stream to the background. Due to momentum conservation, an oblique rarefaction wave is generated at the boundary that propagates back into the stream. This leads to a deceleration of the stream which we discuss in more detail in Section 4.3.

Quite remarkably, the cold stream starts to disrupt on a faster time scale in the MHD simulation with $\beta=1$ than in hydrodynamic simulation. This is a surprising result as magnetic fields inhibit the KHI in the incompressible, subsonic regime (Chandrasekhar, 1961). The faster time scale for growth of the strongly magnetized KHI is predicted by the linear theory in Section 3. We believe that the increased growth rate occurs as a consequence of an increased propagation speed of waves. Physically, an increase in phase speed allows more reflections per unit time. This could explain the increased growth rate, assuming that the energy gain per reflection is roughly constant. This interpretation, is somewhat supported by analytical studies of hydrodynamic body modes where Mandelker et al. (2016) derived that the stability criterion and growth rate for body modes depends on the sum of the sound speed inside and outside of the stream (see equations 32 and 36 in Mandelker et al., 2016). Precisely how this works in the magnetized case is left for future studies.

The middle row of panels in Fig. 5 show the streams at $t=40L/c_{0}$. The streams all have clear signatures of shock waves, with sharp density discontinuities that overlap and form complex structures. The strongly magnetized simulation $\beta=1$ has acquired a bend but has otherwise retained its structural integrity. This is in contrast to the hydrodynamic and $\beta=10$ streams which have started to disrupt.

As evident from Fig. 5, the mode that ends up disrupting the stream is not necessarily the fastest growing one. The reason is two-fold. Firstly, because we excite the instability with white noise and a finite amplitude, slower growing modes might be the first to disrupt the stream if their initial amplitude is large (as this gives them a head start). Secondly, as pointed out in Padnos et al. (2018); Mandelker et al. (2019), the fastest growing mode might saturate and cease exponential growth without breaking up the stream. This happens in the $\beta=1$ simulation where the disturbance seen at $t=20L/c_{0}$ roughly corresponds to the fastest growing mode but the bending at $t=80L/c_{0}$ occurs at a much longer wavelength. Physically, the magnetic field opposes bending of the stream because bent magnetic field lines have an associated magnetic tension. This magnetic tension is larger on small scales than on large scales (the tension force scales with $k^{2}$, see Equation 19). Magnetic tension therefore prevents small wavelength bending of the stream while still allowing large scale bending (for this magnetic field strength).

The small scale mixing of stream material with the CGM is also suppressed by magnetic fields. This is most clearly seen in the insets in the lower row of panels (at $t=80L/c_{0}$). In the MHD simulations the boundary between stream material and CGM remains sharp and clearly defined. In contrast, the boundary in the hydrodynamic simulation is less clearly defined. Instead, wispy filaments of stream material are mixed into the CGM on small scales. The difference in the mixing properties is due to magnetic tension, which becomes important on small scales.

Based solely on the 2D simulations, it seems that magnetic fields could very significantly modify cold stream disruption and mixing, even for a sub-dominant magnetic field ($\beta=10$). We will see that such a strong conclusion cannot be made in 3D, where the extra degree of freedom enhances the mixing rate in both hydrodynamic and MHD simulations. A quantitative analysis of stream mixing in 2D is therefore postponed to Section 4.3 such that we can compare with the 3D simulations.

4.2 3D stream

We extend the study of a dense and cold stream of gas moving with supersonic speed through a hot, dilute background to 3D. The stream has a cylindrical shape with diameter $L$ (radius $R=L/2$) and is cold and dense inside $r<R$. We again assume a smooth transition of velocity and density now given by

[TABLE]

where $a$ is the smoothing length. The stream moves in the $z$-direction, $\bm{\varv}=\varv(r)\bm{e}_{z}$ and the stream direction is aligned with a constant magnetic field, $\bm{B}=B\bm{e}_{z}$.

We show density snapshots from the three (hydro, $\beta=10$ and $\beta=1$) 3D simulations in Fig. 5. The main panels show cuts down the middle of the stream ($xz$-slices with $y=0$) with the perpendicular extent limited to the interesting $x\in[-2.5,4.5]L$. The insets show cross-sections of the stream ($xy$-slices at $z=0$). The cross-section at $t=0$ is shown beneath the color bar on the same scale as the rest of the panels.

The first row of panels compares the three simulations at $t=15L/c_{0}$. The disturbances in density take place at the surfaces of the streams. This differs from the 2D simulations where the disturbances arose as body modes inside the stream. We can understand this difference between 2D and 3D simulations by using results from linear theory. The body modes were dominant in the 2D simulations because surface modes with surface variation along the stream direction are suppressed by supersonic flows. In 3D, surface modes with azimuthal surface variation are however not suppressed, and are predicted to grow faster than the body modes (see Figs. 3 and 10).

The 3D simulations have much faster growth and mixing than the 2D simulations. This is evident in the second row of panels where the streams have begun their disintegration at $t=20L/c_{0}$. The boundary between stream and CGM is much less clearly defined than in the corresponding 2D simulations, even for the $\beta=1$ simulation. The key difference is that the 3D simulations allow azimuthal motions that are unaffected by magnetic tension. As a result, the lower row of panels, at $t=50L/c_{0}$, reveals that even the $\beta=1$ is eventually disrupted and mixed into the CGM. We compare the quantitative differences in stream mixing in the following section.

4.3 Comparison between 2D and 3D

In all cases, some form of the KHI (i.e., arising as surface or body modes, magnetized or unmagnetized) leads to exponential growth of kinetic energy in the perpendicular direction. The streams initially have total, parallel kinetic energy

[TABLE]

in 2D and

[TABLE]

in 3D. We show the time evolution of the perpendicular kinetic energy, relative to this initial parallel kinetic energy, in Fig. 7. For both 2D and 3D simulations, the initial exponential growth is faster in the $\beta=1$ simulations than in the hydrodynamic simulations, as expected from the linear theory.

We seed the simulations with random Gaussian velocity components, which initially decay before growing exponentially. This initial transient in perpendicular kinetic energy, seen at $t\lesssim 10L/c_{0}$, differs between hydrodynamic and MHD simulations, with a more severe reduction in the initial energy in the MHD simulations. This is because a fraction of the kinetic energy of velocity fluctuations is shared with the magnetic field on all length scales via the induction equation, reducing the initial kinetic energy more severely in comparison to the hydrodynamic case, which can dissipate this energy only on small scales. The $\beta=10$ simulation therefore displays a delay in the evolution.999The $\beta=1$ KHI has a larger growth rate than the hydrodynamic KHI and the $\beta=1$ simulation quickly catches up despite the larger initial decay. We explore an alternative seeding mechanism in Appendix A.3 where we use the fastest growing eigenmode to initialize simulations. In order to compare the different simulations on an equal footing, we have shifted the definition of $t=0$ by $3L/c_{0}$ ($5L/c_{0}$) for the 2D (3D) $\beta=10$ simulation in all analysis and figures. These shifts cause the simulations to approximately exhibit the same perpendicular kinetic energy at $t=10L/c_{0}$.

We have seen that the stream disintegrates faster in 3D than in 2D because the extra degree of freedom allows non-axisymmetric fluting modes with $m>0$ to destabilize the stream. Visual inspection of Figs. 5 and 5 also shows that the mixing occurs much faster in 3D and that the suppression of mixing by magnetic fields is less efficient. We now perform a more quantitative analysis in order to understand the consequences for the feeding of cold gas to massive galaxies at high redshift.

The stream initially has density $\rho_{\mathrm{s}}=50\rho_{0}$. A stream fully mixed with the CGM would give a uniform medium with density $\approx 1.24\rho_{0}$ ($\approx 1.001\rho_{0}$) in the 2D (3D) simulations. We study stream mixing with an observationally motivated criterion and calculate the time evolution of the total mass of gas with density above a certain cut-off, $\rho_{\mathrm{c}}$. This quantity is shown in Fig. 7 for various values of the cut-off, i.e., $\rho_{\mathrm{c}}/\rho_{0}=25,20,10,5$ and 2. The values have been normalized to their value at $t=0$, for easier comparison between 2D and 3D.

The 2D results (upper row in Fig. 7) show clear suppression of mixing by magnetic field, with gas with $\rho>10\rho_{0}$ totaling 97 per cent of the mass of the initial stream mass at the end of the $\beta=1$ simulation. With $\beta=10$, this number is 86 per cent, which is still significantly higher than the 57 per cent found for the hydrodynamic simulation.

The 3D simulations have more mixing in general and less suppression of mixing by the magnetic field (lower row in Fig. 7). Here we find that gas with density $\rho>10\rho_{0}$ totals 39 per cent of the initial stream mass at the end of the $\beta=1$ simulation but that this number is 0 per cent in both $\beta=10$ and hydrodynamic simulations. And although the mixing of material with $\rho>10\rho_{0}$ occurs at a slightly slower rate in the $\beta=10$ than in the hydrodynamic simulation, this difference does not seem significant.

There is however a non-trivial dependence on the value of $\rho_{\mathrm{c}}$ considered. For instance, the difference between the hydrodynamic and the $\beta=10$ simulation is quite significant when $\rho_{\mathrm{c}}/\rho_{0}=2$ or 5 is considered. In the $\beta=10$ simulation, gas with $\rho>2\rho_{0}$ still totals more than 50 per cent of the initial stream mass at the end of simulation while this fraction goes to zero already at $t\sim 80L/c_{0}$ in the hydrodynamic simulation.

This survival time criterion is based on the amount of mass that can retain densities above a certain cut-off, $\rho_{\mathrm{c}}$. It is also important to assess whether the stream material can remain cold. Mixing of stream with CGM that reaches a mean density $\rho_{\mathrm{c}}$ yields a mean temperature $T_{\mathrm{c}}=T_{0}\rho_{0}/\rho_{\mathrm{c}}$ (if we assume perfect mixing of the internal energy content and neglect heating/cooling mechanisms). The mixing temperature is lower than the CGM temperature but higher than the initial stream temperature by a factor $\rho_{\mathrm{s}}/\rho_{\mathrm{c}}$. A cutoff of $\rho_{\mathrm{c}}/\rho_{0}=5$ thus also corresponds to a factor 10 increase in temperature if the internal energy is mixed. Our simulations are unfortunately not able to capture this important aspect of the evolution as cooling and heat conduction is not included in the analysis. An important future extension of our work is therefore to include radiative cooling in order to understand the thermal history of the cold stream. The combination of magnetic fields and radiative cooling has been shown to increase the efficacy of thermal instability (Ji et al., 2018).

5 Astrophysical implications

We show the amount of mass with $\rho>5\rho_{0}$ for the 3D simulations in Fig. 8. As a result of the non-linear phase of the KHI the dense stream disintegrates and eventually mixes with the ambient CGM, lowering the mass of the dense stream. For this to happen, the conversion of parallel to perpendicular kinetic energy, $E_{\mathrm{kin},\perp}$, must have saturated. Consequently, while $E_{\mathrm{kin},\perp}$ is still growing ($t\lesssim 20L/c_{0}$) there is no mixing observed in either model. We observe that mixing of cold stream gas with the CGM occurs at an approximately constant rate for $t\gtrsim 20L/c_{0}$ but that the rate at which it proceeds depends on the magnetic field strength.

We perform a linear fit, indicated with black dashed lines, in order to estimate the mixing rates. We find a mixing rate of $-0.023\,c_{0}/L$ for the hydrodynamic simulation, which is $\sim 2$ higher than the one found in the $\beta=10$ simulation ($-0.011\,c_{0}/L$), and $\sim 8$ times higher than the $\beta=1$ simulations ($-0.003\,c_{0}/L$). Magnetic fields can thus inhibit the mixing of cold streams with the CGM. We can convert the mixing rates to half-live survival times, i.e., how long it takes mixing to cause less than half the mass to have $\rho>5\rho_{0}$. We find half-life survival times of $t_{1/2}\sim 20L/c_{0}$, $t_{1/2}\sim 45L/c_{0}$ and $t_{1/2}\sim 175L/c_{0}$, for the hydrodynamic, $\beta=10$ and $\beta=1$ simulations, respectively. The decreased mixing rates thus leads to increased survival times of magnetized streams.

This enables us to convert the estimated survival times into an estimated lower limit on the width of cold streams that can make it to the central galaxy. The cold stream in our model has $\mathcal{M}_{0}=2$, i.e., velocity $V=M\sqrt{\gamma}c_{0}$, such that the distance traveled by the stream during its half-life is

[TABLE]

with $\xi\approx 20$, 45 and 175 for the hydrodynamic, $\beta=10$ and $\beta=1$ simulations, respectively. Equating this estimate for the traversed distance with the virial radius of the halo, we find that streams with

[TABLE]

will have less than half their mass at densities $\rho>5\rho_{0}$ when they reach the central galaxy. In other words, cold streams thinner than this limit will be disintegrated by the KHI before they can reach the central galaxy. The decreased mixing rate of cold streams with the CGM can thus make it possible for thinner streams to reach the central galaxy when magnetic fields are included. We find that a magnetic field with $\beta=10$ has a significant effect while $\beta=1$ has a dramatic effect. We now estimate the required magnetic field strength in physical units.

The magnetic field strength is related to the hydrogen number density, $n_{\mathrm{H}}$, temperature of the hot, dilute CGM gas, and the plasma-$\beta$ by

[TABLE]

For $n_{\mathrm{H}}=10^{-4}\mathrm{cm}^{-3}$ and $T=10^{6}$ K, characteristic values for the CGM (Tumlinson et al., 2017), Equation (12) yields $B\approx 0.8\mu\mathrm{G}$ for $\beta=1$ and $B\approx 0.25\mu\mathrm{G}$ for $\beta=10$. Magnetic field strengths larger than $0.25\mu$G are therefore predicted to modify the mixing rate of cold streams with the CGM, their survival times, and consequently the minimum widths of streams which can make it to the central galaxy.

Whether such magnetic field strengths are found in filaments or the CGM is currently unknown (see also the discussion in the introduction). The magnetic fields present in the CGM are thought to be generated in the galactic disk and subsequently transported into the CGM (Marinacci et al., 2018; Nelson et al., 2019). Observations of radio halos show equipartition magnetic fields of order 7–10 $\mu$G at distances of 6 kpc from the galaxies’ midplane (Stein et al., 2019). While these values are in line with synchrotron modelling of advection-dominated, accelerating galactic winds (Miskolczi et al., 2019), advection-dominated models with a constant wind speed yield values that are smaller by a factor of 2.5 (Heesen et al., 2018). Cosmological simulations of galaxy formation predict magnetic field strengths in the required range at a radial distance of 30 kpc in a Milky Way-like disk galaxy at $z=2$ (see fig. 6 in Pakmor et al. 2014).

The introduction of magnetic fields in our models also required setting a topology. As for the magnetic field strength, to date this cannot be observationally inferred. We have assumed a simple configuration with a uniform magnetic field of constant field strength aligned with the stream both inside and outside of the stream. This choice is supported by cosmological simulations that found that magnetic fields are aligned with filaments in the cosmic web (Brüggen et al., 2005). The choice is not well motivated for the CGM in which the magnetic field is likely to be turbulent and tangled. Whether our results are robust to changes in these assumptions should be the target of future studies. For now we note that McCourt et al. (2015) found a reduced mixing of accelerating clouds containing tangled magnetic fields.

We have argued that the magnetic field strength inside the stream changes its magnetoacoustic response with the result that body modes have an increased growth rate. For the cylindrical stream, we have seen that non-axisymmetric surface modes have the dominant growth rate. The value of the magnetic field strength at the surface of the stream is therefore potentially more important than its value inside the stream. The associated magnetic field tension can have the power to stabilize surface modes of the KHI. While the magnetic field might not have high values everywhere in the CGM, the magnetic field strength could be highly enhanced at the stream surface due to magnetic draping which inevitably occurs when an object is moving super-Alfvénically through a magnetized medium (see e.g. Dursi & Pfrommer 2008; Pfrommer & Dursi 2010; McCourt et al. 2015). Super-Alfvénic motion causes the ambient magnetic field lines to accumulate at the stagnation point on the surface of a cold cloud as the field is unable to adjust to the approaching body for causality reasons. A steady state is reached once the magnetic energy density of the draping layer matches the ram pressure of the incoming wind: the dynamically important magnetic pressure pushes field lines over the body at the rate of which new field lines are accumulated in the draping layer. Provided the Alfvén speed is larger than the gas velocity, the thin magnetic draping layer that separates the two counter-flowing fluids can suppress instabilities on scales significantly larger than its thickness (Dursi, 2007; Dursi & Pfrommer, 2008).

6 Conclusion

Star forming galaxies in the early Universe have been suggested to sustain their star formation via filamentary accretion of cold gas from the cosmic web (Dekel & Birnboim, 2006; Dekel et al., 2009). A key question for this model is whether such cold streams of gas can penetrate all the way down to the central galaxy, where new fuel for star formation is needed, or whether the streams will be disrupted and heated during their propagation through the hot, virialized halo gas. Cosmological simulations have not yet attained the required spatial resolution to address this question directly (Sparre et al., 2019), and such simulations do not agree on whether the cold stream will survive the journey through the CGM (Nelson et al., 2013).

This has motivated idealized studies of cold streams, such as the ones by Mandelker et al. (2016), Padnos et al. (2018) and Mandelker et al. (2019), where instability mechanisms can be more easily understood and the streams can be properly resolved. These studies considered the cold streams using ideal hydrodynamics. Here we extend the analysis to ideal MHD by introducing magnetic fields in idealized simulations of cold streams feeding massive halos in the early Universe.

We find using linear theory that a strong magnetic field ($\beta=1$) inside the stream changes its magnetoacoustic response and enhances the growth of the KHI in the form of reflective waves that are primarily localized inside the stream (so-called reflective or body modes). While a strong magnetic field has long been known to suppress surface modes of the KHI (Chandrasekhar, 1961), a strong magnetic field can also enhance the effective growth rate because body modes are an important instability mechanism for a super sonic stream.

Although a magnetic field can enhance the growth rate of the KHI, we find that a strong magnetic field outside the stream can make it difficult for the stream to be disrupted in the non-linear phase of the KHI. This effect is clearly seen with a uniform magnetic field, everywhere aligned with the stream. In 2D, any perpendicular motion of the gas is associated with a buildup of magnetic tension due to bending of magnetic field lines. We show in Fig. 5 how a stream is able to almost completely retain its integrity when there is an ambient, uniform magnetic field with $\beta=1$. When the magnetic field strength is lower, $\beta=10$, the effect is lowered enough that the stream breaks up. The magnetic field is however still strong enough to significantly enhance the time scale on which the dense gas is fully mixed with the CGM.

In 3D, the picture is slightly changed. The introduction of a third dimension allows disturbances that vary around the surface of a cylindrical stream (and not just along it, as for the 2D version). Since the magnetic field is assumed to be aligned with the stream, these motions are not associated with any magnetic tension. This means that azimuthal surface modes (disturbances localized on the surface of the cylinder, so-called fluting modes) become important in the more realistic 3D treatment. Additionally, and unlike longitudinal surface modes, they are not stabilized by the supersonic flow. The conclusion that azimuthal modes become dominant in 3D thus carries over from the hydrodynamic case studied in Mandelker et al. (2019).

As for the 2D case, the magnetic field is able to increase the growth rates of the instabilities, again, it seems, by changing the magnetoacoustic response of the stream. We find that the most unstable mode for our parameters is the $m=4$ fluting mode and that the growth rate is increased by a factor $\sim 1.4$ with respect to the hydrodynamic case for a strong magnetic field ($\beta=1$).

Although the effect is not as dramatic as in 2D, the magnetic field is also able to suppress mixing of the cold stream with the ambient CGM in 3D (see Fig. 5 for a visualization of the non-linear phases of the KHI and Fig. 7 for a detailed comparison between mixing in 2D and 3D). We find that magnetic fields with $\beta^{-1}=0.1-1$ can increase the survival time of cold streams to be $\sim$2–8 times longer than in our hydrodynamic simulations (see Fig. 8 for a summary of the $\beta$-dependence in 3D). The survival time can be used to estimate the minimum width a stream can have without being disrupted before reaching the central galaxy. We find that streams $\sim 2-8$ times thinner can reach the central galaxy if the magnetic field strength has $\beta^{-1}=0.1-1$, which corresponds to $\sim 0.3-0.8\mu$G for characteristic CGM values of density and temperature (see Section 5 for details).

Our study has assumed a very idealized picture of a cold stream propagating through the CGM. A more realistic treatment would include turbulent magnetic fields in the CGM, radiative cooling, thermal conduction, self-gravity of the stream (Aung et al., 2019) and external gravity from the dark matter halo, which leads to a varying CGM density as the stream penetrates deeper into the halo. Nevertheless, our idealized model of a magnetized cold stream shows that magnetic fields could potentially be very important for the dynamical evolution of cold streams feeding galaxies at high redshift. This motivates including magnetic fields in both future idealized simulations and in cosmological simulations with enhanced CGM resolution.

Acknowledgments

We thank the referee, Yuval Birnboim, for a detailed and insightful report which helped us improve the manuscript. We thank Martin Sparre for sharing his knowledge on CGM physics. TB and CP acknowledge support by the European Research Council under ERC-CoG grant CRAGSMAN-646955.

Appendix A Linear theory and code verification

The equations of ideal MHD (Freidberg, 2014), i.e., the mass continuity, momentum, induction and entropy equations, are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

in SI units. Here $\bm{a}\bm{b}$ is the dyadic product of vectors $\bm{a}$ and $\bm{b}$, $\rho$ is the mass density, $\bm{\varv}$ is the mean fluid velocity, $p$ is the thermal pressure, $\bm{B}$ is the magnetic field with local direction $b$, $\mu_{0}$ is the magnetic permeability and $\gamma=5/3$ is the adiabatic index.

The nonlinear dynamics of cold streams, including fragmentation and mixing with the CGM, is studied in the main body of the paper by solving Equations (13)-(16) numerically with the aid of the MHD code Athena++ (White et al., 2016).

The linear dynamics of cold streams are governed by simpler, linearized equations of ideal MHD. We present these linearized equations in 2D and 3D in Sections A.1 and A.2 and solve them with the aid of Psecas (Berlok & Pfrommer, 2019). The obtained linear solutions are then used in Section A.3 to verify that our Athena++ setup is able to accurately capture the linear dynamics of the KHI.

A.1 2D slab in Cartesian geometry

For the magnetized KHI, the linearized equations have previously been presented in Cartesian geometry (Berlok & Pfrommer, 2019). For our specific setup, where we neglect heat conduction and viscosity as well as spatial variation in pressure and magnetic field strength, Equations 17 to 21 in Berlok & Pfrommer (2019) reduce to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $\delta A$ is the perturbed vector potential which is related to the perturbed magnetic field by $\delta\bm{B}=\bm{\nabla}\bm{\times}(\delta A\bm{e}_{y})$. The background profiles for $\varv$ and $\rho$ are given by Equation (4) and (5), respectively. Equations (17) to (21) is an eigenvalue problem, which we can solve with the rational Chebyshev polynomial grid in Psecas. We assume as boundary conditions that all perturbations approach zero as $z\rightarrow\pm\infty$. We use this procedure to calculate the growth rate of the KHI in the hydrodynamic limit ($\beta^{-1}=0$), for an intermediate field strength ($\beta=10$) and for a strong magnetic field ($\beta=1$). The calculated growth rates are shown in Fig. 9 as a function of wavenumber. At low wavenumbers, where the KHI appears as surface modes, the magnetic field suppresses the instability. This is in line with the intuition obtained from incompressible theory (Chandrasekhar, 1961), who found that the magnetic field appears as a surface tension term and suppresses the instability. At higher wavenumber, the KHI appears as body modes. These are waves that that have their dominant amplitude inside the stream and grow in amplitude by tapping into the kinetic energy of the stream (Payne & Cohn, 1985). For the body modes, the magnetic field increases the growth rate of the KHI with respect to the hydrodynamic case.

A.2 3D cylinder in polar coordinates

The linearized equations for the cylindrical stream with aligned background magnetic field and velocity, i.e., $\bm{\varv}=\varv(r)\bm{e}_{z}$ and $\bm{B}=B(r)\bm{e}_{z}$, can be written in cylindrical coordinates as

[TABLE]

[TABLE]

[TABLE]

where the differential operator $\mathcal{D}_{r}$ is

[TABLE]

and $k$ and $m$ are the longitudinal and azimuthal wavenumbers, respectively. The perturbation to the magnetic field $\delta\bm{B}$, is related to $\delta\bm{b}$ by $\delta\bm{b}=\delta\bm{B}/B$ and we have used the $\nabla\cdot\delta\bm{b}=0$ constraint,

[TABLE]

to simplify Equation (28). The background profiles for $\varv$ and $\rho$ are given in Equations (6) and (7), respectively.

We use psecas with the semi-infinite rational Chebyshev grid to solve Equations (22) to (28) on the semi-infinite interval $r\in[0,\infty]$. This grid has the advantage that it is not necessary to explicitly impose boundary conditions at $r=0$ and $r\rightarrow\infty$ (Boyd, 1987b, 2000). We present the growth rate as a function of longitudinal wavenumber, $k$, for azimuthal wavenumbers $m=0$–$5$ in Fig. 10. We again consider three different magnetizations (hydro, $\beta=10$ and $\beta=1$) and find that the $\beta=1$ cylinder has the highest growth rates. As in 2D this also appears to be related to the presence of body modes. While the lower-order azimuthal modes ($m=0$ and 1) stretches across the radial extent of the stream, the faster growing higher-order azimuthal modes exhibit a finite depth as measured from the surface. With increasing azimuthal wave number $m$ the solutions progressively attain the character of surface modes (see Fig. 3). We also find that the fastest growing eigenmode has $m=4$, independent of magnetization. A compilation of the fastest growing longitudinal wavenumber for each value of $m$ and $\beta=1$ is presented in Table 1. These are the modes shown in Fig. 3.

A.3 Code verification

We use the linear solutions obtained in Sections A.1 and A.2 to verify our Athena++ setup. We consider the fastest growing mode in the $\beta=1$ simulations. For the 2D setup, this mode has $k_{\mathrm{max}}L=3.81069$, $\sigma L/c_{0}=0.59086$ and is shown in Fig. 2. Using the psecas solutions we seed all physical variables in an Athena++ simulation. Details of such a procedure has been previously discussed in Berlok & Pfrommer (2019) for the double periodic KHI of Lecoanet et al. (2016) where simulations where initialized with a Fourier series for each component of the linear solution. Here we instead read in the coefficients for the rational Chebyshev polynomial expansions of the linear solution and use those to construct the initial conditions in the problem generator in Athena++. This yields essentially perfect exponential growth of the instability, as evident in Fig. 11. The growth of perturbations is here tracked by calculating the mean of the absolute difference between the components and their equilibrium values. The growth rate, measured by performing an exponential fit, agrees to within 1 % of the theoretical estimate. The domain size of this simulation is $L_{x}\times L_{z}$ with $L_{x}=2\pi/k_{\mathrm{max}}$ and $L_{z}=20L_{x}$, i.e., it is set up in a way such that the fastest growing mode fits exactly in the longitudinal direction. The simulation has static refinement with 7 levels in the perpendicular direction, corresponding to a resolution of $\sim 310$ cells per stream width inside the central $z\in[-L_{x},L_{x}]$.

We repeat this procedure for a 3D simulation of the KHI in a cylindrical cold stream. We seed the fastest growing $m=4$ mode with $\beta=1$ and choose the domain size to exactly fit the mode and to be much larger in the perpendicular direction, i.e., $L_{z}=2\pi/k_{\mathrm{max}}$ and $L_{x}=L_{y}=20L_{z}$. We use the coefficients of the linear solutions in terms of rational Chebyshev polynomials on a semi-infinite domain to initialize the simulation. The linear solutions are given in cylindrical coordinates while we use a Cartesian grid in Athena++. We convert the solution from cylindrical to Cartesian coordinates using a standard coordinate transformation, i.e., $r=\sqrt{x^{2}+y^{2}}$, $\phi=\mathrm{atan2}(y,x)$, $\delta b_{x}=\cos(\phi)\,\delta b_{r}-\sin(\phi)\,\delta b_{\phi}$, $\delta b_{y}=\sin(\phi)\,\delta b_{r}+\cos(\phi)\,\delta b_{\phi}$ and similarly for $\delta\varv_{x}$ and $\delta\varv_{y}$. The magnetic field is initialized directly and not via a magnetic vector potential. This introduces a numerical divergence of the magnetic field, which could potentially become a problem for the constrained transport scheme in Athena++. The direct initialization however suffices for this linear test, and we obtain excellent exponential growth, see Fig. 12. The 3D simulation has 6 levels of static refinement and $\sim 140$ cells per stream width inside the central $(x,y)\in[-L_{z},L_{z}]\times[-L_{z},L_{z}]$.

We have also verified that the deviations from the background in the simulations agree with the eigenmodes obtained with psecas. We present the background deviations after they have grown by a factor of $\sim 10$ in amplitude in the top rows of Figs. 14 and 14. These are almost indistinguishable from the psecas eigenmodes which are shown in the lower rows of Figs. 14 and 14. We have calculated the discrepancy, using as a measure the mean of the absolute difference between simulation and theory inside the domain shown in the figures. This number is $\sim 1$ % of the amplitude of the eigenmodes. The tests verify that our Athena++ setups are able to accurately model the supersonic version of the magnetized KHI.

A.4 Parameter study of 2D slab

We present a parameter study for the 2D slab stream. The main text considers a dense, cold stream with $\mathcal{M}_{0}=2$ and $\rho_{\mathrm{s}}/\rho_{0}=50$ (that is, $\delta=49$) and compares a hydrodynamic study with MHD studies with $\beta=1$ and 10. In Fig. 15 we consider the same magnetizations but extend the study to three different density contrasts ($\delta=9$, 49 and 99) and two different Mach numbers ($\mathcal{M}_{0}=1$ and 2). In this figure, the middle panel in the lower row is identical to Fig. 9.

We find that magnetic fields suppress the KHI at low density contrasts and flow speeds. For $\mathcal{M}_{0}=1$, $\beta=1$ leads to full suppression when $\rho_{\mathrm{s}}/\rho_{0}=10$ but only to inhibited growth when $\rho_{\mathrm{s}}/\rho_{0}=50$ or 100. At higher flow speeds, $\mathcal{M}_{0}=2$, the $\beta=1$ growth rate is decreased with respect to the hydrodynamic case for $\rho_{\mathrm{s}}/\rho_{0}=10$ only. At higher density contrasts, $\rho_{\mathrm{s}}/\rho_{0}=50$ and 100, the $\beta=1$ calculation instead yields growth rates that are increased with respect to the hydrodynamic case. The magnetic field-induced increase in growth rate is only found at high wavenumbers where the KHI can occur as body modes. At low wavenumbers, where the KHI appears as surface modes, the magnetic field always inhibits the instability via magnetic tension. Hydrodynamic body modes become dominant when $\mathcal{M}_{0}$ and $\delta$ are large (Mandelker et al., 2016), and this is also the regime where MHD body modes become important. Hydrodynamic body modes have a growth rate which depends on the sound wave speeds inside and outside the stream (see equation H8 in Mandelker et al. 2016). We therefore believe the increase in growth rate of body modes occurs because the magnetic field increases the phase speed of waves compared with hydrodynamic sound waves. The dispersion relations for a magnetized slab presented in Hardee et al. (1992) and for a magnetized cylinder presented in Appl & Camenzind (1992), although derived in the vortex sheet approximation and with under-dense jets in mind, might be able to provide additional insights into the behavior of supersonic cold streams and help with this interpretation of our findings. We leave such comparisons for future studies.

Finally, we find that the wavelength at which the body modes are resonant changes with magnetic field strength. This is particularly evident in the solutions for $\beta=10$ and hydrodynamics with $\rho_{\mathrm{s}}/\rho_{0}=10$, $\mathcal{M}_{0}=2$ (the lower left panel) where the $\beta=10$ resonant peaks are slightly shifted to the right. We again speculate that these shifts are due to the difference in propagation speeds of waves in unmagnetized and magnetized fluids.