The density distribution of accreting cosmic filaments as shaped by Kelvin-Helmholtz instability

TL;DR

This study uses high-resolution hydrodynamical simulations to investigate how Kelvin-Helmholtz instability influences the density distribution of cold gas in cosmic filaments, revealing broad, skewed distributions and turbulence effects relevant to galaxy evolution.

Contribution

It demonstrates the impact of Kelvin-Helmholtz instability on filament density distributions and turbulence, providing new insights into gas accretion processes in the circumgalactic medium.

Findings

Density distribution is skewed log-normal with a high-density tail at high Mach numbers.

Maximum densities correlate with accretion velocity, consistent with turbulence models.

Cosmological accretion can induce turbulence and broad density distributions in the CGM.

Abstract

Cosmic filaments play a crucial role in galaxy evolution transporting gas from the intergalactic medium into galaxies. However, little is known about the efficiency of this process and whether the gas is accreted in a homogenous or clumpy way. Recent observations suggest the presence of broad gas density distributions in the circumgalactic medium which could be related to the accretion of filaments. By means of high-resolution hydrodynamical simulations, we explore here the evolution of cold accreting filaments flowing through the hot circumgalactic medium (CGM) of high-z galaxies. In particular, we examine the nonlinear effects of Kelvin-Helmholtz instability (KHI) on the development of broad gas density distributions and on the formation of cold, dense clumps. We explore a large parameter space in filament and perturbation properties, such as, filament Mach number, initial perturbation wavelength, and thickness of the interface between the filament and the halo. We find that the time averaged density distribution of the cold gas is qualitatively consistent with a skewed log-normal probability distribution function (PDF) plus an additional component in form of a high density tail for high Mach-numbers. Our results suggest a tight correlation between the accreting velocity and the maximum densities developing in the filament which is consistent with the variance-Mach number relation for turbulence. Therefore, cosmological accretion could be a viable mechanism to produce turbulence and broad gas density distributions within the CGM.

Tables

Table 1. Table 1: Parameters used for initial conditions in our simulations. λ / R s 𝜆 subscript 𝑅 𝑠 \lambda/R_{s} is the wavelength of the perturbation at the interface between the filament and hot halo. σ / R s 𝜎 subscript 𝑅 𝑠 \sigma/R_{s} is the thickness of the interface between the two media. M b subscript 𝑀 𝑏 M_{b} is the Mach number with respect to the sound speed of the halo gas. In particular, the values for the range 0.5 - 2.1 are M b = subscript 𝑀 𝑏 absent M_{b}= 0.5, 1.0, 1.1, 1.3, 1.5, 2.1. Simulations with parameters above the horizontal line belong to our fiducial runs.

$\lambda /R_s$	$\sigma /R_s$	$M_b$
4	2%	0.5 - 2.1
1	4%	0.5 - 2.1
2	4%	0.5 - 2.1
4	4%	0.5 - 2.1
1	9%	1.1, 2.1
2	9%	1.1, 2.1
4	9%	1.1, 2.1
2	2%	1.1

Equations

$$M_b=\frac{v_{\mathrm{fil}}}{c_{\mathrm{halo}}}$$

$$\begin{split}f(x)&=f_{b}+0.25(f_{s}-f_{b})\times\\ &\bigg{[}1+\rm{tanh}\bigg{(}\frac{R_{s}-x}{\sigma}\bigg{)}\bigg{]}\bigg{[}1+\rm{tanh}\bigg{(}\frac{R_{s}+x}{\sigma}\bigg{)}\bigg{]}.\ \end{split}$$

$$P=A\rm{cos}\bigg{(}\frac{2\pi}{\lambda}y\bigg{)}\times\bigg{[}\rm{exp}\bigg{(}-\frac{(x-R_{s})^{2}}{2\Sigma^{2}}\bigg{)}+\rm{exp}\bigg{(}-\frac{(x+R_{s})^{2}}{2\Sigma^{2}}\bigg{)}\bigg{]},\$$

$$t_{KH}=\frac{\sqrt{\delta}}{M_{b}}\frac{\lambda}{c_{b}}\sim 300\Bigg{(}\frac{\lambda}{R_{s}}\Bigg{)}\Bigg{(}\frac{R_{s}}{4\textrm{kpc}}\Bigg{)}\Bigg{(}\frac{\sqrt{\delta}}{M_{b}}\Bigg{)}\Bigg{(}\frac{T}{10^{6}\textrm{K}}\Bigg{)}^{-1/2}\;\textrm{Myr}$$

$$T_{\mathrm{halo}}^{\prime }=\alpha T_{\mathrm{halo}}$$

$$T_{\mathrm{fil}}^{\prime }=\alpha T_{\mathrm{fil}}$$

$$v_{\mathrm{fil}}^{\prime }={\alpha }^{1/2}{v}_{\mathrm{fil}}$$

$$c_{\mathrm{halo}}^{\prime }={\alpha }^{1/2}{c}_{\mathrm{halo}}$$

$$R_s^{\prime }=\beta R_s$$

$$t^{\prime }=t\beta {\alpha }^{-1/2}$$

$${\rho }_{\mathrm{halo}}^{\prime }=\gamma {\rho }_{\mathrm{halo}}$$

$${\rho }_{\mathrm{fil}}^{\prime }=\gamma {\rho }_{\mathrm{fil}}.$$

$$\text{PDF}(\text{ln}\rho )d\text{ln}\rho$$

$$u$$

$$\lambda$$

$$\textrm{PDF}(\rho)d\rho=\frac{1}{\rho\sigma_{s}\sqrt{2\pi}}\textrm{exp}\Bigg{(}-\frac{(\textrm{ln}(\rho)-\mu)^{2}}{2\sigma_{s}^{2}}\Bigg{)}d\rho,$$

$$j_{\mathrm{m}}=\int {\rho }_{\mathrm{cold}}{v}_{x}dy,$$

$$PD=\max (|{\mathrm{y}}_{\mathrm{interface}}-{\mathrm{y}}_{\rho \le {\rho }_{\mathrm{thresh}}}|)$$

$${\rho }_{\mathrm{thresh}}=\frac{1}{3}{\rho }_{\mathrm{fil}}$$

$$S^2$$

Full text

††pubyear: 2018††pagerange: The density distribution of accreting cosmic filaments as shaped by Kelvin-Helmholtz instability –A

1 Introduction

Our standard cosmological $\Lambda$CDM model predicts that both dark and baryonic matter is arranged into a system of filaments, called the “cosmic web" (e.g. Bond et al., 1996). Recent observations have found extended Ly$\alpha$ emission around quasars at redshifts $z\sim 3$ that is believed to trace the densest parts of this cosmic web (e.g. Cantalupo et al., 2014; Martin et al., 2014; Hennawi et al., 2015; Borisova et al., 2016). Clustering studies (e.g. Eftekharzadeh et al., 2015) suggest that these quasars are hosted by haloes massive enough to be above the critical mass for shock heating (Birnboim & Dekel, 2003; Dekel & Birnboim, 2006) and therefore containing mostly hot gas at temperatures of $T>10^{6}$ K. Theoretical arguments and cosmological simulations (Kereš et al., 2005; Ocvirk et al., 2008; Dekel et al., 2009; Ceverino et al., 2010; Faucher-Giguère et al., 2011; van de Voort et al., 2011) suggest the presence, within these hot haloes, of cold accreting filaments possibly feeding the central galaxy and driving its evolution. However, it is unclear what happens to the cold accreting filaments along their way through the hot CGM. Ly$\alpha$ and metal emission observations of the CGM at $z>2$ suggest that the cold gas (T$\simeq 10^{4}$ K) has a broad density distribution with densities greatly exceeding $n>1\;\rm{cm}^{-3}$ (Cantalupo et al., 2014; Arrigoni Battaia et al., 2015; Borisova et al., 2016; Hennawi et al., 2015; Cai et al., 2017; Cantalupo et al., 2019). What is the origin of such broad gas density distribution within the CGM? Is cosmological accretion of cold filaments responsible for the presence of dense, cold clumps in the CGM? If properly resolved, accreting filaments might not stay laminar, as typically found in cosmological simulations, but could be subject to hydrodynamic and other instabilities that would result in the development of turbulence and fragmentation. In particular, Kelvin-Helmholtz instability (KHI) is expected to occur at the interface between the accreting filament and the hot halo (e.g. Bodo et al., 1998). Recently, a series of studies have focused on the effects of KHI on accreting filaments. These authors have developed a thorough analytical and numerical model on the stability of cold accreting filaments (Mandelker et al., 2016; Padnos et al., 2018; Mandelker et al., 2018). They have found that in some conditions KHI is expected to become highly nonlinear before the filament reaches the central galaxy. Filaments with a size in the range of 0.5 - 5 % of the virial halo radius fully disrupt prior to reaching the central galaxy, leading to powerful turbulence around the stream.

In this paper we study the nonlinear evolution of KHI in cold accreting filaments, using idealised simulations, for a purely adiabatic case. We complement Mandelker et al. (2016)’s and Padnos et al. (2018)’s work by studying the density of the filament, particularly focusing on the high densities, developing due to KHI.

In Section 2 we will describe the methods and initial setup used for the simulations. We present our results in Section 3 followed by a discussion in Section 4. Lastly, we summarise our findings in Section 5.

2 Methods

In this Section, we describe the numerical simulations used to investigate the effect of nonlinear KHI on cold accreting filaments using idealised simulations. Similar to Mandelker et al. (2016) we simulate a gas slab with high density and low temperature, representing a cold accreting filament flowing through a hot, less dense background medium, representing the halo gas, see Fig. 1. To investigate the nonlinear evolution of KHI, we place pressure perturbations on the interface between the two media.

Our simulations are run in 2D with the Eulerian code RAMSES (Teyssier, 2002). The piecewise linear reconstruction used in the Godunov scheme is the MonCen slope limiter (van Leer, 1977), and a HLLC approximate is used for the Riemann solver (Toro et al., 1994). The domain of the simulation lays in the xy plane and is a square box. The filament is centred at the middle of the box and has a radius of $R_{s}=1/12L$ where $L$ is the box size. We use periodic boundary conditions at $y=0$ and $y=L$ and outflow boundary conditions at $x=0$ and $x=L$. Outflow boundaries are chosen in order to minimise waves reflected off the box sides. Gas that crosses the boundary is lost. We have checked that and found the amount of gas lost during simulations is negligible.

The halo gas, at $|x|>R_{s}$ is initially set to have a density $\rho_{\rm{halo}}=1.3\times 10^{-4}\;\rm{cm^{-3}}$ and it is static ($v_{\rm{halo,0}}=0$). However, because our simulations do not include cooling they are scale-free with respect to length, time, temperature and density, see eq. 5 for further details.

The filament gas at $|x|<R_{s}$ is set to a density $\rho_{\rm{fil}}=\delta\rho_{\rm{halo}}$. $\delta=152$ is the density contrast between the halo and the filament and is kept fixed for all our simulations. As commonly done in the literature, we use Mach number instead of a physical velocity for the filament, also because KHI is expected to be scale free. The Mach number $M_{b}$ is unitless, and it is defined as the ratio between the filament velocity and the sound speed in the halo background:

[TABLE]

where $v_{\rm{fil}}$ is the velocity of the filament and $c_{\rm{halo}}$ is the sound speed of the halo. In setting our initial conditions we closely follow the approach of Mandelker et al. (2016).

In order to reduce numerical noise and to approximate real filaments where the boundary is not perfectly discontinuous, we smooth the velocity and density across the interface between the halo and the filament using a ramp function:

[TABLE]

This results in a smooth interface between the two fluids. $\sigma$ represents what we will call the interface thickness. Perturbations of wavelengths of similar size or smaller than $3\sigma$ are suppressed by this smooth interface (Robertson et al., 2010). We have assumed both the filament and the halo background to be ideal gases with adiabatic index $\gamma=5/3$. We initially assume pressure equilibrium between the filament and halo gas. The value of the pressure is given by the temperature and density through the ideal gas law. The temperature of the halo (filament) is $T_{\rm{halo}}=2.3\times 10^{6}\;\rm{K}$ ($T_{\rm{fil}}=1.5\times 10^{4}\;\rm{K}$). Pressure perturbations are placed along the filament interface on both sides:

[TABLE]

where $A=0.05$ is the pressure amplitude, $\Sigma=2\sigma$ is the width of the perturbation and $\lambda$ the perturbation wavelength. We have run further simulations for different values of $\Sigma$ and found no significant differences. Others have also found the exact value of $\Sigma$ has negligible impact as long as it is larger than $\sigma$.

We use a single grid with a fixed resolution. In order to keep possible resolution effects independent of varying perturbation wavelengths, we fix the number of cells per wavelength $\lambda$. The cell size, $\Delta$, is $2^{-13}$, $2^{-12}$, $2^{-11}$ for $\lambda/R_{s}=1,2,4$ respectively.

In all our fiducial calculations we set the filament radius, $R_{s}=4$ kpc and $L=48$ kpc. We have run all our simulations for 1.6 Gyr assuming $T_{\rm{fil}}=1.5\times 10^{4}\;\rm{K}$. Because the adiabatic simulations (and the development of Kelvin-Helmholtz instability) that we are performing are scale free (Chandrasekhar, 1961), this physical time may be rescaled using the relations 5 above. As a comparison, the Kelvin-Helmholtz time $t_{KH}$:

[TABLE]

for the same set of parameters and $M_{b}=1.0$, $\lambda=2R_{s}$ is 165 Myr.

We have run several simulations with different initial conditions shown in table 1 and have chosen all our initial conditions to meet the requirement $\lambda>30\sigma>60\Delta$ stated by Mandelker et al. (2016), for which they find good convergence. The range of $M_{b}$ was chosen to fit likely velocities of the infalling filament through a hot halo. $M_{b}=0.5$ ($M_{b}=2.1$) corresponds to a filament flowing at $1/4$ (1) times the freefall velocity assuming the halo is at the virial temperature.

Because of the scale-free nature of our methods, the results of the simulations will be equivalent for a set of different initial conditions related by the following transformations (at a fixed $M_{b}$, $\delta_{\rm{f}}$ and $\lambda/R_{s}$):

[TABLE]

In the above transformations, $R_{s}$ represents the radius of the filament, $T_{\rm{halo}}$ is the halo temperature, while $\alpha$, $\beta$ and $\gamma$ are rescaling parameters.

3 Results

3.1 Development of high densities during filament accretion

In this Section, we investigate whether accreting filaments into hot haloes develop high density regions during their lifetime (the time they need to reach the central galaxy, $t_{\rm{cross}}$) and if so, which properties of the filament will most likely lead to these structures. In particular, we will address the following questions: i) what are the density distributions of cold accreting filaments, ii) what are the largest densities that develop, iii) how do these high density regions evolve, and iiii) what are the morphological properties of the densest clumps. We will explore these questions by varying the properties of the filament in terms of velocity, initial perturbation wavelength and interface thickness (see Section 2).

Density distribution of accreting filaments

In Fig. 2, we present the time averaged probability distribution functions (PDFs) of the cold gas ($T<5\times 10^{4}$ K) 111 Note, that this breaks the scale invariance, though only with respect to parameter $\alpha$ in eq. 5 densities within the computational volume. The time average is performed from 2 $t_{\rm{KH}}$ until $1.6$ Gyr, that represents a period long enough to include the expected nonlinear development of KHI and is of similar order as the halo crossing time, $t_{\rm{cross}}$ (see Section 2). The densities are expressed as a density contrast ($\delta_{\rm{f}}$) with respect to the initial filament density ($\rho_{\rm{f,init}}$). There are four sets of simulations, each with its own fixed $\lambda/R_{s}$ and $\sigma/R_{s}$ and the same variety of Mach numbers (see Table 1), represented by the different solid lines in the figure. This combination of simulations allows us to understand which parameter drives the development of high densities.

We find that the time averaged PDFs of filaments with $M_{b}\approx 1$ are qualitatively consistent with skewed log-normal distributions that are overlaid for illustrative purposes with dashed lines in Figure 2. These skewed log-normals are defined as in Hopkins (2013):

[TABLE]

where $I_{1}(x)$ is the modified Bessel function of the first kind and PDF $=0$ for $u<0$. When $T\rightarrow 0$ the above distribution converges to a log-normal distribution:

[TABLE]

where $\mu$ is the average of the logarithm of the density. Simulations with other $M_{b}$ are also in qualitative agreement with skewed log-normals but show additional components. In Section 4, we will discuss the possible origin of these distributions. We find three additional components deviating from the skewed log-normal: a low density tail, a high density tail approximating a power law and a peak at $\delta_{\rm{f}}=1$. Most simulations with Mach Numbers $M_{b}>1.3$ show the broadest cold gas density distribution, deviating from the skewed log-normal with both a high and low density tail. In particular, these simulations have cold gas reaching very high (low) values of more (less) than 20 ($10^{-3}$) the initial filament density. The high density tail is most pronounced for $\lambda/R_{s}=4$, $\sigma/R_{s}=2\%$ and $M_{b}=$1.5, 2.1 (lower right figure) and seems to follow a power law. It is present in each simulation with $M_{b}>1.3$, however not always well visible in Fig. 2. This is because PDFs presented in Fig. 2 are time averaged and therefore short-lived high-density regions will not contribute substantially to the PDF. Furthermore, Fig. 3 (bottom right) displays this seemingly missing high density tail for the case of $\lambda/R_{s}=2$, $\sigma/R_{s}=4\%$ at an individual time. Further, all sets of simulations show a peak at $\delta_{\rm{f}}=1$, which corresponds to unperturbed cold gas of the filament. It is present for all simulations with $M_{b}=0.5$, and it disappears for runs with $M_{b}\approx 1.0$, $\lambda/R_{s}>1$ and then reappears when $M_{b}>1.3$. Lastly, we note that the subsonic case ($M_{b}=0.5$) results in the narrowest PDF that never exceeds $\delta_{\rm{f}}=3$ in any of our simulations.

On the left of Fig. 3, we show for illustrative purposes three simulation snapshots showing the density contrast ($\delta_{\rm{f}}$) for $M_{b}=0.5$ (top), $M_{b}=1.3$ (middle) and $M_{b}=2.1$ (bottom) at $t=4\times t_{\rm{KH}}$ for the two upper images and at $t=8\times t_{\rm{KH}}$ for $M_{b}=2.1$. The regions with density contrasts $\delta_{\rm{f}}\geq 3$ ($\delta_{\rm{f}}\geq 2$) are indicated by grey (white) contours. These simulations present very different PDFs at the time of the snapshot (Fig. 3 right). The slower moving filament ($M_{b}=0.5$) shows a flatter density as a function of vertical distance and does not develop high densities. The faster moving filament ($M_{b}=1.3$) develops higher densities at the center and has a broader PDF. The fastest filament with $M_{b}=2.1$ develops the highest densities, which are produced at the interface between the two media.

In addition to the time averaged density distribution we are also interested in the short-lived high density regions because these regions might become long lived high density regions, due to additional physical processes, which will be discussed in Section 4. Indeed, we find these regions exist for some filaments, as shown in 3. In order to capture the development of high-density regions at a given time and the general time evolution of the PDF, we decided to use as a possible proxy for these properties: the maximum density contrast $\rm{max}(\delta_{\rm{f}})$. This parameter may be sensitive to numerical noise but, as we show later on, this is not the case in our simulations (see Fig. 11).

In Fig. 4 we plot $\rm{max}(\delta_{\rm{f}})$ vs time for the same simulations as above, where $\rm{max}(\delta_{\rm{f}})$ is taken at each time step. Simulations with $M_{b}=0.5$ do not reach high densities ($\delta_{\rm{f}}>3$) for the entire duration of the simulations for all tested initial conditions. However, all simulations with $M_{b}\geq 1.3$ show density contrasts of at least 6 independently of initial condition. With individual exceptions, these high density regions are mostly relatively short-lived and vanish after $\sim$ 300 Myr. We notice that the higher the $M_{b}$ the later the high densities seem to develop, see also Fig. 6. However, there are also exceptions to this behaviour. Moreover, there are systematic differences between the sets of simulations, which will be discussed later.

In Fig. 5 we present the maximum density contrast as a function of $M_{b}$, where each symbol and color correspond to one set of parameters. In particular, the color indicates the interface thickness, $\sigma/R_{s}$, and the shape represents the wavelength $\lambda/R_{s}$. There seems to be a clear trend with increasing $M_{b}$ up to $M_{b}=1.5$, independently of the other parameters. For larger $M_{b}$, there is instead a large spread of values for $\rm{max}(\delta_{\rm{f}})$ with no clear correlation with any values of $\lambda/R_{s}$ and $\sigma/R_{s}$. In Section 4 we discuss the possible origin of this trend or lack thereof for $M_{b}>1.5$.

When and where the high densities develop

In Fig. 6 we show at which time the maximum density contrast $t_{\rm{max}(\delta_{\rm{f}})}$ occurs as a function of $M_{b}$ for a large set of parameters. In particular, we have only considered simulations that produced significant over-densities ($\rm{max}(\delta_{\rm{f}})\geq 3$) at any time, i.e. we have excluded models with $M_{b}=0.5$. Symbols and colors represent different combinations of parameters as in Fig. 5. Note that to avoid the overlap of the symbols we have slightly shifted them horizontally. To help the reader, we have joined with solid lines the symbols that correspond to the same set of $\sigma/R_{s}$ and $\lambda/R_{s}$, if simulation results are available for the entire range of $M_{b}$. As is clear from the figure, the maximum over-densities develop faster for $1.0\leq M_{b}\leq 1.3$. Once again, there is a larger spread in $t_{\rm{max}(\delta_{\rm{f}})}$ for $M_{b}=2.1$.

Now, we investigate the effect of the interface thickness $\sigma$ and $\lambda$ on $t_{\rm{max}(\delta_{\rm{f}})}$ fixing $M_{b}=1.1$. This is the typical $M_{b}$ value within the range where over-densities develop faster, as shown in Fig. 6. In particular, Fig. 7 shows how $\rm{max}(\delta_{\rm{f}})$ changes with time for different values of $\sigma$ for a constant $\lambda/R_{s}=2$ and $M_{b}=1.1$. Similarly, in Fig. 8, we show how $t_{\rm{max}(\delta_{\rm{f}})}$ changes with $\sigma$ for three different values of $\lambda$ and $M_{b}=1.1$. As shown in Fig. 7 and 8, the thinner the interface between hot halo and filament, the shorter it takes for high densities to develop. The same trend also applies to $M_{b}=2.1$ (not shown here).

Through a visual inspection of the simulation outputs we have noticed that high densities always start developing at the interface. Afterwards, the majority of these “travel" to the center of the filament and they become less dense. Once they reach the center of the filament, the densities increase again due to the interaction with overdensities coming from the opposite side. For $M_{b}\geq 2.1$ high densities do not always penetrate within the filament and disappear after some time.

In order to quantify this behaviour, we show in Fig. 9 the distance between the position of $\rm{max}(\delta_{\rm{f}})$ and the center of the filament in function of $M_{b}$ for different values of $\sigma$ and $\lambda$. As is clear from the figure, $\rm{max}(\delta_{\rm{f}})$ is either located at the interface or at the center of the filament. In particular, only perturbations with a large wavelength $\lambda/R_{s}=4$ develop the $\rm{max}(\delta_{\rm{f}})$ at the center. For $M_{b}\geq 1.3$ all $\rm{max}(\delta_{\rm{f}})$ are found at the interface. Conversely, simulations with lower $M_{b}$ produce high densities both at the center or at the interface depending on the value of $\lambda$ and $\sigma$. Because we are only tracing the maximum densities, this result does not imply that there are no overdensities in the rest of the filament. We notice that in most cases, there are high densities ($\delta_{\rm{f}}\geq 3$) developing at the center even when $\rm{max}(\delta_{\rm{f}})$ is occurring at the interface, see as an example Fig. 3 (middle).

Morphology of the high-density regions

In this Section we present the results of one simulation with $M_{b}=1.3$, $\lambda/R_{s}=2$ and $\sigma/R_{s}=4\%$ as a representative case for the development of high densities regions forming in the power law tail of the PDF. We use the maximum difference between the density of the cell corresponding to $\rm{max}(\delta_{\rm{f}})$ and the densities of its 8 neighbours (indicated as $\rm{max}(\Delta\delta_{\rm{f}})$ in the remainder of this Section). In Fig. 10 we show both $\rm{max}(\Delta\delta_{\rm{f}})$ and $\rm{max}(\delta_{\rm{f}})$ as a function of time as the red and cyan solid line, respectively. We find that $\Delta(\delta_{\rm{f}})<<(\delta_{\rm{f}})$ at any time. This suggests the seeds developing in our simulations are real and larger than a few simulation cells. Clearly, the time evolution of $\rm{max}(\Delta\delta_{\rm{f}})$ traces the corresponding evolution of $\rm{max}(\delta_{\rm{f}})$ very well.

The resolved nature of the the densest regions made up of numerous pixels can also be verified by a direct visual inspection of the simulation output at the moment in which $\rm{max}(\delta_{\rm{f}})$ is reached, as shown in Fig. 11 (see also Fig. 3). We also note that their morphology appears elongated along the direction of the filament flow. We will discuss the possible origin of this morphology in Section 4.2.

3.2 Filament disruption and accretion rates

How much cold gas from filament accretion reaches the galaxy?

Now, we investigate how much cold gas reaches the galaxy at the center of the halo. To answer this, we investigate the net cold gas mass flux ($j_{\rm{m}}$), defined as

[TABLE]

as a function of time using our fiducial suite of simulations with $\lambda/R_{s}=2$, $\sigma/R_{s}=4\%$.

In Fig. 12, we show the net mass flux of the cold gas ($T<5\times 10^{4}$ K) through a vertical slice of the simulation box as a function of time normalised to the initial value. Because we are using only one slice, the net flux varied strongly over short time scales. In order to show only the average behaviour on large time scales we have binned our results over 0.3 Gyr in time. As it is clear from the figure, the mass flux of cold gas declines monotonically for all tested cases. Filaments that are nearly sonic show the largest decrease in mass flux by 20 percent.

We have seen in Section 3.1 that filaments with $M_{b}$ in the range $1.0-1.3$ produce largest densities faster, possibly indicating a connection between the development of high densities and a reduction in accretion rate.

Do filaments disrupt during accretion?

If accreting filaments were to disrupt, this would likely have an impact on their morphology, which could be observed. In order to study the filament disruption we introduce a quantity called Perturbation Depth (PD) defined as:

[TABLE]

where $y_{\rm{interface}}$ is the initial y-location of the interface (as defined in the initial conditions) and $y_{\rho\leq\rho_{\rm{thresh}}}$ are the y-locations of simulation cells with gas densities below $\rho_{\rm{thresh}}$. In particular, the PD quantifies how deep the surrounding halo gas with a density smaller than a certain threshold penetrates into the filament. This approach is different from Padnos et al. (2018) where tracers in the fluid are used to determine whether the filament has fully disrupted. We tested different values for this threshold and verified that a value of $1/3\rho_{\rm{fil}}$ best captures how far the surrounding halo gas penetrates into the filament without contamination from low density regions in the center. In order to be less sensitive to numerical noise, we re-grid our simulation box into cells with sizes that are nine times larger than the initial ones.

In Fig. 13, we show how the normalised PD evolves with time for different $M_{b}$ and fixing $\lambda/R_{s}=2$ and $\sigma/R_{s}=4\%$. All simulations with $M_{b}\leq 1.3$ reach a PD of 0.8$R_{s}$.

In Fig. 14, we show the maximum Perturbation Depth, $\rm{max}(PD)$ reached throughout the duration of the simulation as a function of $M_{b}$ for a large range of parameters. Simulations with $M_{b}\geq 1.3$ show a small Perturbation Depth ($PD<0.6R_{s}$) for all explored parameters. This result, together with the time evolution of $j_{\rm{m}}$ presented in 12, suggests that filaments with $M_{b}\geq 1.3$ are less likely to be disrupted. We also remark that there is a preferential interface thickness, $\sigma/R_{s}=4\%$ (blue symbols in Fig. 14) where filaments tend to have a larger $\rm{max}(PD)$, e.g. they may disrupt more easily.

4 Discussion

In the previous section we have found that, when properly resolved, accreting filaments into hot haloes develop instabilities that result in: i) broad “cold gas" density distributions that are qualitatively consistent with skewed log-normal Probability Distribution Functions ii) additional components to the skewed log-normal PDF in form of high and low density tails for $M_{b}>1.3$ (where $M_{b}$ indicates the filament velocity with respect to the hot halo sound speed; see Fig. 2); iii) maximum density contrasts up to a factor of 20, depending on the value of $M_{b}$ (see Fig. 5); iv) faster development of high density regions for $1\leq M_{b}\leq 1.3$ and thinner interface regions (see Fig. 6 & 7); v) more likely filament disruption and possibly smaller accretion rates onto galaxies for $1<M_{b}<1.5$ (see Fig. 14). These results suggest that a $M_{b}\geq 1$ is a promising condition for the formation of high density structures within haloes due to filament accretion and a moderate reduction in accretion rate. In this section, we discuss the possible origins and the implications of these results.

4.1 Origin of the broad gas density distribution and fate of the filament

There is a large literature on the origin of log-normal gas density distribution through turbulent processes, traditionally applied to the Interstellar Medium (ISM) (Padoan et al., 1997; Vazquez-Semadeni, 1994; Passot & Vázquez-Semadeni, 1998; Padoan & Nordlund, 2002). In particular, previous studies suggest that a succession of compressive or expansive waves can interact with each other producing density fluctuations on different scales (that are uncorrelated and continuous) which would result in log-normal distributions (e.g. Passot & Vázquez-Semadeni, 1998; Squire & Hopkins, 2017). However, recent studies have found deviations from an exactly log-normal behaviour and argued that a skewed log-normal gas density distribution better reproduces the simulation results (Federrath et al., 2010, 2008; Hopkins, 2013; Molina et al., 2012; Nolan et al., 2015). In particular, these authors argued that a log-normal density distribution on multiple scales would indeed violate mass conservation and that deviation from this shape are therefore to be expected (e.g. Hopkins, 2013). Moreover, Nolan et al. (2015) have found a relation between the sonic Mach number ($M_{s}=v_{s}/c_{s}$) and the variance (S) of the skewed log-normal gas density distribution, in the case of non-isothermal adiabatic turbulence:

[TABLE]

where $\gamma$ is the adiabatic index and $b$ is a constant for the turbulence driving parameter that depends on the mode mixture induced by the turbulent forcing mechanism (e.g. curl-free or divergence-free). Although the precise value of $b$ is unclear when modes are mixed, several studies suggest that $b$ ranges between 1/3 and 1 (see e.g. Federrath et al., 2008, for further details).

Although our simulations are different from isotropic turbulent studies, we will use these results as a possible ansatz for the origin of the skewed log-normal gas density distribution in our case. Instead of externally driving turbulence, our simulations can in principle develop turbulence more naturally due to KHI, i.e., through the disruption of vortices via secondary instabilities, vortex pairing and nonlinear interactions (Lecoanet et al., 2016). Note that $M_{s}$ indicates the velocity of a gas parcel with respect to its own sound speed while $M_{b}$, in our study, represents the filament velocity with respect to the sound speed of the static, hot halo gas.

In order to determine whether turbulence could be the origin of a skewed log-normal PDF also for accreting filaments, we test whether the relation 10 holds in our case. To this end we need to compute the value of $M_{s}$. In particular, we examined a simulation with $M_{b}=1.3$ at a given time (see middle panel of Fig. 3) for which the PDF is well approximated by a skewed log-normal, finding values for $M_{s}$ for cells with temperatures below $5\times 10^{4}$ K. In order to derive these values, we have first computed the velocity dispersion $v_{s}$ for each cell with temperatures below $5\times 10^{4}$ K including all neighbouring cells across a set of of scales ($L_{c}$). These scales range from the thickness of the interface (20 cells) to half the filament radius (340 cells). We have then divided $v_{s}$ of each cell by its own sound speed and obtained in this way a range of $M_{s}$ that spans from 0.02 to 2.1 for all "cold" cells and scales. Note that, as expected, $M_{s}$ increases with $L_{c}$. What are the expected values of S in this case using equation 10? Using the range of $M_{s}$ above, we obtain S$=[0.002,0.3]$ (S$=[0.0013,0.2]$) for $b=1$ ($b=1/3$). This range of S is consistent with the measured dispersion of the skewed log-normal PDF in our simulation, namely $S=0.009$ (see dashed line in Fig. 3), suggesting that this component of the PDF may originate from turbulence, at least for this particular simulation.

How do these results change considering simulations with different $M_{b}$, $\lambda/R_{s}$ and interface thickness $\sigma/R_{s}$? As shown in Fig. 5, we do find a possible linear correlation between max$(\delta_{\rm{f}})$ and $M_{b}$, albeit with a large scatter at $M_{b}=2.1$, from which one could also expect a linear relation between S and $M_{b}$, at least within a range of $M_{b}$. Moreover, this result seems fairly independent of the value of $\lambda/R_{s}$ and $\sigma/R_{s}$, at least in the tested range of values. This suggests that the skewed log-normal PDF could originate from turbulence in the majority of our simulations.

While the skewed-log-normal component is consistent with turbulence as discussed above, the high-density power-law tail cannot be easily explained in this scenario. What is the origin of this component? Previous studies have suggested that KHI are unable to grow at the surface of the filament for $M_{b}>M_{\rm{crit}}=(1+\delta^{-1/3})^{3/2}$ (Mandelker et al., 2016). Note that $M_{\rm{crit}}$ is exactly 1.3 for our simulations, for which $\delta=152$ and this is the same $M_{b}$ above which we notice the appearance of high density tails in the PDFs (see Fig. 2). Although KHI are suppressed at the filament interface above $M_{\rm{crit}}$, we stil expect that instabilities could be present as a consequence of reflected sound waves (also called “body modes", see Mandelker et al., 2016). Such sound waves reflect off the oppositely located interfaces of the filament and eventually compress the gas at the filament surface. Alternatively, “body modes" may slightly alter the overall shape of the gas stream exposing small regions at the interface with the hot gas to a large velocity gradient resulting in extra compression therein. In both cases, these processes could lead to the production of small "chains" of strongly compressed and therefore very dense cold gas, as observed in our simulations at $M_{b}\geq 1.3$ (see Figs. 3, bottom left & 11). In particular, the individual dense gas clumps have the highest density at their center, which decreases with distance, resulting in a continuous density distribution that could resemble a power law.

We note that our fast flowing filaments never develop the highest densities at their center (Fig. 9) and they never disrupt, at least for the range of explored parameters. This is consistent with the picture discussed above, because “body modes" can destroy the filament only if the perturbation length ($\lambda$) is larger than $10R_{s}=\lambda_{\rm{crit}}$ (Padnos et al., 2018), which we did not consider here. On the other hand, filaments with $M_{b}<1.3$, that have sufficiently large perturbation wavelength ($\lambda>2R_{s}$) and with an interface thickness of $\sigma/R_{s}=4\%$ are most likely to disrupt. Are these filaments able to reach the central galaxies before they are disrupted? This depends on the halo crossing time, $t_{\rm{cross}}$ that in our simulations ranges between 0.3 Gyr and 1.3 Gyr for $M_{b}=2.1$, and 0.5 respectively for a halo mass of M_{\rm{vir}}=10^{12.5}$$\textrm{M}_{\odot} at redshift $z\sim 3$. Fig. 13 shows that, if a filament with radius of 4 kpc is subject to disruption, the disruption takes place at $t_{\rm{disrupt}}\sim 0.3-0.4\;\rm{Gyr}<t_{\rm{cross}}$, for this specific set of simulations. This is approximately half the filament sound crossing time, which is $t_{\rm{sound}}\sim 0.6\;\textrm{Gyr}$ for a filament with $R_{s}=4$ kpc and $M_{b}\leq 1.3$. These results are consistent with previous findings (e.g. Padnos et al., 2018). We remark that our choice of $M_{b}$ is motivated by cosmological simulations, which suggest that infalling filaments will have a constant velocity of about the halo virial velocity (Dekel et al., 2009). Depending on assumptions of the virial temperature, the halo mass and the velocity of the filament, various Mach numbers are plausible, ranging from $M_{b}=0.5$ up till 5.0.

Lastly, we comment on possible implications of our results in context of recent observations of Ly$\alpha$ nebulae around quasars that suggest a multiphase, possibly extremely clumpy intergalactic and circumgalactic medium. A possible ansatz for this density distribution is a log-normal PDF of the emitting gas of the Nebulae with dispersion $\sigma_{s}\sim 2$ defined as in eq. 7 (Cantalupo et al., 2019). Although our PDFs have a much smaller dispersion than by these observations, they may be broad enough to trigger a succession of processes which may lead to densities in line with these recent findings. These processes and their effect on the density distribution of the gas will be discussed in the next section.

4.2 Current limitations & future prospects

In this section, we discuss possible processes, not yet included in our idealised experiment, e.g. radiative cooling and gravity, or numerical limitations which could lead to a broadening or narrowing of the gas density distribution or a modification of our results.

In particular, we expect that when high densities are formed (either due to KHI or "body modes" as discussed above), radiative cooling could further enhance them which will then possibly lead to a broadening of the gas density distribution and also extend the cold clumps “lifetime". Furthermore, parts of the filament in which the densities have been enhanced by KHI may be subject to gravitational instabilities, especially when this effect is combined with radiative cooling, leading to a more pronounced broadening of the gas density distribution.

On the other hand, cooling and gravity could also prevent KHI to develop in the first place, at least at the interface, through stabilisation effects (e.g. Vietri et al., 1997). Similarly, KHI may be prevented by thermal conduction through smoothing of small scale inhomogeneities and velocity gradients (Armillotta et al., 2016). However, if KHI grow and “eddies" develop, thermal conduction may avoid subsequent KHI within these regions avoiding their fragmentation and mixing with the hot material. Note, that since we are using a grid code, there is already numerical conduction in our simulation, as determined by the cell size. Lastly, numerous studies have shown that magnetic fields have a stabilising effect on KHI (Vietri et al., 1997, and references therein), however there are little observational constraints on the strength of magnetic fields in the CGM of galaxies at $z\sim 3$.

What is the effect of numerical resolution on our results? Because of limited computational time, it is difficult to assess the numerical convergence of our results by increasing the simulation resolution. However, we have re-run some simulations at lower resolution as shown in the Appendix in order to have at least a qualitative assessment of resolution effects. In particular, we have tested resolution effects on the value of max($\delta_{\rm{f}}$), which could be expected to be particularly prone to changes in resolution, finding very limited changes.

Although performing simulations in 3D might change the results presented here (see also Mandelker et al., 2018), we think that including radiative cooling could be a more important step towards a better characterisation of the real CGM density distributions. We plan therefore to implement radiative cooling in our next set of simulations for future studies.

5 Summary & Conclusions

Our standard cosmological model predicts that cosmic filaments should play a crucial role in galaxy evolution transporting gas from the intergalactic medium into galaxies. In this work, we have explored how hydrodynamical effects such as Kelvin-Helmoltz instabilities (KHI) affect the physical properties of the accreting streams within the circumgalactic medium (CGM) of high-redshift galaxies with a particular focus on the formation and evolution of high density structures. The presence of broad gas density distribution in the cold, Ly$\alpha$ emitting phase of the CGM has been suggested by recent observations of extended emission around high-redshift quasars (e.g. Cantalupo et al., 2014; Borisova et al., 2016; Cantalupo et al., 2019).

In order to achieve the required spatial resolution to study in detail the formation and evolution of KHI, we have performed idealised simulations with the hydrodynamical code RAMSES (Teyssier, 2002) and explored a large parameter space in filament and perturbation properties, such as filament Mach number, initial perturbation wavelength, and thickness of the interface between the filament and the halo. We find that the time averaged density distribution of the cold gas is qualitatively consistent with a skewed log-normal probability distribution function (PDF) plus an additional component in form of a high density tail for simulations with $M_{b}>1.3$. We explored the possible origin of the skewed log-normal PDF and found that the broadness of the PDF is related to the value of $M_{b}$ as qualitatively expected in a turbulence-driven scenario in which the turbulence is produced by the nonlinear development of KHI. We speculate that the high-density tail of the PDF could be produced by the interaction between the hot gas and coherent density waves produced within the filament by the development of reflecting sound waves at high $M_{b}$.

In our simulations the cold gas densities can reach values as high as 20 times the initial filament densities but only for a very limited time before being dissolved. In all explored cases the simulated cold gas density distributions are much narrower than the one implied by the observations, suggesting that additional physical mechanism or processes are needed. However, our simulations suggest that cosmological accretion, in some particular cases, could be a viable mechanisms to produce the ‘seeds’ for the possible growth of higher over-densities in the CGM of high-redshift galaxies.

Acknowledgements

We thank Romain Teyssier for making his code, RAMSES publicly available and gratefully acknowledge support from the Swiss National Foundation grant PP00P2_163824. Also, I want to thank my husband for his consistent and ongoing support.

Appendix A Resolution Study

We have rerun some simulations with lower resolution in order to get an estimate of how resolution may be affecting our results. Fig. 15 (top) shows the evolution of max($\delta_{\rm{f}}$) with time for two resolutions. Although max($\delta_{\rm{f}}$) may be prone to numerical noise, we find the effect of changing resolution to be small. In in Fig. 15 (bottom) we show the PDFs for two different resolutions and likewise, we find the difference of the PDFs to be too small to be significant.

We do not find any significant or systematic dependence on resolution, especially at high densities.