A&A 406, 37-41 (2003)
DOI: 10.1051/0004-6361:20030778

The Kelvin-Helmholtz instability in an expanding universe and its effect on dark matter

K. A. Malik¹ - D. R. Matravers²

1 - GRECO, Institut d'Astrophysique de Paris, CNRS, 98bis Boulevard Arago, 75014 Paris, France
2 - Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 2EG, UK

Received 16 January 2003 / Accepted 12 May 2003

Abstract
We extend the Kelvin-Helmholtz instability to an expanding background. We study the evolution of a non-viscous irrotational fluid and find that for wavelengths much smaller than the Hubble scale small perturbations of the fluid are unstable for wavenumbers larger than a critical value. We then apply this result in the early universe, treating cold dark matter as a classical fluid with vanishing background pressure.

Key words: instabilities - cosmology: theory - cosmology: dark matter

1 Introduction

In the standard cosmological model quantum fluctuations are magnified from microscopic to cosmological scales during a period of exponential expansion (inflation) and source curvature perturbations with a scale-invariant or nearly scale-invariant spectrum (Kolb & Turner 1990; Liddle & Lyth 2000). After inflation these perturbations act as seeds for the large scale structure of the universe: first, beginning during radiation domination, non-baryonic dark matter "falls'' into the potential wells formed during inflation and later, during matter domination but after decoupling, radiation and baryonic matter follow.

In this picture radiation and baryonic matter can "fall'' into the potential wells only after decoupling since before that their Jeans length is larger than the Hubble radius. After decoupling perturbations whose wavelength is larger than the Jeans-length grow, but on scales smaller than the Jeans-length they oscillate, leaving imprints in the CMB that can be observed today. On very small scales structure is wiped out through the tight coupling of the baryons to the photons until decoupling (Silk 1968). Unlike in baryonic matter and radiation, perturbations in the non-baryonic dark matter can start growing as soon as their wavelength is smaller than the horizon, since the Jeans scale is negligible for non-interacting dark matter. The non-baryonic dark matter can therefore cluster and "fall'' into the potential wells much earlier than the other matter.

In fluid dynamics the study of instabilities (Chandrasekhar 1981) is well established and has also been extended to the astrophysical context (see Inogamov 1999, and references therein). The Kelvin-Helmholtz instability has been known for nearly 150 years (Chandrasekhar 1981). In its simplest form two large fluid volumes, separated by an interface, move in opposite directions. Neglecting the stabilizing effects of density and gravity gradients this setup is unstable for all non-zero velocities and wavenumbers.

We study the initial stages in the evolution of small perturbations in an irrotational, non-viscous fluid due to a Kelvin-Helmholtz instability in an expanding universe. We show that even in such a background these instabilities occur. We then apply these results in the early universe. Two large colliding regions of cold dark matter might serve as an example: the Kelvin-Helmholtz instability can arise at the contact surfaces of the two regions. We model the dark matter as a fluid with vanishing background pressure in concordance with the standard cold dark matter models (Kolb & Turner 1990; Liddle & Lyth 2000). Depending on when this happens in the evolution of the background the perturbations arising from the instability grow either exponentially or linearly on scales smaller than a critical wavenumber which we determine. On the basis of these results we argue that it will be interesting to test these perturbations in simulations because they could have relevance to the dark matter cusp and halo problems and to anisotropies in the CMB.

In the next section we give the governing equations for the fluid. In Sect. 3 we specify the model geometry and the general setup. After giving the background solution in Sect. 4, we derive the perturbed equations in Sect. 5 and analyse the stability of the system. In Sect. 6 we then calculate the critical wavenumbers for which instability can occur. We apply our results to standard cold dark matter in the final section and discuss possible implications.

   
2 Governing equations

We are interested in the dynamics of a fluid in the early universe on scales much smaller than the horizon size, i.e. $\lambda \ll H$. We can therefore use a Newtonian approximation on an expanding background where only the background expansion is described by a flat FRW model (Kolb & Turner 1990; Liddle & Lyth 2000).

We assume that the fluid is non-viscous and irrotational, its dynamics therefore governed by the Euler equation. The fluid does not have to be the only or dominant matter component in the model of the universe, we merely require that it does not interact with the other fluids such as radiation. As outlined above we envisage this fluid to describe standard cold dark matter, which has vanishing background pressure. We shall nevertheless keep the discussion in this and the following sections applicable to the general case of a perfect fluid and restrict the analysis to cold dark matter only when we consider a particular example in the final section.

The Euler equation in an expanding background is given by (Peebles 1980)

$$\frac{\partial}{\partial t} u^i+H u^i+u^k\frac{\partial}{a\p... ...tial}{a\partial x^i} P +\frac{\partial}{a\partial x^i}\phi=0 ,$$ (1)

where uⁱ is the three-velocity of the fluid, P is the pressure, $\rho$ is the energy density of the fluid, $\phi$ is the gravitational potential, a=a(t) is the scale factor and $H\equiv\dot a / a$ is the Hubble rate. The energy conservation equation in an expanding background is given by

$$\frac{\partial}{\partial t} \rho +3H \rho + \frac{\partial}{a\partial x^k}\left(\rho u^k\right)=0 .$$ (2)

The condition for the fluid to be irrotational is,

$$\frac{\partial}{\partial x^i} u^i =0 .$$ (3)

   
3 General formulation

We can split the variables describing our problem into a time-dependent background part and time- and space-dependent first order perturbations,

$$\psi=\psi(t)+\delta\psi(t,x^i),$$ (4)

where $\psi\equiv\rho, \phi, P, u^i$.

The geometry of our model is outlined in Fig. 1: two large volumes of fluid are separated by an interface at x³=0 and move in opposite directions with velocities $U_{\rm {a}}$ and $U_{\rm {b}}$above and below the interface, respectively. We allow for a perturbation of the fluid interface separating the two flows which is denoted by $\zeta=\zeta(t, x^1, x^2)$.

The velocity can also be split according to Eq. (4), uⁱ = Uⁱ+vⁱ, where Uⁱ is the background part and vⁱ the perturbation. Without loss of generality we can choose our coordinate system such that the x¹-coordinate is aligned with the direction of the background flow as in Fig. 1. The background flow is therefore separated by the x³=0 plane. The background velocity Uⁱ has then only one non-zero component U¹=U(t), but for the velocity perturbation we still have vⁱ=vⁱ(t,xⁱ).

We neglect gravitational back reaction on the perturbations and therefore set $\delta\phi\equiv 0$, i.e. the gravitational potential $\phi$ only has a time-dependent zeroth order part.

	Figure 1: The configuration: fluid moving to the left with velocity $U_{\rm {b}}$ is separated from fluid moving to the right with velocity $U_{\rm {a}}$by the perturbed interface $\zeta$.
Open with DEXTER

   
4 Background solutions

For the background described in the previous section we get from Eqs. (1) and (2) the simple background Euler and energy conservation equations

  

$$\dot U +HU =0,$$ (5)

$$\dot\rho+3H\rho=0,$$ (6)

which have the solutions

$$U=U_0 \left(\frac{a}{a_0}\right)^{-1} , \qquad \rho=\rho_0 \left(\frac{a}{a_0}\right)^{-3} ,$$ (7)

where U₀ and $\rho_0$ are constants of integration such that at t=t₀ we have U=U₀ and $\rho=\rho_0$.

The background expansion is described by a FRW model. We therefore get a power law solution for the scale factor (Kolb & Turner 1990; Liddle & Lyth 2000)

$$a=a_0~\left(\frac{t}{t_0}\right)^\gamma ,$$ (8)

where the constant of integration is chosen such that a=a₀ at t=t₀. The exponent $\gamma$ characterizes the different epochs of the early and late universe: ${\gamma = 1/2}$ during radiation domination and $\gamma=2/3$ during matter domination. The evolution of the scale factor is governed by the total background energy density. The solution (8) is therefore consistent with our model in the matter dominated regime and in the radiation dominated regime despite the fact that the dark matter fluid always scales as $\rho \propto a^{-3}$.

   
5 Stability analysis

In order to analyze the stability of our model we first linearize the governing equations and then look for growing mode solutions.

Linearizing around the background described in Sect. 3 we get from Eq. (1) the linearized perturbed Euler equation to first order

$$\frac{\partial}{\partial t} v^i+H v^i+U^1\frac{\partial}{a\p... ...^i +\frac{1}{\rho}\frac{\partial}{a\partial x^i} \delta P =0 ,$$ (9)

and the energy conservation Eq. (2) becomes

$$\dot{\delta\rho}+3H\delta\rho+U^1\frac{1}{a}\frac{\partial}{... ...ta\rho +\rho_0\frac{1}{a}\frac{\partial}{\partial x^k}v^k =0 .$$ (10)

The condition for the flow to be irrotational, Eq. (3), becomes

$$\frac{\partial}{\partial x^i} v^i =0 .$$ (11)

Substituting Eq. (11) into Eq. (10), we see that the velocity perturbations decouple from the density perturbations. To analyse the Kelvin-Helmholtz instability it is therefore sufficient to study Eqs. (9) and (11).

The perturbation in the position of the interface, given by $x^3=\zeta(t,x^1,x^2)$, is a first order quantity and has to be included in our analysis. The velocity of the fluid interface in the 3-direction is given by

$$v^3 = a\left[\frac{\partial\zeta}{\partial t} +H~ \zeta +u^1... ...artial x^1 } +u^2\frac{\partial\zeta}{a\partial x^2} \right] ,$$ (12)

where the RHS is simply the material derivative of the physical fluid interface position. Linearizing around the background described in Sect. 3

$$v^3 = a\left[\frac{\partial\zeta}{\partial t} +H~ \zeta +U^1\frac{\partial\zeta}{a\partial x^1 } \right]\cdot$$ (13)

We now make an ansatz for the first order quantities, extending Chandrasekhar (Chandrasekhar 1981),

$$\delta Q(t,x^i) \propto {\rm e}^{i\left( k_1x^1+k_2x^2\right) +\omega(t)} ,$$ (14)

where k₁ and k₂ are the components of the wavevector and $\omega$ is the time dependent frequency. With this ansatz growth of the perturbations and hence possible instability will set in if $\omega(t) >0$. Using the ansatz (14) the linearized Euler Eq. (9) become

 

$$\left(\dot\omega+H+iU^1\frac{k_1}{a}\right)v^{\alpha} -\frac{i}{\rho}\frac{k_{\alpha}}{a} \delta P ,$$ =

$$\left(\dot\omega+H+iU^1\frac{k_1}{a}\right)v^3 -\frac{1}{a\rho} \frac{\partial}{\partial x^3}\delta P ,$$ = (15)

where $\alpha=1,2$, whereas the continuity Eq. (11) and the linearized kinematic condition Eq. (13) give, respectively,

  

$$i\left(k_1v^1+k_2v^2\right)+\frac{\partial}{\partial x^3}v^3 = 0 ,$$

$$v^3=\left(\dot\omega+H+iU^1\frac{k_1}{a}\right)a\zeta.$$ (16)

The system of Eqs. (15) and (16) can be solved for the velocity perturbation v³,

$$\left(\dot\omega+H+iU^1\frac{k_1}{a}\right) \left(\frac{\partial^2}{{\partial x^3}^2}-{\bar k}^2\right) v^3=0 ,$$ (17)

where ${\bar k}^2=k_1^2+k_2^2$. We see from Eq. (17) that the x³ dependent part of the solution to v³ has to be $\propto {\rm e}^{\pm \bar{k} x^3}$. The kinematic condition (16) shows that $v^3/\left(\dot\omega+H+iU^1\frac{k_1}{a}\right)$ has to be continuous across the interface. The solution for v³ is therefore given by

 

$$v^3_{\rm a} A\left(\dot\omega+H+iU^1_{\rm a}\frac{k_1}{a}\right) {\rm e}^{-\bar k x^3} \qquad {\rm {for}} \; x^3 >0 ,$$ =

$$v^3_{\rm b} A\left(\dot\omega+H+iU^1_{\rm b}\frac{k_1}{a}\right) {\rm e}^{\bar k x^3} \qquad {\rm {for}} \; x^3 < 0,$$ = (18)

where we have used the condition that the velocity has to stay finite as $x^3\to \pm \infty$.

We now integrate Eq. (17) across the interface from $\zeta-\epsilon$ to $\zeta+\epsilon$ and then let $\epsilon \to 0$to get

 

$$\left(\dot\omega+H+iU^1_{\rm {a}}\frac{k_1}{a}\right) \frac{\part... ..._{\rm {b}}\frac{k_1}{a}\right) \frac{\partial}{\partial{x^3}}v^3_{\rm {b}} =0 .$$ (19)

Combining Eqs. (18) and (19) finally gives an equation for the time derivative of $\omega$,

$$\dot\omega=-H-\frac{1}{2}\frac{k_1}{a} \Big[ i\left(U_{\rm a... ...}\right) \pm \left\vert U_{\rm b}-U_{\rm a}\right\vert \Big] .$$ (20)

Since we are studying the growing or decaying behavior of the solution we are only interested in the real part of $\omega$. We can therefore neglect the complex term in the square brackets in Eq. (20) above, because it would only contribute an oscillatory part to the solution. We also see that Eq. (20) allows two solutions depending on which sign we choose in front of the $\left\vert U_{\rm b}-U_{\rm a}\right\vert$ term. We only consider the "-'' sign in the following, since it can already be seen that because the Hubble parameter H is always positive the "+'' choice will only lead to damped solutions.

Note that we recover the standard result of Chandrasekhar (Chandrasekhar 1981), i.e. that all wavenumbers are exponentially unstable, if we set a=1.

Equation (20) can now be integrated using the solution for the background velocity Eq. (7) and the scale factor Eq. (8). We have $\left\vert U_{\rm {b}}-U_{\rm {a}}\right\vert =\bar U (t/t_0)^{-\gamma}$ where we have defined $\bar U \equiv \left\vert U_{\rm {b0}}-U_{\rm {a0}}\right\vert$. We choose as initial condition that there is no amplification at the beginning, i.e. $\omega=0$ at t=t₀.

We study the case of a radiation dominated universe, i.e.  ${\gamma = 1/2}$, and a universe with a different exponent $\gamma$ separately.

Integrating the real part of Eq. (20) we get for the case ${\gamma = 1/2}$

$$\omega=\frac{1}{2} \left(\frac{k_1 \bar U t_0}{a_0}-1\right)\ln\tilde t ,$$ (21)

where we have defined $\tilde t\equiv t/t_0$.

In the case where ${\gamma \neq 1/2}$ we get, by integrating the real part of Eq. (20),

$$\omega=-\gamma \ln\tilde t +\frac{k_1}{A\left(2\gamma-1\right)} \left[1-{\tilde t}^{1-2\gamma}\right],$$ (22)

where we have introduced the quantity

$$A\equiv\frac{2a_0}{\bar U t_0} ,$$ (23)

because it plays a crucial role in defining critical wavenumbers in the following section.

   
6 Critical wavenumber

We can now analyze the solutions we have found in the previous section and calculate the critical wavenumber beyond which the perturbations are growing.

Radiation domination, ${\gamma = 1/2}$

In the radiation dominated regime, ${\gamma = 1/2}$, we find from Eq. (21) and our ansatz (14) that the flow is linearly unstable for all wavenumbers with

$$k_1 > k_{\rm crit,rad}\equiv\frac{A}{2} \cdot$$ (24)

From Eqs. (21) and (24) we see that during radiation domination perturbations with wavenumbers $k_1 > k_{\rm crit,rad}$ only grow linearly, i.e, $\delta Q \propto t$.

Case ${\gamma \neq 1/2}$

In the case of ${\gamma \neq 1/2}$the behavior of the perturbations is more complicated than in the radiation dominated case, since here the critical wave number is time dependent.

The solution Eq. (22) has two parts with different behaviors: the logarithmic part will always lead to a decrease in $\omega$, whereas the second term behaves differently depending on whether $\gamma$ is smaller or larger than 1/2. In the case ${\gamma <1/2}$ the second term is always increasing and will dominate after a finite time for large enough wavenumbers. For ${\gamma > 1/2}$ the second term in Eq. (22) will first increase and then decrease but can lead to instability if it dominates for a long enough period. From Eq. (22) the critical wavenumber for which small perturbations begin to grow is

$$k_1> k_{\rm crit,\gamma}\equiv \gamma A \frac{\left(2\gamma-1\right)\ln {\tilde t}}{1-{\tilde t}^{1-2\gamma}}\cdot$$ (25)

Note that the limit $\tilde t \to 1$ of the RHS of Eq. (25) above is well defined, $\lim_{\tilde t\to 1}\ln{\tilde t}/(1-{\tilde t}^{1-2\gamma})=1$.

We will now discuss the two cases ${\gamma <1/2}$ and ${\gamma > 1/2}$ separately.

Case ${\gamma <1/2}$

The case ${\gamma <1/2}$ is particularly simple. Using the relation (Abramowitz & Stegun 1972)

$$\ln x \le n(x^{1/n}-1)~ \qquad {\rm for} \; x,n>0 ,$$ (26)

we see that the time dependent part of Eq. (25) is decreasing and we therefore get, as a conservative value for the unstable wavenumber,

$$k_{\rm crit,L}> \gamma A .$$ (27)

For wave numbers satisfying this condition $\omega>0$ and the perturbations grow exponentially with time as ${\sim} {\exp}({t^{1-2\gamma}})$ and so the flow is exponentially unstable to small perturbations.

Case ${\gamma > 1/2}$

From Eq. (25) we see that for large $\tilde t$ the time dependent term $\ln{\tilde t}/(1-{\tilde t}^{1-2\gamma})$tends to $\ln \tilde t$. Hence, in this case, the critical wave number increases with time. In order to get a conservative estimate of the critical wavenumber we evaluate Eq. (25) at the time when a given perturbation is largest.

The perturbations reach their maximum value at

$$t_{\rm {max}}=t_0 \left( \frac{k_1}{\gamma A} \right)^\frac{1}{2\gamma-1} .$$ (28)

We can now evaluate Eq. (25) at $t=t_{\rm {max}}$ and get an implicit equation for the critical wavenumber for ${\gamma > 1/2}$

$$\frac{k_{\rm crit,G}}{\gamma A}- \ln\left( \frac{k_{\rm crit,G}}{\gamma A} \right)=1.$$ (29)

Solving Eq. (29) numerically shows that the error in neglecting the logarithmic term is $\sim$0.1%. We can therefore neglect it and get as an estimate that the perturbations will grow for wavenumbers

$$k_1 \gtrsim \gamma A .$$ (30)

When the perturbations reach their maximum they have grown as $\delta Q(t_{\rm {max}})=\delta Q(t_0) {\rm e}^{\omega_{\rm {max}}}$, where $\omega_{\rm {max}}=\omega(t_{\rm {max}})$. Note that we have chosen $\omega(t_0)=0$. We therefore find that the perturbations have grown at time $t_{\rm {max}}$ by a factor of ${\rm e}^{\omega_{\rm {max}}}$, where $\omega_{\rm {max}}$ is given by

$$\omega_{\rm {max}}=\frac{\gamma}{2\gamma-1} \left[ \frac{k_1}{\gamma A} -\ln\left(\frac{k_1}{\gamma A}\right) -1\right] .$$ (31)

From Eq. (26) we see that in Eq. (31) the square bracket is dominated by the first term and hence the perturbations are amplified by a factor of ${\sim} {\rm e}^{k_1/A}$.

Our model covers the irrotational, non-viscous case. The inclusion of viscosity is even more complicated: for small Reynolds numbers it tends to dampen the instability, whereas for large Reynolds numbers the development of turbulence becomes more likely. Although the simplest CDM model does not include viscosity there are models, like neutralino dark matter studied in (Hofmann et al. 2001), which are viscous. We note these instabilities will probably not be detected in models based on collision-free Boltzmann distributions and that the fluid model we use will break down at some small scale dependent on the dark matter properties.

   
7 Discussion and conclusions

In standard fluid dynamics the Kelvin-Helmholtz instability in the incompressible case develops for all wavenumbers k>0. However in an expanding universe we have shown that this is no longer the case and that there are critical wave numbers which separate stable and unstable domains. These wave numbers depend on the expansion of the background and the background velocity. For wave numbers larger than the critical wavenumber, $k>k_{\rm {crit}}$, small perturbations will grow.

More specifically, for a scale factor exponent ${\gamma > 1/2}$perturbations will grow for a wavenumber-dependent time. In the case ${\gamma = 1/2}$, i.e. during radiation domination, perturbations grow linearly for all times and for ${\gamma <1/2}$ they grow exponentially for all times.

As an illustration we look at the case of cold dark matter at the beginning of matter domination, $\gamma=2/3$. The physical setting could be two large volumes of cold dark matter colliding or some non-homogeneous velocity flows in the early universe. In the first case the contact surface would be prone to the Kelvin-Helmholtz instability and in the second the instability would occur inside the flow. From Eq. (23) and (30) the physical critical wavenumber at the time of matter and radiation equality is given by

$$\frac{k_{\rm {crit}}}{a_{\rm {eq}}} =\frac{2\gamma}{\bar U t_0}\frac{a_0}{a_{\rm {eq}}} ,$$ (32)

where $a_{\rm {eq}}$ is the scale factor at equality. We choose t₀ to be the present age of the universe, t₀=H₀^-1, where $H_0=100~\rm km~s^{-1}~Mpc^{-1} h$, ${a_0}/{a_{\rm {eq}}}=24~000~\Omega_0h^2$(see Kolb & Turner 1990; Liddle & Lyth 2000), with $\Omega_0\sim 1$ and the characteristic velocity is $\bar U\sim 100~\rm km~s^{-1}$ (see Green 2002; Klypin et al. 2001). With these numbers the critical scale below which perturbations are amplified is

$$\lambda_{\rm {crit}} \sim 1 ~\rm kpc .$$ (33)

Although detailed numerical analysis will be necessary in order to make definite quantitative predictions of the effect of the instability on the density distribution, we can nevertheless speculate that it will lead to a redistribution of power from scales ${\sim} \lambda_{\rm crit}$ to much smaller scales.

A significant feature of this example is that the instability occurs at such a small scale. This is important because it is well known that the cold dark matter model has problems at small scales. It seems to predict higher densities and more small scale structure than is observed. Modifications of the CDM have been tried, e.g. Warm Dark Matter and Self-Interacting Dark Matter but neither of these solves all the problems, see e.g. (Primack 2002; Tasitsiomi 2002) for an outline of the ideas and detailed references. While these discrepancies may be, at least partially, due to numerical and observational resolution and selection effects or the consequences of baryonic physics, there still appears to be a residual problem.

The surprising element of the Kelvin-Helmholtz effect is that it occurs even on expanding backgrounds and depending on the value of $\gamma$ the instabilities can grow exponentially. Our results indicate that further analysis using numerical methods is needed to explore the effect into the nonlinear regime.

Acknowledgements

The authors would like to thank Anne Green, Andrew Liddle, Roy Maartens, Dominik Schwarz, David Wands and David Weinberg for useful discussions. KM is supported by a Marie Curie Fellowship under the contract number HPMF-CT-2000-00981.

References

 
• Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions, 9th Edition (New York: Dover) In the text  
• Chandrasekhar, S. 1981, Hydrodynamic and Hydromagnetic Stability (Dover Publications) In the text  
• Green, A. M. 2002, Phys. Rev. D, 66, 083003 [arXiv:astro-ph/0207366] In the text  
• Hofmann, S., Schwarz, D. J., & Stöcker, H. 2001, Phys. Rev. D, 64, 083507 [arXiv:astro-ph/0104173] In the text  
• Inogamov, N. A. 1999, Astrophys. Space Phys., 10, 1-335 In the text NASA ADS  
• Klypin, A., Hoffman, Y. A., Kravtsov, V., & Gottlöber, S. [arXiv:astro-ph/0107104] In the text  
• Kolb, E. W., & Turner, M. S. 1990, The Early Universe (Addison-Wesley, Frontiers in Physics, 69), Redwood City, USA In the text  
• Liddle, A. R., & Lyth, D. H. 2000, Cosmological Inflation and Large-Scale Structure (Cambridge, UK: Cambridge University Press) In the text  
• Peebles, P. J. E. 1980, The Large Scale Structure of the Universe (Princeton, USA: Princeton University Press) In the text  
• Primack, J. R. [arXiv:astro-ph/0205391] In the text  
• Silk, J. 1968, ApJ, 151, 459 In the text NASA ADS  
• Tasitsiomi, A. [arXiv:astro-ph/0205464] In the text