Classical Waves and Instabilities Using the Minimalist Approach

1. Introduction

The study of perturbations in any physical system sheds light into the mechanisms responsible for the equilibrium, and offers a way to predict its temporal behavior as a response to any variation. Being a difficult task in general, it is important to be able to simplify each case, focusing on the most important ingredients that determine its evolution. Classical waves and instabilities correspond to such simplified structures. For example acoustic compressible perturbations are attributed to sound or magnetosonic waves and can be identified even in systems with complex dynamics. Alfvén waves are important in magnetized plasmas. When gravity or buoyancy forces are important, differences in density drive gravity waves. And by “gravity” we could describe even fictitious forces in non-inertial frames, like the effective gravity in a decelerating system, or the centrifugal. In cases of heavy over light fluid stratifications immediately the Rayleigh-Taylor instability comes to mind. Similarly, a relative velocity between two fluids triggers the Kelvin-Helmholtz instability. These are only a few examples of simplified structures, that are nevertheless crucial in order to understand more complicated and realistic settings. They are often connected to temporal behaviors in laboratory experiments or particular features observed in astrophysical systems.

The full study of the perturbation evolution in most cases can only be done through numerical simulations. The linearization of the system of equations is a simpler way to understand under which conditions the system is unstable, although it does not cover the nonlinear evolution. Normal mode analysis is a further simplification, allowed only in symmetric unperturbed states, and can be used to find the wave frequencies or the instability growth rates. This allows a simpler parametric study, or even to have analytic expressions that directly show the physics of the mechanisms involved.

Since we are discussing classical problems, there are many books, e.g., [1,2,3,4,5], that provide the most basic information. There are also a plethora of related works in the literature over the years, e.g., [6,7,8,9,10,11] to name but a few. Nevertheless, these problems continue to be an active area of research, e.g., [12,13,14,15], and there is always room – and need – to improve our understanding.

The goal of this paper is to present the methodology to find the dispersion relation in each case, following the minimalist approach introduced in Ref. [16]. This novel approach simplifies the procedure as much as possible, since it uses a single first-order ordinary differential equation, dubbed principal equation. We first derive that equation through linearization and discuss the boundary conditions that need to be satisfied (Section 2). Then we apply the method to find the dispersion relation for MHD waves (Section 3), for gravity waves and the Rayleigh-Taylor instability (Section 4), for gravito-acoustic waves (Section 5), and for the Kelvin-Helmholtz instability (Section 6). We analyze the latter more extensively, deriving new analytical expressions for the growth rate in the magnetized case, and discussing the factors that affect most the results.

2. Linear Analysis, the Principal Equation, and Boundary Conditions

The non-relativistic ideal MHD equations in Lorentz-Heaviside units are

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathit{V}\right)=0,\end{array}$$

(1)

$$\begin{array}{c}\hfill \left({\displaystyle \frac{\partial }{\partial t}}+\mathit{V}·\nabla \right)P=c_s^2\left({\displaystyle \frac{\partial }{\partial t}}+\mathit{V}·\nabla \right)\rho ,\end{array}$$

(2)

$$\begin{array}{c}\hfill \rho \left({\displaystyle \frac{\partial }{\partial t}}+\mathit{V}·\nabla \right)\mathit{V}=-\nabla P+\left(\nabla \times \mathit{B}\right)\times \mathit{B}+\rho \mathit{g},\end{array}$$

(3)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathit{B}}{\partial t}}=\nabla \times \left(\mathit{V}\times \mathit{B}\right),\end{array}$$

(4)

$$\begin{array}{c}\hfill \nabla ·\mathit{B}=0.\end{array}$$

(5)

We are interested in exploring the perturbations of a steady-state that depends only on one spatial cartesian coordinate. In particular we assume that the unperturbed state has density ${\rho }_0\left(x\right)$, pressure $P\left(x\right)$, bulk velocity ${\mathit{V}}_0=V_{0z}\left(x\right)\widehat{z}+V_{0y}\left(x\right)\widehat{y}$, magnetic field ${\mathit{B}}_0=B_{0z}\left(x\right)\widehat{z}+B_{0y}\left(x\right)\widehat{y}$, and the acceleration of gravity is constant $\mathit{g}=-g\widehat{x}$. Introducing the total pressure

$$\begin{array}{c}\hfill \mathsf{\Pi }=P+{\displaystyle \frac{B^2}{2}},\end{array}$$

(6)

the zeroth order equations are satisfied provided that

$$\begin{array}{c}\hfill {\displaystyle \frac{d{\mathsf{\Pi }}_0}{dx}}=-{\rho }_{0}g.\end{array}$$

(7)

Adding perturbations of the form

$$\begin{array}{c}\hfill \mathit{V}={\mathit{V}}_0+\left[V_{1z}\left(x\right)\widehat{z}+V_{1x}\left(x\right)\widehat{x}+V_{1y}\left(x\right)\widehat{y}\right]e^{i(k_{y}y+k_{z}z-\omega t)},\end{array}$$

(8)

$$\begin{array}{c}\hfill \mathit{B}={\mathit{B}}_0+\left[B_{1z}\left(x\right)\widehat{z}+B_{1x}\left(x\right)\widehat{x}+B_{1y}\left(x\right)\widehat{y}\right]e^{i(k_{y}y+k_{z}z-\omega t)},\end{array}$$

(9)

$$\begin{array}{c}\hfill \rho ={\rho }_0\left(x\right)+{\rho }_1\left(x\right)e^{i(k_{y}y+k_{z}z-\omega t)},\end{array}$$

(10)

$$\begin{array}{c}\hfill \mathsf{\Pi }={\mathsf{\Pi }}_0\left(x\right)+{\mathsf{\Pi }}_1\left(x\right)e^{i(k_{y}y+k_{z}z-\omega t)},\end{array}$$

(11)

defining the wavevector in the $yz$ plane

$$\begin{array}{c}\hfill {\mathit{k}}_0=k_y\widehat{y}+k_z\widehat{z},\end{array}$$

(12)

and the Doppler-shifted frequency

$$\begin{array}{c}\hfill {\omega }_0=\omega -{\mathit{k}}_0·{\mathit{V}}_0,\end{array}$$

(13)

we linearize the equations as shown in Appendix A; the main steps and related comments follow.

The linearized system reduces to two differential equations for the perturbations of the velocity in the $\widehat{x}$ direction and the total pressure (there are algebraic relations for all the other perturbations, connecting them to $V_{1x}$, ${\mathsf{\Pi }}_1$, and their derivatives). Instead of $V_{1x}$ and ${\mathsf{\Pi }}_1$ it is more convenient to use two other quantities.

The first, replacing $V_{1x}$, is the Lagrangian displacement in the $\widehat{x}$ direction. The relation between the Lagrangian displacement $\mathit{\xi }$ and the velocity perturbation results from the expression of the velocity in the perturbed location of each fluid element $\mathit{V}+(\mathit{\xi }·\nabla )\mathit{V}\approx {\mathit{V}}_0+\delta \mathit{V}+(\mathit{\xi }·\nabla ){\mathit{V}}_0$ which equals ${\mathit{V}}_0+{\displaystyle \frac{d\mathit{\xi }}{dt}}\approx {\mathit{V}}_0+{\displaystyle \frac{\partial \mathit{\xi }}{\partial t}}+({\mathit{V}}_0·\nabla )\mathit{\xi }$, yielding $\delta \mathit{V}={\displaystyle \frac{\partial \mathit{\xi }}{\partial t}}+({\mathit{V}}_0·\nabla )\mathit{\xi }-(\mathit{\xi }·\nabla ){\mathit{V}}_0=-i{\omega }_0\mathit{\xi }-{\xi }_x{\displaystyle \frac{d{\mathit{V}}_0}{dx}}$. The $\widehat{x}$ component gives ${\xi }_x=y_1\left(x\right)e^{i(k_{y}y+k_{z}z-\omega t)}$ with $V_{1x}=-i{\omega }_0{y}_1$.

The second quantity, replacing ${\mathsf{\Pi }}_1$, is the perturbation of the total pressure in the perturbed location of each fluid element $y_2={\mathsf{\Pi }}_1+y_1{\displaystyle \frac{d{\mathsf{\Pi }}_0}{dx}}={\mathsf{\Pi }}_1-{\rho }_{0}g{y}_1$.

The advantage of these replacements is that the new functions $y_1$ and $y_2$ are everywhere continuous, even at locations where the unperturbed state has contact discontinuities.

2.1. The System for $y_1$, $y_2$

The resulting $2\times 2$ system is

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{dx}}\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)+\left(\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& -f_{11}\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)=0,\end{array}$$

(14)

$$\begin{array}{c}\hfill f_{11}={\displaystyle \frac{{\rho }_{0}g}{A}}\left({\displaystyle \frac{{\rho }_0{\omega }_0^2{F}^2}{S}}-{\mathit{k}}_0^2\right),f_{12}={\displaystyle \frac{{\tilde{\kappa }}^2}{A}},f_{21}=-A-{\displaystyle \frac{{\rho }_0^2{g}^2}{A}}\left({\displaystyle \frac{F^4}{S}}-{\mathit{k}}_0^2\right),\end{array}$$

(15)

$$\begin{array}{c}\hfill F={\mathit{k}}_0·{\mathit{B}}_0,A={\rho }_0{\omega }_0^2-F^2,S={\rho }_0(A{c}_s^2+{\omega }_0^2{B}_0^2),{\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-{\mathit{k}}_0^2.\end{array}$$

(16)

To simplify the expressions we define the following quantities that have important physical meaning and effect in the resulting dispersion relations through their appearance in the array elements $f_{ij}$.

The $F={\mathit{k}}_0·{\mathit{B}}_0$ is connected to the angle between the wavevector and the unperturbed magnetic field. Its presence in a dispersion relation represents the influence of the magnetic tension which acts like a spring with restoring force per mass $-{\displaystyle \frac{F^2}{{\rho }_0}}{y}_1$. This is evident in the pure Alfvén waves (in a static homogeneous plasma without gravity) for which the displacement satisfies ${\ddot{y}}_1=-{\omega }^2{y}_1$ with $\omega ={\displaystyle \frac{{\mathit{k}}_0·{\mathit{B}}_0}{\sqrt{{\rho }_0}}}={\displaystyle \frac{F}{\sqrt{{\rho }_0}}}$.

The $A={\rho }_0\left({\omega }_0^2-{\displaystyle \frac{F^2}{{\rho }_0}}\right)$ appears in the denominator of the resulting system of differential equations, and its zeros correspond to the Alfvén waves.

The $S={\rho }_0^2\left[\left({\displaystyle \frac{B_0^2}{{\rho }_0}}+c_s^2\right){\omega }_0^2-{\displaystyle \frac{F^2}{{\rho }_0}}{c}_s^2\right]$ is related to the density perturbation in the perturbed location of each fluid element ${\rho }_1+y_1{\displaystyle \frac{d{\rho }_0}{dx}}={\rho }_0^2{\displaystyle \frac{{\omega }_0^2{y}_2+g{F}^2{y}_1}{S}}$. It is interesting to note that the incompressible limit corresponds to $S\to \infty$.

The $\tilde{\kappa }$ defined through ${\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-{\mathit{k}}_0^2$ represents the local wavenumber in the $\widehat{x}$ direction. Indeed, the latter equation is equivalent to the dispersion relation of the slow/fast magnetosonic waves ${\left({\displaystyle \frac{{\omega }_0^2}{{\mathit{k}}_t^2}}\right)}^2-\left({\displaystyle \frac{{\omega }_0^2}{{\mathit{k}}_t^2}}\right)\left({\displaystyle \frac{B_0^2}{{\rho }_0}}+c_s^2\right)+{\displaystyle \frac{B_0^2}{{\rho }_0}}{c}_s^2{cos}^2{\theta }_0=0$ if we define the total wavevector ${\mathit{k}}_t={\mathit{k}}_0+\tilde{\kappa }\widehat{x}$ and the angle ${\theta }_0$ between ${\mathit{k}}_t$ and ${\mathit{B}}_0$ through $cos{\theta }_0={\displaystyle \frac{F}{|{\mathit{k}}_t|B_0}}$.

Another way to see the connection between $\tilde{\kappa }$ and the wavenumber in the $\widehat{x}$ direction is to consider the homogeneous fluid case without gravity, in which the system (14) (with constant coefficients in that case) has solutions ${\displaystyle \frac{y_1^{\prime }}{y_1}}={\displaystyle \frac{y_2^{\prime }}{y_2}}=i\tilde{\kappa }$. The equation ${\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-{\mathit{k}}_0^2$ gives two opposite solutions for $\tilde{\kappa }$ corresponding to oppositely moving waves in the $\pm \widehat{x}$ direction. Note also that in general $\tilde{\kappa }$ is complex and its imaginary part corresponds to exponential variation of the eigenfunctions, since $e^{i\tilde{\kappa }x}=e^{-\Im \tilde{\kappa }x}{e}^{i\Re \tilde{\kappa }x}$. In the incompressible limit in which $S\to \infty$ and ${\tilde{\kappa }}^2\to -{\mathit{k}}_0^2$, it simplifies to purely exponential dependence $e^{-|k_{0}x|}$ (without sinusoidal dependence on x). More details on these waves will be discussed later; they offer a way to understand the physics of the solutions even in non-homogeneous fluids, and are directly related to the boundary conditions when the fluid extends up to large distances $x\to +\infty$ or $-\infty$.

The cases $A=0$, $S=0$, ${\omega }_0=0$, need to be studied separately, something that can be easily done since they lead to simplified equations. They correspond to waves, and there is no need to be considered in instability studies in which ${\mathit{k}}_0$ is real and $\omega$ complex.

2.2. The Principal Equation

Since the system (14) is linear only the ratio $Y={\displaystyle \frac{y_1}{y_2}}$ is uniquely defined. Using the equations of the system we directly find that Y satisfies the principal equation

$$\begin{array}{c}\hfill {\displaystyle \frac{dY}{dx}}=f_{21}{Y}^2-2f_{11}Y-f_{12}.\end{array}$$

(17)

Following the minimalist approach it is sufficient to work with the ratio Y and solve the principal equation in order to find the dispersion relation of the wave or instability. We just need to integrate this single, first order differential equation, requiring Y to be everywhere continuous and satisfying the correct boundary conditions at the extreme values of x.

Knowing Y we can find all the other perturbations from

$$\begin{array}{c}\hfill {\displaystyle \frac{y_2^{\prime }}{y_2}}=-f_{21}Y+f_{11},y_1=Y{y}_2\end{array}$$

(18)

and the relations (A28)–(A35) given in Appendix A. We emphasize though that these are not needed to find the dispersion relation.

2.3. Boundary Conditions

The linearization was based on the functions $y_1$ and $y_2$ that are everywhere continuous. The same holds for their ratio, so in case of discontinuities of the unperturbed state we simply continue the integration of the principal equation when passing from one medium to the other keeping Y continuous.

In case where the fluid ends at some extreme value of x we always know the value of Y on that end and need to use it as a boundary condition. A first example of this kind is the case of a solid wall at that end, in which case Y vanishes on the wall since the Lagrangian displacement $y_1$ vanishes. A second example is the case when outside the fluid under consideration there exist a medium with constant total pressure, as for example in a hydrodynamic atmosphere with negligible density. In that case $1/Y$ vanishes at the interface since the perturbation of the total pressure $y_2$ vanishes. (We arrive at the same result if we consider nonzero perturbation in the atmosphere, solve the problem in both sides requiring continuity of Y, and at the end take the limit of zero density outside.)

The nontrivial boundary conditions that need closer examination correspond to cases where the fluid extends to theoretically infinite distances $x\to +\infty$ or $x\to -\infty$. At these distances the unperturbed fluid is homogeneous and the solutions of the principal equation are given in Appendix B. They correspond to two oppositely moving waves in the $\widehat{x}$ direction whose amplitude has an exponential dependence on x, one increasing and the other decreasing. In order to avoid the wave with the diverging amplitude1 the solution of the principal equation should also approach one of the constant values $Y_{\pm }={\displaystyle \frac{f_{11}\pm \sqrt{f_{11}^2+f_{12}{f}_{21}}}{f_{21}}}={\displaystyle \frac{f_{12}}{-f_{11}\pm \sqrt{f_{11}^2+f_{12}{f}_{21}}}}$ that correspond to complex wavenumbers $\pm K=\pm i\sqrt{f_{11}^2+f_{12}{f}_{21}}$, respectively. Assuming that $\sqrt{f_{11}^2+f_{12}{f}_{21}}$ is the principal value of the root (with positive real part) we choose the upper sign at $x\to +\infty$ and the lower at $x\to -\infty$.

In case that $\Im K=0$ (the $f_{11}^2+f_{12}{f}_{21}$ is negative real number) the amplitude of the wave remains constant at infinity, and we can choose the solution whose sign of $\pm \Re K$ corresponds to the desired propagation.

The equations needed to apply the minimalist approach are summarized in Table 1.

Note that $\sqrt{f_{11}^2+f_{12}{f}_{21}}=\sqrt{-{\tilde{\kappa }}^2+{\displaystyle \frac{{\rho }_0^2{g}^2{\mathit{k}}_0^2}{S}}}=\left|{\mathit{k}}_0\right|\sqrt{1-{\displaystyle \frac{{\omega }_0^4/{\mathit{k}}_0^2-g^2}{\left(B_0^2/{\rho }_0+c_s^2\right){\omega }_0^2-c_s^2{F}^2/{\rho }_0}}}$. In the presence of gravity this is not constant sine the unperturbed total pressure is variable with its gradient balancing gravity. However, in cases where the fluid extends to theoretically infinite distances (in which the total pressure also reaches theoretically infinite values) it approaches a constant $\sqrt{f_{11}^2+f_{12}{f}_{21}}\to \left|{\mathit{k}}_0\right|$ corresponding to the incompressible limit.

In the absence of gravity $\sqrt{f_{11}^2+f_{12}{f}_{21}}=\sqrt{-{\tilde{\kappa }}^2}$ and without loss of generality we can choose the sign such that $\tilde{\kappa }=\pm K=\pm i\sqrt{-{\tilde{\kappa }}^2}$ (taking the principal value of the root and the upper/lower sign at $x\to \pm \infty$, respectively).

3. MHD Waves

3.1. Slow/Fast Magnetosonic Waves

For homogeneous unperturbed state and zero gravity $f_{11}=0$, $f_{12}={\displaystyle \frac{{\tilde{\kappa }}^2}{A}}$, $f_{21}=-A$, and the principal equation becomes ${\displaystyle \frac{dY}{dx}}=-A{Y}^2-{\displaystyle \frac{{\tilde{\kappa }}^2}{A}}$.

Its constant solutions are $Y=i{\displaystyle \frac{\tilde{\kappa }}{A}}$ (both signs are allowed). According to equation (18) ${\displaystyle \frac{y_2^{\prime }}{y_2}}=i\tilde{\kappa }$, i.e., $\tilde{\kappa }$ is the $\widehat{x}$ wavenumber. The solutions corresponds to the continuous spectrum of the slow/fast magnetosonic waves, and as already stated, the equation ${\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-{\mathit{k}}_0^2$ is equivalent to the dispersion relation of these waves ${\displaystyle \frac{{\omega }_0^4}{{\mathit{k}}_0^2+{\tilde{\kappa }}^2}}-{\omega }_0^2\left({\displaystyle \frac{B_0^2}{{\rho }_0}}+c_s^2\right)+c_s^2{\displaystyle \frac{F^2}{{\rho }_0}}=0$.

The variable solution of the principal equation $Y={\displaystyle \frac{\tilde{\kappa }}{A}}cot(\tilde{\kappa }x+{\varphi }_0)$ corresponds to a superposition of waves as discussed in Appendix B.

3.2. Alfvèn Waves

The continuous spectrum of Alfvèn waves ${\omega }_0^2=F^2/{\rho }_0$ corresponds to the singular case of the principal equation $A=0$. In that case $y_1$ could be any function of x, $y_2=0$ (and $Y\to \infty$), with finite ${\displaystyle \frac{y_2}{A}}=-{\displaystyle \frac{y_1^{\prime }}{{\tilde{\kappa }}^2}}$ and ${\displaystyle \frac{1}{AY}}=-{\displaystyle \frac{y_1^{\prime }}{{\tilde{\kappa }}^2{y}_1}}$. (Note that for the Alfvén waves $\tilde{\kappa }$ is given by ${\tilde{\kappa }}^2={({\mathit{k}}_0·{\mathit{B}}_0/B_0)}^2-{\mathit{k}}_0^2=-{({\mathit{k}}_0\times {\mathit{B}}_0/B_0)}^2$, and does not represent a wavenumber.)

3.3. Two Semi-Infinite Incompressible Plasmas

As another example on how to use the minimalist approach lets consider two semi-infinite incompressible, homogeneous, static, magnetized plasmas in the regimes $x>0$ (subscript "1") and $x<0$ (subscript "2"), and zero gravity.

With $S\to \infty$ in both regimes we have $f_{11}=0$, $f_{12}=-{\displaystyle \frac{k_0^2}{A}}$, $f_{21}=-A$. The principal equation is ${\displaystyle \frac{dY}{dx}}={\displaystyle \frac{k_0^2}{A}}-A{Y}^2$ and has constant solutions $Y={\displaystyle \frac{f_{11}\pm \sqrt{f_{11}^2+f_{12}{f}_{21}}}{f_{21}}}=\mp {\displaystyle \frac{k_0}{A}}$. According to the boundary conditions given in Table 1 (or directly checking the sign of ${\displaystyle \frac{y_2^{\prime }}{y_2}}=-f_{21}Y+f_{11}=\mp k_0$) the upper sign should be used for $x>0$ and the lower for $x<0$. Thus $Y=-{\displaystyle \frac{k_0}{{\rho }_1{\omega }^2-F_1^2}}$ for $x>0$ and $Y=+{\displaystyle \frac{k_0}{{\rho }_2{\omega }^2-F_2^2}}$ for $x<0$.

The continuity of Y gives the dispersion relation ${\omega }^2={\displaystyle \frac{F_1^2+F_2^2}{{\rho }_1+{\rho }_2}}$ representing a stable Alfvén wave whose amplitude drops exponentially as we move away from the contact discontinuity at $x=0$. (We remind that $F_1={\mathit{k}}_0·{\mathit{B}}_{01}$, $F_2={\mathit{k}}_0·{\mathit{B}}_{02}$, and that all the unperturbed states considered in the paper have the form described in Section 2.)

3.4. Effect of Finite Depth

Following the previous example, suppose now that the bottom plasma has finite depth H (there is a solid wall at $x=-H$). Then in the bottom part the solution of the principal equation ${\displaystyle \frac{dY}{dx}}={\displaystyle \frac{k_0^2}{A_2}}-A_2{Y}^2$ that vanishes at $x=-H$ is $Y={\displaystyle \frac{k_0}{A_2}}tanh\left[k_0(x+H)\right]$. (In the upper part the accepted solution corresponding to vanishing amplitude at $x\to +\infty$ continues to be $Y=-{\displaystyle \frac{k_0}{A_1}}$ as before.)

The continuity of Y at $x=0$ gives the dispersion relation ${\omega }^2={\displaystyle \frac{F_1^2+F_2^{2}coth\left(k_{0}H\right)}{{\rho }_1+{\rho }_{2}coth\left(k_{0}H\right)}}$.

4. Gravity Waves and Rayleigh-Taylor Instability

Consider two semi-infinite incompressible, static, magnetized plasmas, the first (subscript "1") in the region $x>0$ and the second (subscript "2") in the region $x<0$, inside uniform gravity $-g\widehat{x}$, as in the left panel of Figure 1.

With $S\to \infty$ in both regimes we have $f_{11}=-g{\rho }_0{\displaystyle \frac{k_0^2}{A}}$, $f_{12}=-{\displaystyle \frac{k_0^2}{A}}$, $f_{21}=-A+g^2{\rho }_0^2{\displaystyle \frac{k_0^2}{A}}$. The principal equation is ${\displaystyle \frac{dY}{dx}}={\displaystyle \frac{k_0^2{(1+g{\rho }_{0}Y)}^2}{A}}-A{Y}^2$ and has constant solutions $Y={\displaystyle \frac{f_{12}}{-f_{11}\pm \sqrt{f_{11}^2+f_{12}{f}_{21}}}}={\displaystyle \frac{-k_0}{{\rho }_{0}g{k}_0\pm A}}$. (Although these parts are not homogeneous, since there must be a gradient of the total pressure balancing gravity, the incompressibility assumption allows to have constant solutions, provided that ${\rho }_0$ and A are constant in both parts.)

Choosing the correct sign in each part following Table 1 (or directly using ${\displaystyle \frac{y_2^{\prime }}{y_2}}=-f_{21}Y+f_{11}=\mp k_0$) we conclude that $Y={\displaystyle \frac{-k_0}{{\rho }_1{\omega }^2-F_1^2+{\rho }_{1}g{k}_0}}$ for $x>0$ and $Y={\displaystyle \frac{-k_0}{-{\rho }_2{\omega }^2+F_2^2+{\rho }_{2}g{k}_0}}$ for $x<0$.

The continuity of Y gives the dispersion relation ${\omega }^2=-{\displaystyle \frac{{\rho }_1-{\rho }_2}{{\rho }_1+{\rho }_2}}g{k}_0+{\displaystyle \frac{F_1^2+F_2^2}{{\rho }_1+{\rho }_2}}$.

This relation includes as subcases the surface gravity waves corresponding to zero magnetic fields and ${\rho }_1\ll {\rho }_2$, for which ${\omega }^2=g{k}_0$, and the classical Rayleigh-Taylor instability corresponding to zero magnetic fields and heavier fluid on top of lighter, for which the growth rate is $\omega =i\sqrt{{\displaystyle \frac{{\rho }_1-{\rho }_2}{{\rho }_1+{\rho }_2}}g{k}_0}$. It also shows the stabilizing effect of the magnetic field through its magnetic tension manifest in the dispersion relation in the presence of $F_1^2$ and $F_2^2$. It always reduces the growth rate, and if it is sufficiently strong it makes ${\omega }^2$ positive, converting the disturbance from an instability to an Alfvén wave.

Effect of Finite Depth

If the bottom plasma has finite depth H (there is a solid wall at $x=-H$, see the right panel of Figure 1) we need to consider the variable solution of the principal equation $Y={\displaystyle \frac{f_{12}}{-f_{11}-Kcot(Kx+{\varphi }_0)}}$ with $K=i\sqrt{f_{11}^2+f_{12}{f}_{21}}=i{k}_0$ (given in Appendix B). Requiring Y to vanish at $x=-H$ we find $Y={\displaystyle \frac{-k_0}{{\rho }_{2}g{k}_0-A_{2}coth(k_{0}x+k_{0}H)}}$ for $x<0$.

As a specific example suppose the upper part is a hydrodynamic medium with negligible density. In that case $Y=\infty$ for $x>0$ and the continuity of Y gives the dispersion relation for gravity waves ${\omega }^2={\displaystyle \frac{F_2^2}{{\rho }_2}}+g{k}_{0}tanh\left(k_{0}H\right)$, modified by the presence of the magnetic field. (The same can be found considering the constant solution $Y={\displaystyle \frac{-k_0}{{\rho }_1{\omega }^2+{\rho }_{1}g{k}_0}}$ in the upper part, requiring the continuity of Y at $x=0$, and then taking the limit ${\rho }_1/{\rho }_2\to 0$.) In the shallow water approximation $k_{0}H\ll 1$ we get ${\omega }^2={\displaystyle \frac{F^2}{{\rho }_2}}+gH{k}_0^2$.

5. Gravito-Acoustic Waves

In the examples presented up to now the fluids where homogeneous at least partially, allowing to easily solve the principal equation. In inhomogeneous cases we rely on the numerical integration of the principal equation. An exception to that, where it is possible to solve analytically the differential equation although the fluid is not homogeneous, is the case of gravito-acoustic waves. In addition it serves as a nice example showing the interplay between acoustic and gravity waves.

For a hydrodynamic fluid between two solid walls at $x=0$ and $x=H$, inside homogeneous gravity $\mathit{g}=-g\widehat{x}$, with unperturbed density ${\rho }_0={\rho }_{00}{e}^{-\alpha x}$, pressure $P_0={\displaystyle \frac{{\rho }_{0}g}{\alpha }}$ satisfying the hydrostatic equilibrium condition $P_0^{\prime }=-{\rho }_{0}g$, and ${\displaystyle \frac{1}{\alpha }}={\displaystyle \frac{c_s^2}{\gamma g}}$ the constant scaleheight, we get $f_{11}=-{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}$, $f_{12}={\displaystyle \frac{1/c_s^2-k_0^2/{\omega }^2}{{\rho }_{00}}}{e}^{\alpha x}$, $f_{21}={\rho }_{00}(g^2{k}_0^2/{\omega }^2-{\omega }^2)e^{-\alpha x}$, and the principal equation can be written as

$${\displaystyle \frac{d\left(Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}\right)}{dx}}=\left({\displaystyle \frac{g^2{k}_0^2}{{\omega }^4}}-1\right){\left(Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}\right)}^2+\left(2{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}-\alpha \right)\left(Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}\right)+k_0^2-{\displaystyle \frac{{\omega }^2}{c_s^2}}.$$

Its general solution is $\left({\displaystyle \frac{g^2{k}_0^2}{{\omega }^4}}-1\right)Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}={\displaystyle \frac{\alpha }{2}}-{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}+qtan(qx+C)$ with $q=\sqrt{-{\displaystyle \frac{{\alpha }^2}{4}}+{\displaystyle \frac{{\omega }^2}{c_s^2}}-k_0^2+{\displaystyle \frac{k_0^2{N}^2}{{\omega }^2}}}$, and $N^2=g\left({\displaystyle \frac{1}{\gamma P}}{\displaystyle \frac{dP}{dx}}-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d\rho }{dx}}\right)=\alpha g-{\displaystyle \frac{g^2}{c_s^2}}={\displaystyle \frac{\gamma -1}{\gamma }}\alpha g$ the square of the Brunt-Väisälä (or buoyancy) frequency. The boundary condition ${Y|}_{x=0}=0$ specifies C and gives $\left({\displaystyle \frac{g^2{k}_0^2}{{\omega }^4}}-1\right)Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}={\displaystyle \frac{q^2+{\left({\displaystyle \frac{\alpha }{2}}-{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}\right)}^2}{qcot\left(qx\right)+{\displaystyle \frac{\alpha }{2}}-{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}}}={\displaystyle \frac{\left({\displaystyle \frac{g^2{k}_0^2}{{\omega }^4}}-1\right)\left(k_0^2-{\displaystyle \frac{{\omega }^2}{c_s^2}}\right)}{qcot\left(qx\right)+{\displaystyle \frac{\alpha }{2}}-{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}}}$. The boundary condition ${Y|}_{x=H}=0$ requires $qH=n\pi$, $n=\pm 1,\pm 2,\dots$ and the dispersion relation is $-{\displaystyle \frac{{\alpha }^2}{4}}+{\displaystyle \frac{{\omega }^2}{c_s^2}}-k_0^2+{\displaystyle \frac{k_0^2{N}^2}{{\omega }^2}}=q^2\iff {\displaystyle \frac{{\omega }^4}{c_s^2}}-\left(k_0^2+q^2+{\displaystyle \frac{{\alpha }^2}{4}}\right){\omega }^2+k_0^2{N}^2=0$.

The equation for $y_1$ is ${\displaystyle \frac{y_1^{\prime }}{y_1}}={\displaystyle \frac{g{k}_0^2}{{\omega }^2}}+{\displaystyle \frac{k_0^2-{\omega }^2/c_s^2}{Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}}}={\displaystyle \frac{\alpha }{2}}+qcot\left(qx\right)$ with solution $y_1\propto e^{\alpha x/2}sin\left(qx\right)$.

(The cases ${\omega }^2=k_0^2{c}_s^2$ and ${\omega }^4=k_0^2{g}^2$ need to be examined separately. The former gives $q^2=-{(\alpha /2-g/c_s^2)}^2$ and the latter gives $\left(2{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}-\alpha \right)\left(Y{\rho }_{00}{\omega }^2{e}^{-\alpha x}\right)+k_0^2-{\omega }^2/c_s^2\propto exp\left[\left(2{\displaystyle \frac{g{k}_0^2}{{\omega }^2}}-\alpha \right)x\right]$. Neither can satisfy the boundary conditions, so they cannot be accepted).

6. Kelvin-Helmholtz Instability

6.1. Incompressible Limit

For the incompressible case and two semi-infinite plasmas, similarly to the Rayleigh-Taylor instability (Section 4) we find that $Y={\displaystyle \frac{-k_0}{{\rho }_{1}g{k}_0+A_1}}$ for $x>0$ and $Y={\displaystyle \frac{-k_0}{{\rho }_{2}g{k}_0-A_2}}$ for $x<0$. The only difference is that $A_1={\rho }_1{(\omega -{\mathit{k}}_0·{\mathit{V}}_{01})}^2-F_1^2$ and $A_2={\rho }_2{(\omega -{\mathit{k}}_0·{\mathit{V}}_{02})}^2-F_2^2$ now depend on the velocities.

The continuity of Y gives the dispersion relation ${\rho }_1{(\omega -{\mathit{k}}_0·{\mathit{V}}_{01})}^2+{\rho }_2{(\omega -{\mathit{k}}_0·{\mathit{V}}_{02})}^2=F_1^2+F_2^2-({\rho }_1-{\rho }_2)g{k}_0$ with solutions

$$\begin{array}{c}\hfill \omega ={\displaystyle \frac{{\rho }_1{\mathit{V}}_{01}+{\rho }_2{\mathit{V}}_{02}}{{\rho }_1+{\rho }_2}}·{\mathit{k}}_0\pm i\sqrt{{\displaystyle \frac{{\rho }_1{\rho }_2{\left[({\mathit{V}}_{01}-{\mathit{V}}_{02})·{\mathit{k}}_0\right]}^2}{{({\rho }_1+{\rho }_2)}^2}}+{\displaystyle \frac{{\rho }_1-{\rho }_2}{{\rho }_1+{\rho }_2}}g{k}_0-{\displaystyle \frac{{({\mathit{k}}_0·{\mathit{B}}_{01})}^2+{({\mathit{k}}_0·{\mathit{B}}_{02})}^2}{{\rho }_1+{\rho }_2}}}.\end{array}$$

(19)

The three terms within the square root have obvious physical meaning. The first is due to relative velocity and always leads to instability. The second is due to gravity and leads to stability/instability in the case of stable/unstable stratification (lighter fluid above heavier or the opposite). The third term is due to the magnetic field stress and is always stabilizing.

6.2. Compressible Kelvin-Helmholtz Instability

Taking into account compressibility and neglecting gravity we have $f_{11}=0$, $f_{12}={\displaystyle \frac{{\tilde{\kappa }}^2}{A}}$, $f_{21}=-A$ and the constant solutions of the principal equation are $Y={\displaystyle \frac{f_{11}\pm \sqrt{f_{11}^2+f_{12}{f}_{21}}}{f_{21}}}=i{\displaystyle \frac{\tilde{\kappa }}{A}}$ with $\tilde{\kappa }=\pm i\sqrt{-{\tilde{\kappa }}^2}$. According to Table 1 we need to choose the upper sign in the $x>0$ and the lower in the $x<0$. Thus ${\tilde{\kappa }}_1=i\sqrt{-{\tilde{\kappa }}_1^2}$ and $Y=i{\displaystyle \frac{{\tilde{\kappa }}_1}{A_1}}={\displaystyle \frac{\sqrt{-{\tilde{\kappa }}_1^2}}{-A_1}}$ for $x>0$, ${\tilde{\kappa }}_2=-i\sqrt{-{\tilde{\kappa }}_2^2}$ and $Y=i{\displaystyle \frac{{\tilde{\kappa }}_2}{A_2}}={\displaystyle \frac{\sqrt{-{\tilde{\kappa }}_2^2}}{A_2}}$ for $x<0$. The dispersion relation is ${\displaystyle \frac{{\tilde{\kappa }}_1}{A_1}}={\displaystyle \frac{{\tilde{\kappa }}_2}{A_2}}$, or,

$$\begin{array}{c}\hfill {\displaystyle \frac{\sqrt{{\mathit{k}}_0^2-{\displaystyle \frac{{\rho }_1{(\omega -{\mathit{k}}_0·{\mathit{V}}_{01})}^4}{({\rho }_1{c}_{s1}^2+B_{01}^2){(\omega -{\mathit{k}}_0·{\mathit{V}}_{01})}^2-{({\mathit{k}}_0·{\mathit{B}}_{01})}^2{c}_{s1}^2}}}}{-{\rho }_1{(\omega -{\mathit{k}}_0·{\mathit{V}}_{01})}^2+{({\mathit{k}}_0·{\mathit{B}}_{01})}^2}}={\displaystyle \frac{\sqrt{{\mathit{k}}_0^2-{\displaystyle \frac{{\rho }_2{(\omega -{\mathit{k}}_0·{\mathit{V}}_{02})}^4}{({\rho }_2{c}_{s2}^2+B_{02}^2){(\omega -{\mathit{k}}_0·{\mathit{V}}_{02})}^2-{({\mathit{k}}_0·{\mathit{B}}_{02})}^2{c}_{s2}^2}}}}{{\rho }_2{(\omega -{\mathit{k}}_0·{\mathit{V}}_{02})}^2-{({\mathit{k}}_0·{\mathit{B}}_{02})}^2}}.\end{array}$$

(20)

(We remind that we need to take the principal values of both square roots, i.e., the ones with positive real parts.)

6.3. On the Physics of the Kelvin-Helmholtz Instability

A relatively simple case for which we can find analytical results will help to understand the physics of the instability. We consider two hydrodynamic fluids with the same unperturbed characteristics (same densities ${\rho }_0$, same pressure $P_0$ as required from the equilibrium at the contact discontinuity, and same sound speeds $c_s$) and work in the frame where the two fluids move with opposite velocities $\pm V_0\widehat{z}$ (with positive $V_0$ and the upper sign corresponding to $x>0$).

The dispersion relation becomes ${\displaystyle \frac{\sqrt{1-M^2{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}-1\right)}^2}}{-{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}-1\right)}^2}}={\displaystyle \frac{\sqrt{1-M^2{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}+1\right)}^2}}{{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}+1\right)}^2}}$, where $M={\displaystyle \frac{V_0{k}_{0z}/k_0}{c_s}}$. Unstable modes (with $\Im \omega >0$) exist for $M<\sqrt{2}$, with purely imaginary $\omega =i{k}_{0z}{V}_0\sqrt{{\displaystyle \frac{2-M^2}{1+M^2+\sqrt{1+4M^2}}}}$. One way to find this result is to substitute $\omega =i{k}_{0z}{V}_{0}tan{\displaystyle \frac{\mu }{2}}$ when the dispersion relation reduces to ${cos}^2\mu +cos\mu =M^2$ with solutions $cos\mu ={\displaystyle \frac{\sqrt{1+4M^2}-1}{2}}$. Actually the angle $\mu$ has an important meaning connected to the argument of the resulting complex wavenumbers in the $\widehat{x}$ direction, which are ${\tilde{\kappa }}_{1,2}=-k_0{e}^{\mp i\mu }=-k_0{\displaystyle \frac{\sqrt{1+4M^2}-1}{2}}\pm i{k}_0\sqrt{{\displaystyle \frac{\sqrt{1+4M^2}+1-2M^2}{2}}}$.

Although the results apply for any angle between ${\mathit{k}}_0$ and ${\mathit{V}}_0$, to simplify the expressions from now on we restrict ourselves to the case ${\mathit{k}}_0\Vert {\mathit{V}}_0$, i.e., we consider a disturbance with ${\mathit{k}}_0=k_0\widehat{z}$ (with positive $k_0$). The characteristics of the unstable mode, the Lagrangian displacement, and the perturbations of the pressure, density, and velocity (written in a way to show the phase difference between each quantity with $y_1$) are

$$\begin{array}{c}\hfill \omega =i{k}_0{V}_{0}tan{\displaystyle \frac{\mu }{2}},{\tilde{\kappa }}_{1,2}=-k_0{e}^{\mp i\mu },\mu =arccos{\displaystyle \frac{\sqrt{1+4M^2}-1}{2}},M={\displaystyle \frac{V_0}{c_s}},\end{array}$$

(21)

$$\begin{array}{c}\hfill y_1{e}^{i(k_{z}z-\omega t)}\propto e^{k_0{V}_{0}ttan(\mu /2)}{e}^{\pm k_{0}xsin\mu }{e}^{i{k}_0(z-xcos\mu )},\end{array}$$

(22)

$$\begin{array}{c}\hfill P_1=y_2={\displaystyle \frac{{\rho }_0{k}_0{V}_0^2}{{cos}^2(\mu /2)}}{e}^{i\pi /2}{y}_1,{\rho }_1={\displaystyle \frac{P_1}{c_s^2}},\end{array}$$

(23)

$$\begin{array}{c}\hfill V_{1x}={\displaystyle \frac{k_0{V}_0}{cos(\mu /2)}}{e}^{\pm i(\pi /2-\mu /2)}{y}_1,V_{1z}={\displaystyle \frac{k_0{V}_0}{cos(\mu /2)}}{e}^{\mp i(\pi /2-\mu /2)}{y}_1.\end{array}$$

(24)

Suppose we initially disturb the interface between the two fluids by ${\Delta x|}_{x=0,t=0}=D{e}^{i{k}_{0}z}$. The Lagrangian displacement along $\widehat{x}$ at any later time and for any other fluid element has the form ${\xi }_x=y_1{e}^{i(k_{0}z-\omega t)}=D{e}^{i(k_{0}z+{\tilde{\kappa }}_{1,2}x-\omega t)}$; the goal is to find $\omega$ and ${\tilde{\kappa }}_{1,2}$. Using the Lagrangian displacement $\mathit{\xi }$ and the Lagrangian perturbations of the velocity, pressure, density, which are $\Delta \mathit{V}=\dot{\mathit{\xi }}$, $\Delta P$, $\Delta \rho ={\displaystyle \frac{\Delta P}{c_s^2}}$, respectively, the equations that will lead to that goal are the momentum ${\rho }_0\ddot{\mathit{\xi }}=-\nabla (\Delta P)$ and the continuity $\nabla ·\mathit{\xi }=-{\displaystyle \frac{\Delta \rho }{{\rho }_0}}$ (coming from $\nabla ·\mathit{V}=-{\displaystyle \frac{\dot{\rho }}{\rho }}$). Since all perturbations are proportional to $e^{i(k_{0}z+\tilde{\kappa }x-\omega t)}$ a time derivative has the meaning ${\displaystyle \frac{\partial }{\partial t}}+{\mathit{V}}_0·\nabla =-i(\omega \mp k_0{V}_0)$, and the divergence $\nabla =i{\mathit{k}}_{t1,2}$ with ${\mathit{k}}_{t1,2}={\tilde{\kappa }}_{1,2}\widehat{x}+{\mathit{k}}_0$ the total wavenumber in the two fluids (their $\widehat{z}$ components are the same, but the $\widehat{x}$ components differ).

The momentum equation along $\widehat{x}$ connects the displacement with the pressure perturbation $\Delta P={\displaystyle \frac{{\rho }_0{(\omega \mp k_0{V}_0)}^2}{i{\tilde{\kappa }}_{1,2}}}{y}_1$. Requiring $\Delta P$ and $y_1$ to be continuous at the interface between the two fluids we arrive at a first relation between the unknowns ${\displaystyle \frac{{\tilde{\kappa }}_1}{A_1}}={\displaystyle \frac{{\tilde{\kappa }}_2}{A_2}}$.

The pressure gradient along the interface is connected to the corresponding Lagrangian displacement ${\xi }_z$ through the momentum equation along $\widehat{z}$, which gives ${\xi }_z={\displaystyle \frac{k_0}{\tilde{\kappa }}}{y}_1$. (The latter means that $\mathit{\xi }\Vert {\mathit{k}}_t$, a characteristic of longitudinal waves.)

The continuity also connects the two components of the displacement. It gives $k_0{\xi }_z+\tilde{\kappa }{y}_1=i{\displaystyle \frac{\Delta P}{{\rho }_0{c}_s^2}}$. Substituting $\Delta P={\displaystyle \frac{{\rho }_0{(\omega \mp k_0{V}_0)}^2}{i{\tilde{\kappa }}_{1,2}}}{y}_1$ and ${\xi }_z={\displaystyle \frac{k_0}{\tilde{\kappa }}}{y}_1$ we arrive at the second relation between the unknowns (one for each fluid) ${\tilde{\kappa }}_{1,2}^2={\displaystyle \frac{{(\omega \mp k_0{V}_0)}^2}{c_s^2}}-k_0^2$.

The perturbation essentially consists of two sound waves in the two fluids with wavenumbers ${\mathit{k}}_{t1,2}={\tilde{\kappa }}_{1,2}\widehat{x}+{\mathit{k}}_0$. In the fluid rest frames the two waves move with different velocities satisfying ${(\omega -{\mathit{k}}_0·{\mathit{V}}_0)}^2=c_s^2{\tilde{\kappa }}_{1,2}^2$, a relation equivalent to ${\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-k_0^2$ with $S={\rho }_0^2{c}_s^2{\omega }_0^2$. The two waves meet at the contact discontinuity and as a consequence they have the same ${\mathit{k}}_0$ and $\omega$ such that their phases ${\mathit{k}}_t·(\mathit{r}-{\mathit{V}}_{0}t)-(\omega -{\mathit{k}}_0·{\mathit{V}}_0)t={\mathit{k}}_t·\mathit{r}-\omega t$ are continuous at the interface. As already presented, for the case of fluids with the same density, in a frame in which they move with opposite velocities $\pm V_0\widehat{z}$, and ${\mathit{k}}_0=k_0\widehat{z}$ is parallel to the relative velocity, the resulting unstable mode has $\omega =i{k}_0{V}_{0}tan{\displaystyle \frac{\mu }{2}}$ and ${\tilde{\kappa }}_{1,2}=-k_0{e}^{\mp i\mu }$, with $\mu =arccos{\displaystyle \frac{\sqrt{1+4M^2}-1}{2}}$, $M={\displaystyle \frac{V_0}{c_s}}$.

The value of the complex wavenumber $\tilde{\kappa }$ (in each part) is directly connected to compressibility through S. In the incompressible limit $S\to \infty$ it is purely imaginary, since in that limit $M=0$, $\mu =\pi /2$ and $\tilde{\kappa }=\pm i{k}_0$. (Note that in the incompressible limit the continuity $\nabla ·\mathit{\xi }=0$ gives ${\mathit{k}}_t·\mathit{\xi }=0$. Since the two vectors are complex, this does not mean ${\mathit{k}}_t\perp \mathit{\xi }$. It rather means ${\mathit{k}}_t^2=0$, using ${\mathit{k}}_t\Vert \mathit{\xi }$ from the momentum equation. In addition, the relation ${\mathit{k}}_t^2=0$ does not mean that the vector ${\mathit{k}}_t$ is zero.) If we consider cases with decreasing S, i.e., decreasing $c_s$ keeping $V_0$ the same, in which the compressibility becomes more and more important, the first part of the expression of ${\tilde{\kappa }}^2={\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{S}}-k_0^2$ affects the value of $\tilde{\kappa }$ in two ways. Firstly the imaginary part decreases, meaning that the perturbation survives at longer distances from the interface. Secondly the $\Re \tilde{\kappa }$ increases contributing more to the phase of the perturbation which is $k_0(z-xcos\mu )$. Both effects are expressed through the angle $\mu$, which decreases with decreasing S (increasing M). The growth rate is also connected to S through $\mu$. It decreases as the result of compressibility, from $k_0{V}_0$ in the limit $S\to \infty$ to zero when $M=\sqrt{2}$, corresponding to the minimum S for which the perturbation is unstable.

The vorticity is a key quantity in the Kelvin-Helmholtz instability. Its undisturbed value ${\mathit{\zeta }}_0=-2V_0\delta \left(x\right)\widehat{y}$ is the reason behind the development of the instability. Even with the perturbation included it is nonzero only in the interface between the two fluids. Nevertheless, the motion of the fluids along the interface redistributes the vorticity compared to the unperturbed state. The related velocity inside the upper/lower fluid, just above/below the interface ${\xi }_x{{|}_{x=0}=y_1|}_{x=0}{e}^{i(k_{z}z-\omega t)}$, is $V_{1z}{{|}_{x=0^{\pm }}=k_0{V}_0[tan(\mu /2)\mp i]{\xi }_x|}_{x=0}$. The mean value $k_0{V}_{0}tan(\mu /2){\xi }_x{|}_{x=0}$ shows that fluid accumulates near the positions of lower pressure, where the vorticity increases (in the $-\widehat{y}$ direction) by $-i{k}_0{\xi }_x{|}_{x=0}{\mathit{\zeta }}_0$. This accumulation of vorticity further increases the displacement ${\xi }_x{|}_{x=0}$ leading to instability. It becomes stronger for larger $\mu$, i.e., smaller M, since the mean value of $V_{1z}$ increases with $\mu$.

Three example solutions are shown in Figure 2. The upper panel corresponds to an incompressible case with $M\ll 1$ for which $\omega \approx i(1-M^2)k_0{V}_0$ and ${\tilde{\kappa }}_{1,2}\approx (-M^2\pm i)k_0$. (For $M\ll 1$ the approximate expression for $\mu$ is $\mu \approx {\displaystyle \frac{\pi }{2}}-M^2$.) The exponential decrease of the perturbation as we move away from the interface is evident. The perturbed interface is shown, as also the circulation around points of minimum pressure (the mean velocity of the fluids on the interface points toward these minimum pressure points).

The two other panels correspond to cases in which compressibility is important. The solution in the middle panel has $M=1$, $\omega =i0.486k_0{V}_0$ and ${\tilde{\kappa }}_{1,2}=(-0.618\pm i0.786)k_0$ and in the lower panel $M=1.4$, $\omega =i0.082k_0{V}_0$ and ${\tilde{\kappa }}_{1,2}=(-0.987\pm i0.163)k_0$. The $\Re {\mathit{k}}_t=\Re \tilde{\kappa }\widehat{x}+k_0\widehat{z}$ has now a nonnegligible $\widehat{x}$ part and as a result the iso-phase planes are $z-xcos\mu =$ constant, tilted as shown in the panels.

For cases approaching the maximum value of M for which the perturbation is unstable (the $M=\sqrt{2}$) the values of $\mu$, $\omega$ and $\Im \tilde{\kappa }$ approach zero, and the perturbation practically consists of two standing sound waves (one in each fluid) with real wavenumbers. (For $M^2-2\to 0^{-}$ the approximate expression for $\mu$ is $\mu \approx \sqrt{{\displaystyle \frac{4-2M^2}{3}}}$.)

6.4. Influence of Magnetic Field

For simplicity we consider again two homogeneous fluids with the same unperturbed characteristics, and work in the frame where the fluids move with opposite velocities $\pm V_0\widehat{z}$. Now there is also a constant magnetic field ${\mathit{B}}_0=B_{0y}\widehat{y}+B_{0z}\widehat{z}$ in the unperturbed state, the same in both fluids.

Inspecting the dispersion relation (20) we see that the magnetic field enters in two ways. Firstly through its pressure $B_0^2/2$, or equivalently the square of the Alfvén velocity ${\mathit{v}}_{\mathrm{A}}={\displaystyle \frac{{\mathit{B}}_0}{\sqrt{{\rho }_0}}}$. Secondly through its tension manifest in the terms ${\mathit{k}}_0·{\mathit{B}}_0$, i.e. its component parallel to ${\mathit{k}}_0$, or equivalently the component of the Alfvén velocity $v_{\mathrm{A}\Vert }={\displaystyle \frac{{\mathit{B}}_0·{\mathit{k}}_0/k_0}{\sqrt{{\rho }_0}}}$.

The magnetic pressure $B_0^2/2$ enters in the expression of S and increases its value. Thus it moves the dynamics toward the incompressible limit, and according to the discussion of the previous section destabilizes (the second terms inside the square roots equal ${\displaystyle \frac{{\rho }_0^2{\omega }_0^4}{k_0^{2}S}}$, so increase of $v_{\mathrm{A}}^2$ leads to the square roots being closer to unity). Actually if we include magnetic filed normal to ${\mathit{k}}_0$ only, the dispersion relation is exactly equivalent to its hydrodynamic analogue, with the only difference that M represents now the fast magnetosonic Mach number $M={\displaystyle \frac{V_0{k}_{0z}/k_0}{\sqrt{c_s^2+v_{\mathrm{A}}^2}}}$. Figure 3 shows the resulting growth rates. Notably it includes as subcases the purely hydrodynamic case (${\mathit{v}}_A=0$) considered in the previous section, and the cold case ($c_s=0$) with only $B_{0y}$ component of the magnetic field.

The magnetic tension enters in the dispersion relation through $F={\mathit{k}}_0·{\mathit{B}}_0$ (essentially the component of the magnetic field along ${\mathit{k}}_0$) in two places: inside A and S. It always has a stabilizing effect since the tension is a restoring force; we have seen that in the previous examples of Alfvén waves and the Rayleigh-Taylor instability, even in the incompressible limit. It also affects incompressibility through its appearance inside S. Since it decreases S it moves the dynamics away from the incompressible limit, something that also in general stabilizes.

Another related connection is that the magnetic tension affects the perturbation of the vorticity inside each fluid, which, using the expressions of the velocity perturbations given in Appendix A, can be expressed as $\Delta \mathit{\zeta }=-i{\displaystyle \frac{{\omega }_{0}F{\mathit{k}}_t\times {\mathit{B}}_0}{{\rho }_{0}A}}\Delta \rho$.

In the following we present the methodology to find numerical results.

The dispersion relation (20) in the case under consideration can be written as

$$\begin{array}{c}\hfill {\displaystyle \frac{\sqrt{1-{\displaystyle \frac{{\displaystyle \frac{V_0^2{k}_{0z}^2/k_0^2}{c_s^2+v_{\mathrm{A}}^2}}{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}-1\right)}^4}{{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}-1\right)}^2-{\displaystyle \frac{c_s^2}{c_s^2+v_{\mathrm{A}}^2}}{\displaystyle \frac{v_{\mathrm{A}\Vert }^2}{V_0^2}}}}}}{-{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}-1\right)}^2+{\displaystyle \frac{v_{\mathrm{A}\Vert }^2}{V_0^2}}}}={\displaystyle \frac{\sqrt{1-{\displaystyle \frac{{\displaystyle \frac{V_0^2{k}_{0z}^2/k_0^2}{c_s^2+v_{\mathrm{A}}^2}}{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}+1\right)}^4}{{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}+1\right)}^2-{\displaystyle \frac{c_s^2}{c_s^2+v_{\mathrm{A}}^2}}{\displaystyle \frac{v_{\mathrm{A}\Vert }^2}{V_0^2}}}}}}{{\left({\displaystyle \frac{\omega }{k_{0z}{V}_0}}+1\right)}^2-{\displaystyle \frac{v_{\mathrm{A}\Vert }^2}{V_0^2}}}},\end{array}$$

(25)

and can be solved numerically. However, for the case of fluids with the same characteristics considered here it is possible to proceed analytically. As shown in Appendix C the dispersion relation can be transformed to the following quartic polynomial equation for the square of the growth rate ${\displaystyle \frac{{(\Im \omega )}^2}{k_0^2(c_s^2+v_{\mathrm{A}}^2)}}=M^2{tan}^2{\displaystyle \frac{\lambda }{2}}$

$$\begin{array}{c}\hfill M^8{tan}^8{\displaystyle \frac{\lambda }{2}}+2(1+2M^2)M^6{tan}^6{\displaystyle \frac{\lambda }{2}}\\ \hfill +[2M^2(1+3M^2)+2(2{cos}^2\eta +1){cos}^2\mathsf{\Theta }-(2{cos}^2\eta +1){cos}^4\mathsf{\Theta }]M^4{tan}^4{\displaystyle \frac{\lambda }{2}}\\ \hfill +\left\{2M^4(2M^2-1)-4M^2(2{cos}^2\eta -1){cos}^2\mathsf{\Theta }\right.\\ \hfill \left.+2\left[{cos}^4\eta -2(M^2-1){cos}^2\eta -M^2\right]{cos}^4\mathsf{\Theta }-2{cos}^2\eta {cos}^6\mathsf{\Theta }\right\}{M}^2{tan}^2{\displaystyle \frac{\lambda }{2}}\\ \hfill +M^6(M^2-2)+2M^4(2{cos}^2\eta +1){cos}^2\mathsf{\Theta }\\ \hfill -M^2[2{cos}^4\eta +2(M^2+2){cos}^2\eta +M^2]{cos}^4\mathsf{\Theta }+2(M^2+{cos}^2\eta ){cos}^2\eta {cos}^6\mathsf{\Theta }=0,\end{array}$$

(26)

using the parametrization

$$\begin{array}{c}\hfill \omega =i{k}_{0z}{V}_{0}tan{\displaystyle \frac{\lambda }{2}},M={\displaystyle \frac{V_0{k}_{0z}/k_0}{\sqrt{c_s^2+v_{\mathrm{A}}^2}}},{\displaystyle \frac{v_{\mathrm{A}}}{c_s}}=tan\eta ,{\displaystyle \frac{v_{\mathrm{A}\Vert }}{V_0}}={\displaystyle \frac{cos\mathsf{\Theta }}{M}}.\end{array}$$

(27)

(Note that the angle $\mathsf{\Theta }$ is connected to the angle between the magnetic field and the wavenumber since $cos\mathsf{\Theta }={\displaystyle \frac{{\mathit{B}}_0·{\mathit{k}}_0}{B_0{k}_0}}{\displaystyle \frac{k_{0z}}{k_0}}sin\eta$.)

The left panel of Figure 4 shows the growth rate for cases where the magnetic field has only component along ${\mathit{k}}_0$. In general the field decreases the growth rate and if it is sufficiently strong completely suppresses the instability.

Similar behavior is shown in the right panel of Figure 4 for the cold case. For all strengths of the magnetic field, if its orientation is sufficiently close to the wavevector the magnetic tension completely suppresses the instability.

6.5. Range of Instability

It is interesting to explore the regions of M for which the Kelvin-Helmholtz instability is present, shown in Figure 4. The question is for which cases the dispersion relation (25) has purely imaginary roots. For simplicity we consider disturbance with ${\mathit{k}}_0=k_0\widehat{z}$ parallel to the velocities, the results can be easily generalized.

One could naively think that the extreme values of these instability regions can be found by putting $\omega =0$ in the dispersion relation (25) and require both numerators to vanish (since the denominators are opposite real numbers). However this means $\tilde{\kappa }=0$ and there is no x-dependence of the disturbance. That case corresponds to magnetosonic waves in the frame of each fluid with total wavevectors ${\mathit{k}}_t=k_0\widehat{z}$ and ${\omega }_0=\mp k_0{V}_0$ with $V_0^4-V_0^2(v_{\mathrm{A}}^2+c_s^2)+c_s^2{v}_{\mathrm{A}\Vert }^2=0$ (see Section 3.1). The absence of x dependence makes these cases unrelated to the unstable modes.

However, there are two other possibilities. The vanishing of both denominators of the dispersion relation (25), for $\omega =0$, needs to be considered as a possible limiting case. It corresponds to Alfvén waves with ${\omega }_0=\mp k_0{V}_0$ with $V_0^2=v_{\mathrm{A}\Vert }^2$ (see Section 3.2). Note that $\tilde{\kappa }$ does not enter in the dispersion relation. Nevertheless its value is nonzero ${\tilde{\kappa }}^2=k_0^2\left(1-{\displaystyle \frac{v_{\mathrm{A}\Vert }^2}{v_{\mathrm{A}}^2}}\right)$. The x dependence of the disturbance allows for a possible connection with unstable modes.

A third possibility is to realize that the limiting values may correspond to bifurcations of the dispersion relation. It is actually evident from inspection of Figure 3 that the slope ${\displaystyle \frac{\partial \omega }{\partial M}}$ becomes infinity when $\omega =0$. A dispersion relation in general depends on various parameters, and the slope is meant as a derivative with respect to one of them keeping constant all the others. For a dispersion relation of the form $f(\omega ,p_1,p_2,\cdots )=$ constant its differential ${\displaystyle {\displaystyle \frac{\partial f}{\partial \omega }}d\omega +\sum _n{\displaystyle \frac{\partial f}{\partial p_n}}d{p}_n=0}$ shows that bifurcation corresponds to ${\displaystyle \frac{\partial f}{\partial \omega }}=0$. Thus, extreme values of the instability regions may be connected to the condition ${\left.{\displaystyle \frac{\partial f}{\partial \omega }}\right|}_{\omega =0}=0$.

In our case there are three parameters, and we can choose $p_1=M^2={\displaystyle \frac{V_0^2}{c_s^2+v_{\mathrm{A}}^2}}$, $p_2={\displaystyle \frac{c_s^2}{c_s^2+v_{\mathrm{A}}^2}}$, $p_3={\displaystyle \frac{v_{\mathrm{A}\Vert }^2}{V_0^2}}$. The algebra gives that ${\left.{\displaystyle \frac{\partial f}{\partial \omega }}\right|}_{\omega =0}=0$ leads to $p_1-2+p_1{p}_3-2p_1{p}_2{p}_3^2-2p_2^2{p}_3^2+4p_2{p}_3=0$, or the following cubic for $V_0^2$

$$\begin{array}{c}\hfill V_0^6-(2c_s^2+2v_{\mathrm{A}}^2-v_{\mathrm{A}\Vert }^2)V_0^4+{\displaystyle \frac{2c_s^2{v}_{\mathrm{A}\Vert }^2}{c_s^2+v_{\mathrm{A}}^2}}(2c_s^2+2v_{\mathrm{A}}^2-v_{\mathrm{A}\Vert }^2)V_0^2-{\displaystyle \frac{2c_s^4{v}_{\mathrm{A}\Vert }^4}{c_s^2+v_{\mathrm{A}}^2}}=0,\end{array}$$

(28)

or, $M^6+({cos}^2\mathsf{\Theta }-2)(M^4-2{cos}^2\eta {cos}^2\mathsf{\Theta }{M}^2)=2{cos}^4\eta {cos}^4\mathsf{\Theta }$.

For the left panel of Figure 4 we have $cos\mathsf{\Theta }=sin\eta$, ${\displaystyle \frac{v_{\mathrm{A}}}{c_s}}=tan\eta$, and the extreme values of M are given by $M^6-(2-{sin}^2\eta )M^4+2{cos}^2\eta {sin}^2\eta (2-{sin}^2\eta )M^2-2{cos}^4\eta {sin}^4\eta =0$. The roots are $M=sin\eta$, $M=cos\eta \sqrt{1-\sqrt{cos\left(2\eta \right)}}$, $M=cos\eta \sqrt{1+\sqrt{cos\left(2\eta \right)}}$, and the corresponding wavenumbers ${\tilde{\kappa }}^2=0$, $k_0^{2}cos\left(2\eta \right)$, $k_0^{2}cos\left(2\eta \right)$, respectively. Thus the value $M=sin\eta$ corresponds to a magnetosonic wave without x dependence (unrelated to instability), while the other two solutions $M=cos\eta \sqrt{1-\sqrt{cos\left(2\eta \right)}}$, $M=cos\eta \sqrt{1+\sqrt{cos\left(2\eta \right)}}$ are the extreme values of M related to the unstable modes.

The result is that for ${\mathit{B}}_0\Vert {\mathit{k}}_0$ the Kelvin-Helmholtz instability occurs only if $v_{\mathrm{A}}<c_s$ and for velocities in the interval $c_s\sqrt{1-\sqrt{{\displaystyle \frac{c_s^2-v_{\mathrm{A}}^2}{c_s^2+v_{\mathrm{A}}^2}}}}<V_0<c_s\sqrt{1+\sqrt{{\displaystyle \frac{c_s^2-v_{\mathrm{A}}^2}{c_s^2+v_{\mathrm{A}}^2}}}}$. As we approach the limits the wavenumbers $\tilde{\kappa }$ approach real values and the instability is transformed to magnetosonic waves with constant amplitudes and wavenumbers ${\mathit{k}}_t=k_0\widehat{z}\pm k_0\sqrt{{\displaystyle \frac{c_s^2-v_{\mathrm{A}}^2}{c_s^2+v_{\mathrm{A}}^2}}}\widehat{x}$.

For the right panel of Figure 4 we have $\eta =\pi /2$ and ${\displaystyle \frac{B_{0z}}{B_{0y}}}=cot\mathsf{\Theta }$ and the maximum value of M corresponds to bifurcation for which $M^6-(1+{sin}^2\mathsf{\Theta })M^4=0$. The nontrivial root is $M=\sqrt{2-{\displaystyle \frac{B_{0z}^2}{B_0^2}}}$, and the corresponding wavenumbers ${\tilde{\kappa }}^2=k_0^2\left(1-{\displaystyle \frac{B_{0z}^2}{B_0^2}}\right)$.

The lower limit of the instability region corresponds to Alfvén waves $V_0^2=v_{\mathrm{A}\Vert }^2$ so $M={\displaystyle \frac{|B_{0z}|}{B_0}}$.

The result is that for the cold case $c_s=0$ the Kelvin-Helmholtz instability occurs only if ${\displaystyle \frac{|B_{0z}|}{\sqrt{{\rho }_0}}}<V_0<\sqrt{{\displaystyle \frac{2B_0^2-B_{0z}^2}{{\rho }_0}}}$. As we approach the limits the wavenumbers $\tilde{\kappa }$ approach real values and the instability is transformed to waves (Alfvén wave in the lower limit and magnetosonic in the upper) with constant amplitudes and wavenumbers ${\mathit{k}}_t=k_0\widehat{z}\pm k_0\sqrt{1-{\displaystyle \frac{B_{0z}^2}{B_0^2}}}\widehat{x}$.

6.6. Bifurcations in the General Case

The dispersion relation in the general case is $f={\displaystyle \frac{A_2{\tilde{\kappa }}_1}{A_1{\tilde{\kappa }}_2}}=1$, and taking its logarithmic derivative we conclude that bifurcations occur when ${\displaystyle \frac{\partial {\tilde{\kappa }}^2}{{\tilde{\kappa }}^2\partial \omega }}-2{\displaystyle \frac{\partial A}{A\partial \omega }}$ is continuous at the interface. Equivalently, since ${\displaystyle \frac{\tilde{\kappa }}{A}}$ is also continuous, ${\displaystyle \frac{\partial {\tilde{\kappa }}^2}{A^2\partial \omega }}-2{\tilde{\kappa }}^2{\displaystyle \frac{\partial A}{A^3\partial \omega }}$ is continuous. Substituting S and $A={\rho }_0{\omega }_0^2-F^2$ this quantity equals

${\displaystyle \frac{4{\rho }_0{\omega }_0}{A^3}}{\mathit{k}}_0^2-{\displaystyle \frac{4{\rho }_0{\omega }_0^3{F}^2}{A^3(A{c}_s^2+{\omega }_0^2{B}_0^2)}}-{\displaystyle \frac{2{\rho }_0{\omega }_0^5({\rho }_0{c}_s^2+B_0^2)}{A^2{(A{c}_s^2+{\omega }_0^2{B}_0^2)}^2}}$.

In the hydrodynamic case the latter simplifies to $2{\displaystyle \frac{2{\mathit{k}}_0^2{c}_s^2-{\omega }_0^2}{{\rho }_0^2{c}_s^2{\omega }_0^5}}$, while in the cold case to $2{\omega }_0{\displaystyle \frac{2{\mathit{k}}_0^2{v}_{\mathrm{A}}^2-k_0^2{v}_{\mathrm{A}\Vert }^2-{\omega }_0^2}{{\rho }_0^2{v}_{\mathrm{A}}^2{({\omega }_0^2-k_0^2{v}_{\mathrm{A}\Vert }^2)}^3}}$.

7. Summary

The main goal of the paper is to present the minimalist approach in stability problems with planar geometry that can be described using cartesian coordinates. It shows the power of the approach in finding the dispersion relation by integrating a single first-order differential equation, the principal equation.

The mathematical formalism is quite cumbersome, but can be made easier by defining intermediate quantities with important physical meaning. For example the complex $\tilde{\kappa }$ has a direct connection to wave propagation and simultaneous amplitude variation. The function S is connected to the compressibility that depends on both, the thermal and magnetic pressure, and shows how the wave propagation and the growth rate of the instability depend on them. The function F is the signature of the stabilizing nature of the magnetic tension. All these are discussed and reviewed when applying the method to the classical instabilities.

A more extensive analysis was done on the Kelvin-Helmholtz instability, for which new results were also found. Namely analytical solutions of the dispersion relation in certain cases, as well as the study of bifurcations and their connection to the ranges of instability. These tools make much easier the parametric study, which is far from being considered complete even in classical instabilities.

Detailed studies of specific problems in hydrodynamics or magnetohydrodynamics hopefully will benefit from the examples and the formalism presented, but also from the ideas how to explore possible existence of analytical solutions and specify the range of instability. Applications in other geometries and more general theoretical frameworks are also possible and will be presented in other connections.