The Rayleigh–Taylor instability, or RT instability (after Lord Rayleigh and G. I. Taylor), is an instability of an interface between two fluids of different densities, which occurs when the lighter fluid is pushing the heavier fluid.
Drazin (2002) pp. 50–51. This is the case with an interstellar cloud and shock system. The equivalent situation occurs when gravity is acting on two fluids of different density — with the dense fluid above a fluid of lesser density — such as water balancing on light oil.
Consider two completely plane-parallel layers of immiscible fluid, the heavier on top of the light one and both subject to the Earth's gravity. The equilibrium here is unstable to certain perturbations or disturbances. An unstable disturbance will grow and lead to a release of potential energy, as the heavier material moves down under the (effective) gravitational field, and the lighter material is displaced upwards. This was the set-up as studied by Lord Rayleigh. The important insight by G. I. Taylor was, that he realised this situation is equivalent to the situation when the fluids are accelerated (without gravity), with the lighter fluid accelerating into the heavier fluid. This can be experienced, for example, by accelerating a glass of water downward faster than the Earth's gravitational acceleration.
As the instability develops, downward-moving irregularities ('dimples') are quickly magnified into sets of inter-penetrating Rayleigh–Taylor fingers. Therefore the Rayleigh–Taylor instability is sometimes qualified to be a fingering instability. The upward-moving, lighter material is shaped like mushroom caps.

This process is evident not only in many terrestrial examples, from salt domes to weather inversions, but also in astrophysics and electrohydrodynamics. RT fingers are especially obvious in the Crab Nebula, in which the expanding pulsar wind nebula powered by the Crab pulsar is sweeping up ejected material from the supernova explosion 1000 years ago.
Note that the RT instability is not to be confused with the Plateau-Rayleigh instability (also known as Rayleigh instability) of a liquid jet. This instability, sometimes called the hosepipe (or firehose) instability, occurs due to surface tension, which acts to break a cylindrical jet into a stream of droplets having the same volume but lower surface area.

Linear stability analysis

Base state of the Rayleigh–Taylor instability. Gravity points downwards.

The inviscid two-dimensional Rayleigh–Taylor (RT) instability provides an excellent springboard into the mathematical study of stability because of the exceptionally simple nature of the base state.Drazin (2002) pp. 48–52. This is the equilibrium state that exists before any perturbation is added to the system, and is described by the mean velocity field $U(x,z)=W(x,z)=0,\backslash ,$ where the gravitational field is $\backslash textbf\{g\}=-g\backslash hat\{\backslash textbf\{z\}\}.\backslash ,$ An interface at $z=0\backslash ,$ separates the fluids of densities $\backslash rho\_G\backslash ,$ in the upper region, and $\backslash rho\_L\backslash ,$ in the lower region. In this section it is shown that when the heavy fluid sits on top, the growth of a small perturbation at the interface is exponential, and takes place at the rate

$\backslash exp(\backslash gamma\backslash ,t)\backslash ;,\; \backslash qquad\backslash text\{with\}\backslash quad\; \backslash gamma=\{\backslash sqrt\{\backslash mathcal\{A\}g\backslash alpha\}\}\; \backslash quad\backslash text\{and\}\backslash quad\; \backslash mathcal\{A\}=\backslash frac\{\backslash rho\_\{\backslash text\{heavy\}\}-\backslash rho\_\{\backslash text\{light\}\}\}\{\backslash rho\_\{\backslash text\{heavy\}\}+\backslash rho\_\{\backslash text\{light\}\}\},\backslash ,$

where $\backslash gamma\backslash ,$ is the temporal growth rate, $\backslash alpha\backslash ,$ is the spatial wavenumber and $\backslash mathcal\{A\}\backslash ,$ is the Atwood number.

The perturbation introduced to the system is described by a velocity field of infinitesimally small amplitude, $(u\text{'}(x,z,t),w\text{'}(x,z,t)).\backslash ,$ Because the fluid is assumed incompressible, this velocity field has the streamfunction representation

$\backslash textbf\{u\}\text{'}=(u\text{'}(x,z,t),w\text{'}(x,z,t))=(\backslash psi\_z,-\backslash psi\_x),\backslash ,$

where the subscripts indicate partial derivatives. Moreover, in an initially stationary incompressible fluid, there is no vorticity, and the fluid stays irrotational, hence $\backslash nabla\backslash times\backslash textbf\{u\}\text{'}=0\backslash ,$. In the streamfunction representation, $\backslash nabla^2\backslash psi=0.\backslash ,$ Next, because of the translational invariance of the system in the x-direction, it is possible to make the ansatz

$\backslash psi\backslash left(x,z,t\backslash right)=e^\{i\backslash alpha\backslash left(x-ct\backslash right)\}\backslash Psi\backslash left(z\backslash right),\backslash ,$

where $\backslash alpha\backslash ,$ is a spatial wavenumber. Thus, the problem reduces to solving the equation

$\backslash left(D^2-\backslash alpha^2\backslash right)\backslash Psi\_j=0,\backslash ,\backslash ,\backslash ,\backslash \; D=\backslash frac\{d\}\{dz\},\backslash ,\backslash ,\backslash ,\backslash \; j=L,G.\backslash ,$

The domain of the problem is the following: the fluid with label 'L' lives in the region c, which in turn determiens the stability properties of the system.

The first of these conditions is provided by details at the boundary. The perturbation velocities $w\text{'}\_i\backslash ,$ should satisfy a no-flux condition, so that fluid does not leak out at the boundaries $z=\backslash pm\backslash infty.\backslash ,$ Thus, $w\_L\text{'}=0\backslash ,$ on $z=-\backslash infty\backslash ,$, and $w\_G\text{'}=0\backslash ,$ on $z=\backslash infty\backslash ,$. In terms of the streamfunction, this is

$\backslash Psi\_L\backslash left(-\backslash infty\backslash right)=0,\backslash qquad\; \backslash Psi\_G\backslash left(\backslash infty\backslash right)=0.\backslash ,$

The other three conditions are provided by details at the interface $z=\backslash eta\backslash left(x,t\backslash right)\backslash ,$.
Continuity of vertical velocity: At $z=\backslash eta$, the vertical velocities match, $w\text{'}\_L=w\text{'}\_G\backslash ,$. Using the streamfunction representation, this gives

$\backslash Psi\_L\backslash left(\backslash eta\backslash right)=\backslash Psi\_G\backslash left(\backslash eta\backslash right).\backslash ,$

Expanding about $z=0\backslash ,$ gives

$\backslash Psi\_L\backslash left(0\backslash right)=\backslash Psi\_G\backslash left(0\backslash right)+\backslash text\{H.O.T.\},\backslash ,$

where H.O.T. means 'higher-order terms'. This equation is the required interfacial condition.
The free-surface condition: At the free surface $z=\backslash eta\backslash left(x,t\backslash right)\backslash ,$, the kinematic condition holds:

$\backslash frac\{\backslash partial\backslash eta\}\{\backslash partial\; t\}+u\text{'}\backslash frac\{\backslash partial\backslash eta\}\{\backslash partial\; x\}=w\text{'}\backslash left(\backslash eta\backslash right).\backslash ,$

Linearizing, this is simply

$\backslash frac\{\backslash partial\backslash eta\}\{\backslash partial\; t\}=w\text{'}\backslash left(0\backslash right),\backslash ,$

where the velocity $w\text{'}\backslash left(\backslash eta\backslash right)\backslash ,$ is linearized on to the surface $z=0\backslash ,$. Using the normal-mode and streamfunction representations, this condition is $c\; \backslash eta=\backslash Psi\backslash ,$, the second interfacial condition.
Pressure relation across the interface: For the case with surface tension, the pressure difference over the interface at $z=\backslash eta$ is given by the Young–Laplace equation:

$p\_G\backslash left(z=\backslash eta\backslash right)-p\_L\backslash left(z=\backslash eta\backslash right)=\backslash sigma\backslash kappa,\backslash ,$

where σ is the surface tension and κ is the curvature of the interface, which in a linear approximation is

$\backslash kappa=\backslash nabla^2\backslash eta=\backslash eta\_\{xx\}.\backslash ,$

Thus,

$p\_G\backslash left(z=\backslash eta\backslash right)-p\_L\backslash left(z=\backslash eta\backslash right)=\backslash sigma\backslash eta\_\{xx\}.\backslash ,$

However, this condition refers to the total pressure (base+perturbed), thus

$\backslash left[P\_G\backslash left(\backslash eta\backslash right)+p\text{'}\_G\backslash left(0\backslash right)\backslash right]-\backslash left[P\_L\backslash left(\backslash eta\backslash right)+p\text{'}\_L\backslash left(0\backslash right)\backslash right]=\backslash sigma\backslash eta\_\{xx\}.\backslash ,$

(As usual, The perturbed quantities can be linearized onto the surface z=0.) Using hydrostatic balance, in the form

$P\_L=-\backslash rho\_L\; g\; z+p\_0,\backslash qquad\; P\_G=-\backslash rho\_G\; gz\; +p\_0,\backslash ,$

this becomes

$p\text{'}\_G-p\text{'}\_L=g\backslash eta\backslash left(\backslash rho\_G-\backslash rho\_L\backslash right)+\backslash sigma\backslash eta\_\{xx\},\backslash qquad\backslash text\{on\; \}z=0.\backslash ,$

The perturbed pressures are evaluated in terms of streamfunctions, using the horizontal momentum equation of the linearised Euler equations for the perturbations,

$\backslash frac\{\backslash partial\; u\_i\text{'}\}\{\backslash partial\; t\}\; =\; -\; \backslash frac\{1\}\{\backslash rho\_i\}\backslash frac\{\backslash partial\; p\_i\text{'}\}\{\backslash partial\; x\}\backslash ,$ with $i=L,G,\backslash ,$

to yield

$p\_i\text{'}=\backslash rho\_i\; c\; D\backslash Psi\_i,\backslash qquad\; i=L,G.\backslash ,$

Putting this last equation and the jump condition on $p\text{'}\_G-p\text{'}\_L$ together,

$c\backslash left(\backslash rho\_G\; D\backslash Psi\_G-\backslash rho\_L\; D\backslash Psi\_L\backslash right)=g\backslash eta\backslash left(\backslash rho\_G-\backslash rho\_L\backslash right)+\backslash sigma\backslash eta\_\{xx\}.\backslash ,$

Substituting the second interfacial condition $c\backslash eta=\backslash Psi\backslash ,$ and using the normal-mode representation, this relation becomes

$c^2\backslash left(\backslash rho\_G\; D\backslash Psi\_G-\backslash rho\_L\; D\backslash Psi\_L\backslash right)=g\backslash Psi\backslash left(\backslash rho\_G-\backslash rho\_L\backslash right)-\backslash sigma\backslash alpha^2\backslash Psi,\backslash ,$

where there is no need to label $\backslash Psi\backslash ,$ (only its derivatives) because $\backslash Psi\_L=\backslash Psi\_G\backslash ,$
at $z=0.\backslash ,$
Solution

Now that the model of stratified flow has been set up, the solution is at hand. The streamfunction equation $\backslash left(D^2-\backslash alpha^2\backslash right)\backslash Psi\_i=0,\backslash ,$ with the boundary conditions $\backslash Psi\backslash left(\backslash pm\backslash infty\backslash right)\backslash ,$ has the solution

$\backslash Psi\_L=A\_L\; e^\{\backslash alpha\; z\},\backslash qquad\; \backslash Psi\_G\; =\; A\_G\; e^\{-\backslash alpha\; z\}.\backslash ,$

The first interfacial condition states that $\backslash Psi\_L=\backslash Psi\_G\backslash ,$ at $z=0\backslash ,$, which forces $A\_L=A\_G=A.\backslash ,$ The third interfacial condition states that

$c^2\backslash left(\backslash rho\_G\; D\backslash Psi\_G-\backslash rho\_L\; D\backslash Psi\_L\backslash right)=g\backslash Psi\backslash left(\backslash rho\_G-\backslash rho\_L\backslash right)-\backslash sigma\backslash alpha^2\backslash Psi.\backslash ,$

Plugging the solution into this equation gives the relation

$Ac^2\backslash alpha\backslash left(-\backslash rho\_G-\backslash rho\_L\backslash right)=Ag\backslash left(\backslash rho\_G-\backslash rho\_L\backslash right)-\backslash sigma\backslash alpha^2A.\backslash ,$

The A cancels from both sides and we are left with

$c^2=\backslash frac\{g\}\{\backslash alpha\}\backslash frac\{\backslash rho\_L-\backslash rho\_G\}\{\backslash rho\_L+\backslash rho\_G\}+\backslash frac\{\backslash sigma\backslash alpha\}\{\backslash rho\_L+\backslash rho\_G\}.\backslash ,$

To understand the implications of this result in full, it is helpful to consider the case of zero surface tension. Then,

$c^2=\backslash frac\{g\}\{\backslash alpha\}\backslash frac\{\backslash rho\_L-\backslash rho\_G\}\{\backslash rho\_L+\backslash rho\_G\},\backslash qquad\; \backslash sigma=0,\backslash ,$

and clearly

• If , $c^2>0\backslash ,$ and c is real. This happens when the

lighter fluid sits on top;

• If $\backslash rho\_G>\backslash rho\_L\backslash ,$, and c is purely imaginary. This happens

when the heavier fluid sits on top.

Now, when the heavier fluid sits on top, , and

$c=\backslash pm\; i\; \backslash sqrt\{\backslash frac\{g\backslash mathcal\{A\}\}\{\backslash alpha\}\},\backslash qquad\; \backslash mathcal\{A\}=\backslash frac\{\backslash rho\_G-\backslash rho\_L\}\{\backslash rho\_G+\backslash rho\_L\},\backslash ,$

where $\backslash mathcal\{A\}\backslash ,$ is the Atwood number. By taking the positive solution, we see that the solution has the form

$\backslash Psi\backslash left(x,z,t\backslash right)=Ae^\{-\backslash alpha|z|\}\backslash exp\backslash left[i\backslash alpha\backslash left(x-ct\backslash right)\backslash right]=A\backslash exp\backslash left(\backslash alpha\backslash sqrt\{\backslash frac\{g\backslash tilde\{\backslash mathcal\{A\}\}\}\{\backslash alpha\}\}t\backslash right)\backslash exp\backslash left(i\backslash alpha$

x-\alpha|z|\right)\,
and this is associated to the interface position η by: $c\backslash eta=\backslash Psi.\backslash ,$ Now define $B=A/c.\backslash ,$

Hydrodynamics simulation of a single "finger" of the Rayleigh–Taylor instability Note the formation of Kelvin–Helmholtz instabilities, in the second and later snapshots shown (starting initially around the level $y=0$). As well as the formation of a "mushroom cap" at a later stage in the third and fourth frame in the sequence.

The time evolution of the free interface elevation $z\; =\; \backslash eta(x,t),\backslash ,$ initially at $\backslash eta(x,0)=\backslash Re\backslash left\backslash \{B\backslash ,\backslash exp\backslash left(i\backslash alpha\; x\backslash right)\backslash right\backslash \},\backslash ,$ is given by:

$\backslash eta=\backslash Re\backslash left\backslash \{B\backslash ,\backslash exp\backslash left(\backslash sqrt\{\backslash mathcal\{A\}g\backslash alpha\}\backslash ,t\backslash right)\backslash exp\backslash left(i\backslash alpha\; x\backslash right)\backslash right\backslash \}\backslash ,$

which grows exponentially in time. Here B is the amplitude of the initial perturbation, and $\backslash Re\backslash left\backslash \{\backslash cdot\backslash right\backslash \}\backslash ,$ denotes the real part of the complex valued expression between brackets.

In general, the condition for linear instability is that the imaginary part of the "wave speed" c be positive. Finally, restoring the surface tension makes c² less negative and is therefore stabilizing. Indeed, there is a range of short waves for which the surface tension stabilizes the system and prevents the instability forming.

Late-time behaviour

The analysis of the previous section breaks down when the amplitude of the perturbation is large. The growth then becomes non-linear as the spikes and bubbles of the instability tangle and roll up into vortices. Then, as in the figure, numerical simulation of the full problem is required to describe the system.

See also

• Richtmyer-Meshkov instability

• Kelvin–Helmholtz instability

• Mushroom cloud

• Plateau-Rayleigh instability

• Salt fingering

• Hydrodynamic stability

• Kármán vortex street