Fluid Instabilities

Gravitational Instabilities in a Static vs. Expanding Medium

Let’s recall two of the dispersion relations we have derived:

Free Wave Dispersion Relation: Sound Waves

ω=vsk

Surface Water Waves

Here, we had gravity (constant), which gives a dispersive relation:

ω2=gktanh(kH)

where H is the depth of the sea. Also note that here, we had Φ=g.

And also, recall our equations for linear perturbation theory:

ρ=ρ0+ρ1
P=P0+P1
v=v0+v1
Φ=Φ0+Φ1

Note that ρ1ρ0, and so forth. We also had things like the adiabatic sound speed:

v2s=(Pρ)const. S=P1ρ1

Static, Uniform Medium

For a static, uniform medium, we have ρ0=const., P0=constant, and v0=0. The potential is interesting, though, which we will return to. The fluid equations are:

Zeroth Order Continuity Equation

ρ0t+(ρ0v0)=0

This is trivially true.

Zeroth Order Euler Equation

The LHS is 0 since v0=0. The RHS of the Euler equation tells us that:

1ρ0P0Φ0=0Φ0=0

Zeroth Poisson Equation

We have:

2Φ0=4πGρ0

Since we are radially symmetric, we have:

Φ0=4π3Gρ0r

Hey! Wait a minute! The two zeroth order equations tell us that ρ0 has to be 0. This isconsistency is actually called the Jeans-Swindle. The reason we have an inconsistency is because we need our medium to be infinitely large. Thus, the Jean’s Swindle is a reflection of the unrealistic background we are considering.

We now go to the perturbed equations:

First Order Continuity

ρ1t+ρ0v=0

First Order Euler

v1t=v2sρ0P1Φ1

First Order Poisson

2Φ1=4πGρ1

We can combine the first two equations, taking the time derivative of the first. Doing so, we get:

2ρ1t2+ρ0(v2sρ02ρ12Φ1)=0

Now use the Poisson equation (equation 3), to give a wave’s equation with a driving term:

2ρ1t2v2s2ρ14πGρ0ρ1=0

We use Fourier analysis, turning a differential into algebraic equation, since the modes decouple. We do this by saying ρ(r,t)ei(krωt)c(k)d3k . However, if we had non-linear terms like ρ21, all Hell breaks loose in Fourier space. This is because we have convolutions when things are non-linear, and that is ugly.

Anyways, we get:

ω2+v2sk24πGρ0=0

Dispersion Relation for Jean’s Instability

ω2=v2sk24πGρ0

This is a dispersive wave (since we are non-linear in k). We can also immediately see that ω might be imaginary! If pressure wins, we have waves/oscillatory behavior! If gravity wins, we get exponential collapse.

We can write this in a bit different way:

ω2=v2s(k2k2J)

where:

kJ4πGρ0v2s

is the Jean’s wavenumber. We also often see:

λJ2πkJ

And Jean’s mass (the mass enclosed within λJ):

MJ=4π3λ3Jρ0

And we often say things like ``what is the Jean’s mass/Jean’s length at a time in cosmic history?’’ Here, we compare the density to the sound speed, more or less, through the definition of Jean’s wavenumber.

Note that pre-CMB, the Jean’s Length is huge. After free electrons combine with ionized hydrogen (forming neutral hydrogen), photons free stream. Thus, the Jean’s Mass drops tremendously. The Jean’s Mass post-recombination is of order 106 solar masses. This is when structure formation takes off!

Two Regimes

Regime A

We consider here k>kJ or λ<λJ. In this case, we have ω2>0, and thus ω is real. ρ1 oscillates like sound waves. This is a stable solution.

Regime B

Here ω2<0, and thus ω is imaginary. This implies:

ρ1e±|ω|t

The growing solution is exponential!

An Side

Note that the dispersion relation for Jean’s instability is very similar to the plasma frequency:

ω2=v2sk2+4πnee2me>0

Thus the plasma frequency is always oscillatory!

An Expanding Medium

In this case, we have:

ρ0=ρ0(t)

We can define a δ such that:

δρ1ρ0=ρρ0ρ0=(δρ)

The overdensity can be as large as can be, but the underdensity can only get as negative as 1. This negative overdensity regions of the Universe are called voids. Thus (1δ<).

Also note that our fluid velocity is no longer 0. Instead, we have v0=Hr (Hubble expansion). This is also v=˙aar. By convection, we choose a=1 today, and a<1 before today. Note that taking:

v0=3˙aa

And note that

v1= peculiar velocity

We now go to our fluid equations:

Zeroth Order Continuity

ρ0t+ρ0v0=ρ0t+3˙aaρ0=0ρ0(t)1a(t)31volume

This is nice and in accordance with our intuition.

First Order Continuity

ρ1t+ρ0v+(ρ1v)=0

We’ll now want to go to δ coordinates in our first order continuity equation:

ρ1=ρ0δρ0˙δ+δ˙ρ0+ρ0v1+ρ0v0δ+ρ0δv=0

Also note that ˙ρ0=3˙aa from the continuity equation, and we know that v0=3˙aa. This cancels nicely:

Dividing through by ρ0, we get the linearized Continuity Equation for an expanding Universe:

˙δ+v1+˙aarδ=0

where we have replace v0 with Hubble’s Law.

Zeroth Order Euler Equation

Again, v0 is not zero, and so we just get a constra:

\partitalv0t+(v))v0=Φ0

First Order Euler Equation

v1t+(v0)+(v1)v0=v2sρ0ρ1Φ1

Let’s simplify the ugly term:

(v1)v0=˙aa(v1)r

And note that:

v1=v1xx+v1yy+v1zz

And so acting on r, we get: (v1)r=˙aav1.

This makes our First Order Euler equation;

v1t+˙aa(r)v1+˙aav1=v2sδΦ1

And finally, our Poisson:

Zeroth and First Order Poission

We get:

2Φ0=4πGρ0
2Φ1=4πGρ0δ

We no longer have Jean’s Swingle, by the way!


Here, we continue from last week. Let’s recall a few things:

Static Case: Pressure vs. Gravity

In this case, we derived:

ω2=v2sk24πGρ0

or

ω2=v2s(k2k2J)

where kJ is the Jeans Wavenumber, defined as:

kJ4πGρ0v2s

We got this by perturbing our fluid equations in the limit where the sub-0 pieces do not depend on time, nor space. We then moved into the contracting or expanding medium.

Expanding Case

The linearized fluid equations were dervied last time. These were:

From the continuity equation: $˙δ+v1+˙aa(r)δ=0$

From the Euler equation:

v1t+˙aa(r)v1+˙aav1=v2sδΦ

Note that we might have had this expression wrong last time. We had an extra ρ0 somewhere.

And from the Poisson equation:

2Φ1=4πGρ0δ

To solve these equations, we go into Fourier space. Also, we will go into the co-moving coordinate system which will cancel out many of the terms. First, doing the Fourier Transform:

To go to co-moving coordinates:

r=a(t)x

where r is physical and x are the co-moving coordinates. Note that a(today)=1. Correpsondingly, our co-moving wave number k will be paried with x and not r. Go to Fourier Space:

δ(r,t)
v1(r,t)
Φ1(r,t)

These three quanities are:

=d3keikra(t)fk(t)

where fk(t) are the Fourier transformed functions. In one shot, we will do a lot of algebra here. We will calculate ˙δ in real space:

˙δ(r,t)=d3keikra(t)(˙δk˙aaikra(t)δk)

And

δ(r,t)=ika(t)δk

And

v1(r,t)=ikavk

Now, with these, our CONTINUITY EQUATION becomes:

˙δk˙aaikraδk+ikvka+˙aa(rika)δk=0

We get cancellations, and thus:

˙δk+ikvka=0

Our Euler equation is:

˙vk˙aaikravk+˙aa(rika)vk+˙aavk=v2sikaδkikaΦk

We get cancellations and combinations (with total derivatives), giving the linearized euler equation!

t(avk)+iv2skδk+ikΦk=0

And our Poisson equation:

k2a2Φk=4πGρ0δk

And we now combine our three boxed equations in the most efficient way possible for an equation for δ. How do we do that?

1. We begin by getting rid of v. We do this by multiplying our first equation by a2. Doing so, this becomes:

a2˙δk+ik(avk)=0

2. And we now take t, and use the second equation to replace the first term. This gives, using the second equation above:

a2¨δk+2a˙a˙δk+ik(iv2skδkikΦk)=0

3. Divide by a2:

¨δk+2˙aa˙δk+v2sk2a2δk4πGρ0δkusing Poisson from above=0

And now we have one equation!

¨δk+2˙aa˙δk+v2sa2(k2k2J)δk=0

where

k2J4πGρ0a2v2s

Note that this is the co-moving Jeans wavenumber, and since r=ax, the a cancels out in the physical Jeans wavenumber. The key difference in the above equation is the ˙a/a term, by the way, which makes sense for our intuition. Let’s look at some regimes.

Examine k < k_J regime:

In this case, we have:

¨δk+2˙aa˙δk4πGρ0δk=0

In order to go further, we need to specify our Universe. The simplest case is a flat Universe (Einstein-de Sitter) model. In this case, Ωm=1 (non-relativistic matter, no curvature, no Cosmological Constant, etc). Recall (or derive), and remember that at2/3 for EdS Universe:

4πGρ=32H2(t)=˙a2a2=23t2

And lastly, our equation for ¨δ is:

¨δk+43t˙δk23t2δk=0

Our Ansatz is: δktα. This makes our expression:

α(α1)+43α23=0α2+α3=23
(α23)(α+1)=0α=23,1

What did we find here? We found that δt2/3a or δt1. These represents the growing and decaying mode respectively. We can compare these with the static Universe, in which δ is exponential e±ωt, and not a power law. The thing that slows the exponential growth is the drag term due to the expansion of the Universe.

The expansion slows to a power-law form, from an exponential form in the static case. All these different modes, by the way, decouple and have the same behavior. Density perturbations grow like the scale factor in the flat Universe case.

Take a look into the dark matter power spectrum, by the way! We derived the time-growth of the power spectrum! The power spectrum of dark matter grows with a2.

Note that we only solved these expressions for the pertubration growth, but we can do the same with the v and Φ.

Gravitational Instability of a Rotating (Self-Gravitating) Disk

Consider a uniformly rotating sheet of zero thickness with constant angular velocity Ω=Ωˆz. We assume a constant surface density Σ0. We also assume the sheet is in the (x,y) plane at z=0. Usually, it’s easier to go into the moving frame. Thus, we go into the rotating frame.

In the frame that rotates with the unperturbed sheet, consider small perturbations in the (x,y) plane.

Σ(x,y,t)=Σ0+Σ1(x,y,t)

And note that our 3D density, knowing that our disk is razor thin, is:

ρ=Σδ(z)

We have a fluid moving in the disk with velocity:

v(x,y,t)=v0+v1(x,y,t)

In the rotating frame, v0=0, and thus:

v(x,y,t)=v1(x,y,t)

As a consequence of going into the rotating frame, we will have to deal with fictitious forces later.

Our pressure is:

P(x,y,t)=P0+P1(x,y,t)

Note that, because we are in 2D, our sound speed is defined a bit differently, such that:

v2s=PΣ=P1Σ1

Which is a statement saying that pressure acts only in the plane, and not upward or donward.

And potential:

Φ(x,y,t)=Φ0+Φ1(x,y,t)

The fluid equations take the following form.

Continuity in 2javascript

\frac{\partial \Sigma}{\partial t} + \vec \nabla (\Sigma \vec v)= 0 \rightarrow \text{ Perturbed } \rightarrow \boxed{\frac{\partial \Sigma_1}{\partial t} + \Sigma_0 \vec \nabla \cdot \vec v = 0 }

Euler Equation

vt+(v)v=1ΣPΦ2Ω×vCoriolis+Ω2(xˆx+yˆy)centrifugal force

Note that this is where we have to pay the price of going to the rotating frame. Those extra terms might look ugly, but they’re not all that bad.

We can write down our perturbed pieces:

v1t=v2sΣ0Σ1Φ12Ω×v

Poisson

2Φ=4πGΣδ(z)

Perturbing this, we find:

2Φ1=4πGΣ1δ(z)

Note that the boxed equations are the perturbed equations. This is pretty similar to what we had before, but we have a nasty coriolis term which prevents our usual tricks.


Let’s continue. Recall from last time that:

Σ1t+Σ0v1=0
v1t=v2sΣ0Σ0Φ12Ω×v1
2Φ1=4πGΣ1δ(z)

from the Continuity equation, Euler’s equation, and Poisson equation, respectively. The coriolis force restricts us from doing our normal trick. Still, though, we will do a Fourier Transform since we have linearized equations. In this case, we have (x,y)(kx,ky). But, we can choose our orientation such that k=kˆx, i.e., ky=kz=0. Additionally, we have: kr=kx.

Doing so:

Σ1(x,y,t)=d3kΣ1(k)ei(kxωt)

And the same form for velocities v1x and v1y:

v1x(x,y,t)=d3kv1x(k)ei(kxωt)
v1y(x,y,t)=d3kv1y(k)ei(kxωt)
Φ1(x,y,t)=d3kΦ1(k)ei(kxωt) at z=0
Φ1(x,y,z,t)=d3kΦ1(k)ei(krωt) at z not 0

We will start with the Poisson equation first, keeping in mind it depends on where in z we are:

At z=0, we have:

2Φ1=4πGΣ1

At z0:

2Φ1=0

This is the Laplace equation, which we have solved before when solving for water waves. For z0, that is – above or below the plane – we have:

Φ1(x,y,z,t)ei(kxωt)f(z)Φ1(k)

Recall that the vertical structure f(z), we had:

k2f+d2fdz2=0f(z)e±kz

To have f be 0 at infinity, we need:

fek|z|

For z0: we can integarte the Poisson equation from slightly below ϵ to slightly above +ϵ:

lim

And thus:

\lim_{\epsilon \rightarrow 0} \left(\frac{\partial \Phi_1}{\partial z} \Bigg\vert_{z=+\epsilon} - \frac{\partial \Phi_1}{\partial z} \Bigg\vert_{z=-\epsilon}\right) = 4\pi G\Sigma_1

We now go into k space:

\lim_{\epsilon \rightarrow 0} \left(-ke^{-k|z|} \Phi_1\Big\vert_{z=+\epsilon} - k e^{k|z|}\Phi_1 \Big\vert_{z=-\epsilon}\right) = 4\pi G \Sigma_1
-2k \Phi_1 = 4\pi G \Sigma_1

We thus just showed that:

\boxed{\Phi_1 = -\frac{2\pi G}{k} \Sigma_1} \text{ in } k \text{ space}

We have used equation 3 from above to get this result. We now return to equations 1 and 2 from the Continuity and Euler’s equation.

From the first equation (continuity), we have:

- i \omega \Sigma_1 + \Sigma_0 i k v_{1x} = 0 \rightarrow \boxed{v_{1x} = \frac{\omega}{k}\frac{\Sigma_1}{\Sigma_0}}

From the Euler’s equation, let’s look at \hat{y} first:

\hat{y}: -i\omega v_{1y} = -2(\vec \Omega \times \vec v_1)_{y} = -2\Omega v_{1x} \rightarrow v_{1y} = \frac{2\Omega}{i\omega} v_{1x}

Note that the \Sigma_1 and \Phi_1 terms vanished since we chose \vec k = k \hat{x}. Our Fourier Transform by choice has no y-dependence, so taking the gradient gives 0 in non-x directions.

And now the \hat{x} component:

\hat{x} : -i \omega v_{1x} = -\frac{v_s^2}{\Sigma_0} i k\Sigma_1 - ik\Phi_1 - \underbrace{2(\vec \Omega \times \vec v_1)_{x}}_{-\Omega v_{1y}}

Now, let’s combine our three linearized equations, getting rid of \Phi_1 and the v’s, writing in terms of \Sigma_1. We will start with the epxression for \hat{x} direction:

-i \frac{\omega^2}{k}\frac{\Sigma_1}{\Sigma_0} = -\frac{ik v_s^2 \Sigma_1}{\Sigma_0} + i2\pi G \Sigma_1 + 2 \Omega (- i\frac{2\Omega}{k}\frac{\Sigma_1}{\Sigma_0})

Simplifying things a bit, we get:

\boxed{\omega^2 = v_s^2 k^2 - 2\pi G k \Sigma_0 + 4\Omega^2}

Dispersion Relation for Self-Gravitating, Rotating Disk

\omega^2 = \underbrace{v_s^2 k^2}_\text{pressure (stable)} - \underbrace{2\pi G k \Sigma_0}_\text{gravity (unstable)} + \underbrace{4\Omega^2}_\text{rotation (stable)}

Note that statements of stability are based on the sign of the term.

Let’s compare the above result with the spherical Jean’s dispersion relation:

\omega^2 = v_s^2 k^2 - 4\pi G\rho_0

Some Comments

If gravity dominates:

\omega^2 \approx -2 \pi G k \Sigma_0 < 0 \rightarrow \text{ always unstable}

Gravity and Pressure dominate (\Omega \sim 0):

\omega^2 = v_s^2 k^2 - 2\pi G k \Sigma_0 \rightarrow \omega^2 < 0 \text{ iff } \boxed{k<\frac{2\pi G\Sigma_0}{v_s^2} \equiv k_J}

Long-\lambda perturbations are unstable.

Gravity and Rotation dominate (v_s^2 \sim 0):

\omega^2 < 0 \text{ iff } \boxed{k > \frac{2\Omega^2}{\pi G \Sigma_0}}

Hey! Look at that! The two regimes have different equality signs. In this case, short-\lambda perturbations are unstable. This kinda of makes sense – its the overall rotation that stabilizes the system.

In general

Let’s complete the square

\omega^2 = v_s^2 \left(k - \frac{\pi G \Sigma_0}{v_s^2}\right)^2 + 4\Omega^2 - \left(\frac{\pi G \Sigma_0}{v_s}\right)^{2}

When do we have stable conditions? This is when \omega > 0 if:

2\Omega > \frac{\pi G \Sigma_0}{v_s}

Or, we can combine things:

\boxed{\frac{v_s \Omega}{G \Sigma_0} > \frac{\pi}{2} \sim 1.571}

Stability Criterion for Self-Gravitating Disk

\frac{v_s \Omega}{G \Sigma_0} > \frac{\pi}{2}

For a thin disk of (collisionless) stars rather than (collisional) fluid…

In this case, the sound speed v_s is replaced by the characterstic speed of the stars (the velocity dispersion). In this case, we have a stable disk of stars if:

\boxed{\frac{\sigma \Omega}{G \Sigma_0} > 1.68}

for a Maxwellian velocity distribution of the stars.

This is the Toomre (1964) relation. Thus,

For a disk with a more realistic vertical structure (say, exponentially stratified instead of razor thin)

Then you can show that the gas disk is stable if:

\boxed{\frac{v_s \Omega}{G \Sigma_0} > 1.06}

This also assumes an isothermal equation of state.

For a differentially rotating disk of gas…

In this case, \Omega is not a constant but has a radial profile \Omega(R). Often, we replace \Omega(R) with \Omega with the epicycle frequency:

\kappa^2 = R \frac{\partial (\Omega)^2}{\partial R} + 4 \Omega^2

This makes the dispersion relation:

\omega^2 = k^2 v_s^2 - 2\pi G k \Sigma_0 + \kappa^2

And the disk is stable if:

Toomre Q Parameter

\boxed{Q \equiv \frac{v_s \kappa}{\pi G \Sigma_0} > 1 }

Note that we use \sigma_R if we are talking about stars (radial velocity dispersion). Also note that we can teaase out a timescale, which is bascially the dynamical time. Think in terms of free-fall timescale.

This is called the Toomre’s Q parameter:

If you’re curious, check out Binney & Tremaine, Section 3.2. Epicycling is actually important in axisymmetric systems. This is because slightly elliptical orbits can be approximated as epicycling (to first-order Taylor expansions).

Examples of kappa:
\Omega(R) = \frac{V_{circ}(R)}{R} \rightarrow V_{circ} \equiv \sqrt{\frac{GM(<R)}{R}}

1. Point Mass:

Here, we have Keplerian orbits, and so:

V_{circ} \propto \sqrt{1/R} \rightarrow \Omega \propto R^{-3/2} \rightarrow \kappa^2 = \Omega^2

Thus, in the Keplerian case, there is just one frequency – the orbital frequency!

2. Flat Rotation Curves

For velocity to be constant, we need M(<R)\propto R, and thus \rho \propto r^{-2}. Note that V_{c} = R \Omega \rightarrow \Omega \propto R^{-1}, and finally, \kappa = \sqrt{2}\Omega.

3. Solid Body Rotation

Here, we have a rigid sheet with \rho = \text{ constant}. In this case, M \propto R^{3} and V_{c} \propto R. This also makes \Omega = \text{ constant}, and we recover that \kappa = 2\Omega, which we started with.

Thus \kappa ranges from roughly \Omega to 2\Omega for reasonable galaxy shapes.

An Aside on the Jeans Instability

We are continuing from last time, and we are starting with a back of the envelope calculation for the Jean’s length. The key result is the dispersion relation:

\omega^2 = k^2 v_s^2 - 4\pi G \rho_0 = v_s^2 (k^2 - k_J^2)

Well, we know that:

k_J \equiv \sqrt{\frac{4\pi G\rho_0}{v_s^2}}

Dimensionally, we can do the same thing to tease out G\rho / v_s^2 dependence.

Back of the Envelope Analysis

The underyling physics is the balance of pressure and gravity. Let’s compare the forces, or more concretely, the timescale over which they act. Instability occurs if gravity wins over pressure.

1. Option 1: Comparing the Forces (per volume)

Consider a fluid of density \rho , total mass M, and radius R. The fluid elements in this fluid has volume V. Thus:

\text{Gravity} \sim \frac{F_g}{V} \sim \frac{GM \rho}{r^2} \sim G\rho^2 R

Outward pressure is:

\text{Pressure} \sim \frac{F_p}{V} \sim \nabla P \sim v_s^2 \nabla \rho \sim \frac{v_s^2 \rho}{R}

Gravity wins if:

G \rho^2 R > \frac{v_s^2 \rho}{R} \rightarrow \boxed{R > \sqrt{\frac{v_s^2}{G \rho}}} \sim \lambda_J

2.Option 2: Comparing the Timescales

We have the timescale of gravity, the free-fall timescale, which is the time it takes to fall R:

R \sim a t_g^2 \sim \frac{GM}{R^2} t_g^2 \sim G\rho R t_g^2 \rightarrow t_g \sim \sqrt{\frac{R}{a}} \sim \boxed{\frac{1}{\sqrt{G\rho}}}

And the timescale for pressure is the sound-crossing time:

t_p \sim \frac{R}{v_s}

Gravity wins when t_g < t_P:

\frac{1}{\sqrt{G\rho}} < \frac{R}{v_s^2} \rightarrow R > \sqrt{\frac{v_s^2}{G\rho}}

which is exactly the same as we found above!

An Aside: The Coriolis Force

A quick physical insight! Recall that the Coriolis force is:

-2 \vec \Omega \times \vec v

Consider walking outward radially from a spinning disk and the forces!

Some final thoughts

We continue here from last time. Recall that, for a self-gravitating disk, we have:

\omega^2 = k^2 v_s^2 - 2\pi G k \Sigma_0 + 4 \Omega^2

Note that rotation and pressure are stabilizing, whereas gravity causes instabilities.

Rayleigh-Taylor Instability

fig15

This instability occurs when we have two fluids in a consatnt gravitational field \vec{g}. We will call the top fluid \rho_+ and the bottom fluid as \rho_{-}. What happens to perturbations traveling along the boundary?

fig13

Recall too that we defined \xi_{z} which is a displacement of the perturbed fluid velocities.

In the top fluid, we have \rho_+, v_+, and \vec v_+. With subscript minus for the lower fluid. We begin with our zeroth order equations:

Zeroth Order Equations

In the single fluid case (surface water waves), we had:

\vec v_0 = 0
\frac{1}{\rho_0} \frac{\partial P_0}{\partial z} = - g

We found, for water waves, that:

P_0 = -\rho_0 g z + \mathcal{C}

We can easily generalize this to two fluids! In the top fluid, we will have:

P_{+0} = - \rho_{+0} g z + \mathcal{C}_{+}

In the bottom fluid, we have

P_{-0} = - \rho_{-0} g z + \mathcal{C}_{-}

Our first order equations are:

First Order Equations

For the single fluid case, recall that we had:

\vec \nabla \cdot \vec v_1 = 0 \text{ (incompressible)}

And the Euler’s equation gave us:

\frac{\partial \vec v_1}{\partial t} = -\frac{1}{\rho_0} \vec \nabla P_1

Which, together, gave: \nabla^2 P = 0.

The perturbed pressure became:

P_1(x,z,t) = Ae^{i(kx-\omega t)} \left(e^{kz} + B e^{-kz}\right)

And the perturbed velocity was:

v_{1z} = \frac{Ak}{\rho_0 i \omega} e^{i(kx-\omega t)} \left(e^{kz} - B e^{-kz}\right)

Again, let us generalize to the two fluid case. The top fluid has to be finite, so:

P_{+1} \propto e^{-kz} \text{ z}>0

The top fluid has:

P_{-1} \propto e^{+kz} \text{ z}<0

in order to be finite.

\frac{\partial v_{\pm 1 z}}{\partial t} = - i \omega v_{\pm 1z}

We now use the second equation in this note block. This gives, for the top fluid:

-i\omega v_{+1z} = k \frac{P_{+1}}{\rho_{+0}}

And the bottom fluid:

-i\omega v_{-1z} = -k \frac{P_{-1}}{\rho_{-0}}

We will now apply the boundary conditions at the surface of the two fluids.

Boundary Conditions

1. Fluid displacements are continuous at z=0. Remember that we defined the fluid displacement such that:

\frac{d\vec \xi}{dt} = \vec v_1 \rightarrow \frac{d\xi_z}{dt} = v_{1z} \propto e^{-i\omega t}

We integrate this, giving:

\xi_z = \frac{i}{\omega} v_{1z}

and this is true for both \pm fluids at z=0 since this is continuous.

And we bring all of our work from above together for boundary condition 2:

2. Force of Top Fluid on Bottom = Force of Bottom Fluid on Top

a. without surface tension:

P_{+1,L} = P_{-1,L} \leftarrow \text{ Lagrangian perturbation}

b. with surface tension:

P_{+1,L} \underbrace{- T \frac{d^2 \xi_z}{dx^2}}_\text{downward} = \underbrace{P_{-1,L}}_\text{upward} \text{ at z}=0

Combining things…

Recall we had:

P_{1L} = P_1 + \vec \xi \cdot \vec \nabla P_0

in the single fluid case. In the top fluid:

P_{+1L} = P_{+1} + \xi_z \frac{\partial P_{+0}}{\partial z} = -\frac{i\omega}{k} \rho_{+0} v_{+1z} + \frac{i}{\omega}v_{+1z} \cdot (-g\rho_{+0}) = -i v_{+1z} \rho_{+0}\frac{\omega}{k}\left(1+\frac{gk}{\omega^2}\right)

And the bottom fluid has:

P_{-1L} = P_{-1} + \xi_z \frac{\partial P_{-0}}{\partial z} = +\frac{i\omega}{k} \rho_{-0} v_{-1z} + \frac{i}{\omega}v_{1z} \cdot (-g\rho_{+0}) = i v_{1z} \rho_{-0}\frac{\omega}{k}\left(1-\frac{gk}{\omega^2}\right)

And lastly, we use the equation where we consider surface tension. From this, we get:

T \frac{\partial^2 \xi_z}{\partial z^2} = \frac{i}{\omega}T \frac{\partial^2 v_{1z}}{\partial x^2} = -\frac{ik^2}{\omega} T v_{1z}

Using BC#2, we get:

-\rho_{+0} \frac{\omega}{k}(1+\frac{gk}{\omega^2}) + \frac{k^2}{\omega} T = \rho_{-0}\frac{\omega}{k}(1-\frac{gk}{\omega^2})

Let us re-arrange, multiplying by \omega k:

\boxed{\omega^2 \left(\rho_{+0} + \rho_{-0}\right) = gk\left(\rho_{-0} - \rho_{+0}\right) + k^3 T}

Rayleigh-Taylor Dispersion Relation

The dispersion relation for two fluids in a gravitational field for waves traveling on the boundary is thus:

\boxed{\omega^2 \left(\rho_{+0} + \rho_{-0}\right) = gk\left(\rho_{-0} - \rho_{+0}\right) + k^3 T}

We will now look at some different cases.

For denser fluids on bottom…

Here, \rho_{-0} > \rho_{+0}, \omega^2>0 always, and thus the perturbations/waves oscillate and travel as e^{-i\omega t} waves.

For extremely dense fluids on top…

Here, \rho_{-0} \gg \rho_{+0} (like ocean waves!). In this limit:

\omega^2 = gk+ \frac{T}{\rho_{-0}}k^3 > 0

And thus we have stable waves! More importantly, let’s recall the water wave lecture result which has \omega^2 = (gk + \frac{Tk^3}{\rho_0})\tanh(kH). For kH \ll 1, \tanh \rightarrow 1 (deep-water limit). Thus, we recover the deep-water limit!

When the bottom fluid is less dense…

This is the quintessential Rayleigh-Taylor instability. We want to know when \omega^2<0. This implies:

\boxed{k^2 < \frac{g (\rho_{+0} - \rho_{-0})}{T} \equiv k_{crit}^2}

The perturbation is exponentially-unstable. Wavelengths larger than the critical wavelength are unstable. Longer wavelengths are more unstable in the Rayleigh-Taylor instability. Surface tension helps to stabilize the instability.

When surface tension is negligible…

In this case, T \approx 0, and:

\omega^2 = gk\frac{\rho_{-0} - \rho_{+0}}{\rho_{-0} + \rho_{+0}}

which has \omega < 0 is always unstable if we are top-heavy.

In astrophysics…

In astrophysics, a downward \vec g is the same as an upward \vec a. In the case of a supernova remnant, this might be an astrophysical wind moving into the ISM. When the blast wave happens, less dense material moves into a more dense material, and we get RT instabilities!


Kelvin-Helmholtz Instability

fig16

We have a very similar setup to last time. This is the picture to have in mind:

fig14

The key difference, however, is that the two fluids are moving at different speeds, v_+ and v_{-}. This relative velocity gives rise to a vertical sheer. Modifying the Rayleigh-Taylor disperion relation, we get:

(\omega - k v_{+0})^2 \rho_{+0} + (\omega - kv_{-0})^2\rho_{-0} = gk \left(\rho_{-0}-\rho_{+0}\right) + k^3 T

And we skip some math to show when \omega is imaginary. Assuming that T is negligible, \omega is imaginary (i.e., unstable waves) if:

k > \frac{(\rho_{-0}^2 - \rho_{+0}^2) g }{\rho_{-0}\rho_{+0}(v_{+0} - v_{-0})^2} \equiv k_{crit}

There are some important things to note:

Notes

1. Modes with short wavelengths are thus Kelvin-Helmholtz unstable. If surface tensions T is important, T helps to stabilize short wavelength modes.

2. We can also turn off gravity by setting g=0. In this case, k_{crit} is 0, and all modes are unstable. Gravity is helping us maintain stability!

3. If the densities are very disparate, i.e., \rho_{-0} \gg \rho_{+0}, then k_{crit} is very large. Thus, only very small wavelenghts are unstable.

4. If \rho_{-0} \gtrsim \rho_{+0}, k_{crit} \approx 0, and KH is very common in these systems. There are many applications of this case, but on important one is fast moving dry air over slower moving, water-laden air over the ocean. This gives rise to billow clouds.

fig17

In Astrophysics:

AGN jets, or stellar jets, moving at high speed relative to the surrounding medium often display KH instabilities.

Thermal Instabilities

First, heat conduction – the transportation of heat (i.e., thermal energy) from hot to cold regions. In the simplest case, there is no mass motion, unlike convection. An empirical relation has been found, but it’s hard to derive from first principles. This relation is called Fourier’s Law of Heat Conduction.

\boxed{\underbrace{\vec F}_\text{heat flux} = -\underbrace{k}_\text{thermal conductivity} \underbrace{\vec \nabla T}_\text{temperature gradient}}

Note that heat flux has units of \text{energy}/(\text{Area} \cdot \text{ time}), and thus thermal conductivity has units of \text{Watts}/(\text{m}\cdot \text{ Kelvin}). Here are some example values:

  • Air : k = 0.024 \text{ W } \text{ m}^{-1} \text{K}^{-1}

  • Aluminum has k=200

  • Copper has k=400

  • Diamond has k=1000 or so…

  • Nano tubes/graphene sheets are a bit more special.

For our fluid equations, we will:

1. Turn on Thermal Conduction & Cooling/Heating. 2. Ignore gravity.

Our fluid equations become:

\frac{\partial \rho}{\partial t} + \vec \nabla \cdot (\rho \vec v) = 0
\frac{\partial \vec v}{\partial t} + (\vec v \cdot \vec \nabla) \vec v = -\frac{1}{\rho} \vec \nabla P
\frac{\partial \varepsilon}{\partial t} + \vec \nabla \cdot \vec F = \dot{Q}

where \dot{Q} represents other heating rates (per volume) and \varepsilon is the thermal energy density. The third equation can be re-written:

\frac{\partial \varepsilon}{\partial t} = \vec \nabla (k \vec \nabla T) + \dot{Q}

We can perturb the three equations above, which we won’t do here (see Pringle & King Ch. 8, or Clarke & Carswell Ch. 9, or Field (1965, ApJ)).

Skipping the math, we can quote the result:

For thermal instability to occur, we need:

\left(\frac{\partial \dot{Q}}{\partial T}\right)\Bigg\vert_\text{const. P} > 0 \rightarrow \text{ The Field Criterion}

This makes sense since, as T increases slightly, net \dot{Q} increases, meaning that we heat the material even more! This is the instability.

An Example

Consider the interstellar medium (ISM). We will write heating/cooling rates as:

\dot{Q} = \underbrace{\dot{Q}_{+}}_\text{heating} - \underbrace{\dot{Q}_{-}}_\text{cooling}

The typical cooling curve looks something like this:

NEED IMAGE FROM NOTES HERE

We thus have three phases of the ISM. When the ISM is near T_{1} (cool phase) in some regions. T_{3} (hot phase) is stable in other regions. This gives rise to multiphase ISM!

Comments

1. There is an E&M analogy to Fourier’s Law – this is Ohm’s Law!:

\vec J = \sigma \vec E

since the electric field can be written as the gradient of the electric potential.

2. Thermal Energy Conservation:

\frac{\partial \varepsilon}{\partial t} + \vec \nabla \cdot \vec F = 0

where \varepsilon is the energy density. and \vec F is the heat flux we introduced above. Note that \varepsilon \propto T most of the time, and thus we have:

\frac{\partial T}{\partial t} \propto \nabla^2 T

which is a diffusion equation. We diffuse heat transfer due to e^{-}-nuclei collisions.

3. This is hard to derive from first principles since this is a non-equilibrium process. But, people work on that!