Waves and instabilities
Starting from the simplest
Jiˇr´ı Krtiˇcka
Masaryk University
Sound waves
Sound waves
Let us assume that the hydrodynamical equations
∂ρ
∂t
+ ∇ · (ρv) = 0,
ρ
∂v
∂t
+ ρv · ∇v = −∇p + ρg,
have static solution ρ0 = const. with g = 0.
1
Sound waves
Let us search for a small perturbation δρ ≪ ρ0 of the static solution in a
form of ρ = ρ0 + δρ and v = δv, which fullfills the hydrodynamical
equations:
∂(ρ0 + δρ)
∂t
+ ∇ · ((ρ0 + δρ)δv) = 0,
(ρ0 + δρ)
∂δv
∂t
+ (ρ0 + δρ)δv · ∇δv = −∇(p0 + δp).
Neglecting second-order terms we derive
∂δρ
∂t
+ ρ0∇ · δv = 0,
ρ0
∂δv
∂t
+ ∇δp = 0.
Derivating the first equation with respect to t, inserting from the second
one and rewritting δp = dp
dρ δρ ≡ a2
δρ we arrive at the wave equation
∂2
δρ
∂t2
− a2
∇2
δρ = 0.
2
The sound speed
The constant in the wave equation is the sound speed:
a =
dp
dρ
.
For isothermal perturbations we derive from the perfect gas equation of
state
a =
dp
dρ
=
dp
dρ T
=
kT
µmH
,
where µ is the mean molecular weight and mH is the mass of hydrogen
atom. For fully ionized hydrogen µ = 1
2 and a = 2kT/(mH).
For adiabatic perturbations we have
a =
dp
dρ
=
dp
dρ S
=
κkT
µmH
,
where κ is the specific heat ratio. For fully ionized hydrogen
a = 10kT/(5mH).
3
Sound waves revisited: acoustic cutoff
Let us again study the sound waves in an atmosphere, but this time the
atmosphere is statified in a homogeneous gravitational field directioned
along the z axis with g as gravity acceleration. Corresponding
hydrodynamical equations are
∂ρ
∂t
+
∂
∂z
(ρv) = 0,
∂v
∂t
+ v
∂v
∂z
= −
1
ρ
∂p
∂z
+ g,
In an isothermal atmosphere p = a2
ρ. The static density distribution
follows the barometric law ρ0 ∼ e−z/H
with the density scale-height is
H = a2
/g. Again, we shall consider small perturbations of the stationary
solution in the form of ρ = ρ0 + δρ and v = δv. However, this time we
will assume that the wavelength of perturbations can be comparable with
H. Therefore, we shall also take care of vertical derivative of the density.
4
Sound waves revisited: the dispersion relation
Assuming for simplicity just the vertical variations gives
δ ˙ρ + ρ0δv′
−
ρ0δv
H
= 0,
ρ0δ ˙v = −a2
δρ′
− gδρ.
Derivating the Euler equation with respect to time and inserting from the
continuity equation (twice) gives
δ¨v = a2
δv′′
+ gδv′
.
For g = 0 we recover the ordinary sound waves. Inserting the harmonic
waves δv ∼ ei(ωt−kz)
gives the dispersion relation
ω2
− a2
k2
+ igk = 0.
5
Acoustic cutoff frequency
The dispersion relation shall be solved assuming ω ∈ R and k ∈ C.
Setting the imaginary part of the dispersion relation to zero gives
Im(k) =
g
2a2
=
1
2H
.
The real part of the dispersion relation now requires
ω2
− a2
Re2
(k) − ω2
a = 0,
where ωa = a/(2H) is the accoustic cutoff frequency. For Earth’s
atmosphere 2π/ωa ≈ 7 min. Inserting this in the expression for δv gives
δv ∼ e
z
2H .
Therefore, the velocity perturbation grows with height eventually
steepening into the shock. This contributes to the heating of the upper
atmosphere. Moreover, from the real part of the dispersion relation
follows that waves can propagate vertically only for
ω > ωa.
6
Characteristics of differential equations
Characteristic direction
Let us assume that f = f (x, y). Then a linear combination afx + bfy
(where fx = ∂f /∂x) is a directional derivative of f along the direction
dx : dy = a : b. If (x(σ), y(σ)) is a curve parameterized by σ,
xσ : yσ = a : b, then afx + bfy is a directional derivative along the curve.
Let us consider system of 2 equations for two functions u(x, y), v(x, y):
L1 ≡ A11ux + B11uy + A12vx + B12vy + C1 = 0,
L2 ≡ A21ux + B21uy + A22vx + B22vy + C2 = 0.
We ask for a linear combination
L = λ1L1 + λ2L2
so that in the differential expression L the derivatives of u and v combine
to derivatives in the same direction. Such direction is characteristic.
7
Characteristic relations
Suppose that the characteristic direction is given by the above ratio
xσ : yσ. Then the condition that u and v are differentiated in L in the
same direction is
λ1A11 +λ2A21 : λ1B11 +λ2B21 = λ1A12 +λ2A22 : λ1B12 +λ2B22 = xσ : yσ.
This gives the system of equations for λ1 and λ2
ˆM
λ1
λ2
= 0, where ˆM =
A11yσ − B11xσ A21yσ − B21xσ
A12yσ − B12xσ A22yσ − B22xσ
leading to characteristic relations. The system has a non-trivial solution if
det ˆM = 0. This gives equation in a form of
ay2
σ − 2bxσyσ + cx2
σ = 0.
For ac − b2
> 0, this cannot be satisfied by any direction. Such equations
are called elliptic. For ac − b2
< 0 we have two characteristic directions.
Such systems are called hyperbolic. There are two sets of equations
dy
dx = ξ+ and dy
dx = ξ− defining two sets of characteristics C+ and C−.
8
Characteristic relations for 1D flow
For 1D flow ρ = ρ(x, t) ≡ u and v = v(x, t) the corresponding system is
L1 ≡ ρt + vρx + ρvx = 0,
L2 ≡ vt + v vx +
a2
ρ
ρx = 0.
This gives (t ≡ y)
ˆM =
vtσ − xσ
a2
ρ tσ
ρtσ vtσ − xσ
.
From det ˆM = 0 the characteristic relation is (vtσ − xσ)2
− a2
t2
σ = 0, or
(v ± a)tσ = xσ.
The characteristics correspond to sound waves.
9
Characteristic relations for 1D flow
The relation between λ1 and λ2 can be derived, e.g., from the first
equation of the system ˆMλ = 0
(vtσ − xσ)λ1 +
a2
ρ
tσλ2 = 0,
which, after insterting the characteristic relation, simplifies to
λ1 = ±
a
ρ
λ2.
Therefore, the linear combination of hydrodynamical equations is
L = λ1L1 + λ2L2 =
a
ρ
[±ρt + (a ± v)ρx ] + vt + (v ± a)vx = 0,
where we further selected λ2 = 1. As we can see, ρ and v are
differentiated in the same direction.
10
Transformation to ordinary differential equation
Because ρσ = ρttσ + ρxxσ = [ρt + (v ± a)ρx ] tσ from the characteristic
relation, by doing the linear combination we have transformed the original
system of partial differential equations to a system of ordinary differential
equation with new variables α ≡ σ (for + root) and β ≡ σ (for − root)
a
ρ
ρα + vα = 0,
−
a
ρ
ρβ + vβ = 0,
11
Domain of dependence
The solution of hydrodynamical equations define two sets of
characteristics
(v ± a)tσ = xσ.
Let us assume that the initial
conditions are given on curve J .
There are two characteristics that go
throught a selected point P. The
line AB intercepted by the two
characteristics is called domain of
dependence of P. This can be
utilized for a numerical integration.
12
Range of influence
The solution of hydrodynamical equations define two sets of
characteristics
(v ± a)tσ = xσ.
Let us assume that the initial
conditions are given on curve J .
The range of influence of a point Q
is the totality of points which are
influenced by the initial data at the
point Q. The range of influence of
the point Q consists of all points P
whose domain of dependence
contains Q. The range of influence
of influence is defined by two
characteristic drawn through Q.
13
Expansion of a gas: withdrawing piston
Let us study a tube filled with a gas bounded by a piston withdrawing
subsonically.
The piston starts at O and recedes
towards left causing an expansion of
the gas. The gas adjanced to a
piston moves with the same velocity
as the piston. Only one set of
characteristics drawn from piston
propagates into the flow. The flow
in a zone I is not influenced by a
moving piston. The flow in a zone II
is a rarefaction wave.
14
Compression of a gas: advancing piston
Let us study a tube filled with a gas bounded by a piston advancing
subsonically.
The piston starts at O and moves
towards righ causing an compression
of the gas. The gas adjanced to a
piston moves with the same velocity
as the piston. Only one set of
characteristics drawn from piston
propagates into the flow. The flow
in a zone I is not influenced by a
moving piston. Intersecting
characteristic form an envelope. The
solution is not unique at the
intersection. This leads to a
formation of a shock wave.
15
Characteristics of more than two equations
We consider n differential equations
Li ≡ Aij
∂uj
∂x
+ Bij
∂uj
∂y
+ Cj , i = 1, · · · , n.
We ask for a linear combination
L = λi Li
so that in the differential expression L the derivatives of uj
combine to
derivatives in the same direction. This gives the conditions
λi Aij : λi Bij = xσ : yσ, j = 1, · · · , n,
or
λi (Aij yσ − Bij xσ) = 0, j = 1, · · · , n.
This system has a non-trivial solution if
det |Aij yσ − Bij xσ| = 0.
16
Application: Charactersistics including the energy equation
For isentropic 1D flow ρ = ρ(x, t) ≡ u1
, v = v(x, t) ≡ u2
, and
s = s(x, t) ≡ u3
the corresponding system of equations is
L1 ≡ ρt + vρx + ρvx = 0,
L2 ≡ vt + v vx +
a2
ρ
ρx = 0,
L3 ≡ st + vsx = 0.
ˆM =



vtσ − xσ
a2
ρ tσ 0
ρtσ vtσ − xσ 0
0 0 vtσ − xσ


 .
From det ˆM = 0 the first two characteristic relations are the same as
without energy equation
(v ± a)tσ = xσ,
vtσ = xσ.
This corresponds to sound waves propagating at the sound speed and
entropy wave propagating at zero speed (with respect to the flow). 17
Gravity waves
Gravity waves
We will study the propagation of waves in an atmosphere, which is in
hydrostatic equilibrium given by external gravitational field.
We will study the movement of a
blob with density ρ in hydrostatic
equilibrium with outside medium
with density ρ0(ξ). We will assume
the density gradient in the
atmosphere dρ
dz
at
and neglect the
heat exchange between the blob and
the atmosphere: the processes are
adiabatic. The equation of motion
including the buoyancy is:
ρ
d2
ξ
dt2
= −g(ρ − ρ0).
0000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000
00000000000000000000000000
000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
111111111111111111111111111111111111111
11111111111111111111111111
111111111111111111111111111111111111111
111111111111111111111111111111111111111
11111111111111111111111111
11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
0 ρ
ξ
0
g
000000000000000
00000000000000000000
111111111111111
11111111111111111111
ρ
18
Gravity waves
In the equation of motion,
ρ
d2
ξ
dt2
= −g(ρ − ρ0),
we shall use the Taylor expansion to derive the buoyancy term,
ρ0(ξ) = ρ0(0) +
dρ
dz at
ξ,
ρ(ξ) = ρ(0) +
dρ
dz ad
ξ.
Because the blob is initially in equilibrium, ρ0(0) = ρ(0), the equation of
motion
d2
ξ
dt2
= −ω2
BVξ
describes an oscillatory motion, so-called gravity waves. The frequency of
oscillations,
ω2
BV =
g
ρ
dρ
dz ad
−
dρ
dz at
is the Brunt-V¨ais¨al¨a frequency. 19
Gravity waves: the case of ω2
BV > 0
For dρ
dz
ad
> dρ
dz
at
, i.e., for dρ
dz
ad
< dρ
dz
at
we have ω2
BV > 0.
The initial perturbation results in oscillations ξ(t) = ξ0e±i|ωBV|t
.
Gravity waves in Earth’s atmosphere:
20
Gravity waves: the case of ω2
BV > 0: Earth’s atmosphere
21
Gravity waves: the case of ω2
BV > 0: Earth’s atmosphere
22
Gravity waves: solar differential rotation
Angular velocity as a function of radius in the Sun from accoustic modes
in heliseismic observations (Turck-Chi´eze). The lack of differential
rotation in the radiative zone is due to angular momentum transport by
gravity waves (Charbonnel & Talon 2005).
23
Gravity waves: the case of ω2
BV < 0
For dρ
dz ad
< dρ
dz at
, i.e., for dρ
dz ad
> dρ
dz at
we have ω2
BV < 0.
The initial perturbation results in instability ξ(t) = ξ0e±|ωBV|t
.
The instability leads to convection.
24
Schwarzschild stability criterion
The stability criterion can be recast in another intuitive form. From the
ideal gas equation of state,
dρ
dz at
=
ρ
p
dp
dz at
−
ρ
T
dT
dz at
.
The convective plumes are in hydrostatic equilibrium with the
sorrounding environment meaning that
dp
dz at
≡
dp
dz ad
= γ
p
ρ
dρ
dz ad
.
Therefore, the stability criterion is
(1 − γ)
dρ
dz ad
> −
ρ
T
dT
dz at
.
From the adiabatic equation follows that (dρ/dz)ad =
= 1/ (γ − 1) ρ/T (dT/dz)ad, which yields
dT
dz ad
<
dT
dz at
,
which is Schwarzschild stability criterion. 25
Temperature distribution in a convective atmosphere
The convective motions are typically slower than the sound waves
maintaining the hydrostatic equilibrium, therefore one can use
dp
dz
= −ρg
to determine the temperature gradient. The pressure is
p = a2
ρ = kTρ/µ, which gives
dT
dr
+
T
ρ
dρ
dr
= −
µg
k
.
The adiabatic equation Tρ1−γ
= const. gives
T
ρ
dρ
dr
=
1
γ − 1
dT
dr
,
which yields for the temperature gradient
dT
dr
= −
γ − 1
γ
µg
k
.
This predicts the temperature gradient of −10 Kkm−1
for the atmosphere
of our Earth and about 5 × 106
K R−1
⊙ for the envelope of Sun. 26
Simulation of convection
27
Solar granulation
28
Kelvin-Helmholtz & Rayleigh-Taylor
K-H & R-T instabilities: the initial setup
Consider a shear flow with velocity V1 and density ρ1 in the upper half
plane and V2 and ρ2 in the lower half plane. We expect instability would
occur within crossing time scale of the flow over the characteristic length
scale. The surplus kinetic and potential energies proportional to
(V2 − V1)2
and ρ1 − ρ2 are the energy sources of turbulence.
The hydrodynamical equations are
∂ρ
∂t
+ div(ρv) = 0,
ρ
∂v
∂t
+ρv·∇v−ρg = −∇p = −a2
∇ρ.
We will assume an initial state with
v = (V (y), 0, 0) and ρ = ρ0(y).
V1
V2
y
x
2
1
ρ
ρ
g
29
K-H & R-T instabilities: perturbing the initial state
We will assume the perturbed quantities in the form of
v = (V (y) + δ˜vx , δ˜vy , 0) and ρ = ρ0(y) + δ˜ρ, where δ˜vx, y ≪ V (y) and
δ˜ρ ≪ ρ0. The perturbations are assumed to obey harmonical expansion
δ˜vx,y = δvx,y (y) exp(ikx + iωt),
δ˜ρ = δρ(y) exp(ikx + iωt).
Therefore, we shall substitute ∂/∂t → iω and ∇ → (ik, ∂/∂y, 0).
The linearized hydrodynamical equations are
∂ρ
∂t
+ div(ρv) = 0 → ωδρ + Vkδρ + ρ0kδvx − iρ0δv′
y = 0,
ρ
∂vx
∂t
+ρvi
∂vx
∂xi
= −a2 ∂ρ
∂x
→ ρ0ωδvx +ρ0Vkδvx −iρ0δvy V ′
= −a2
kδρ,
ρ
∂vy
∂t
+ρvi
∂vy
∂xi
= −a2 ∂ρ
∂y
−ρg → ρ0ωδvy +ρ0Vkδvy = ia2
δρ′
+iδρg,
where prime (′
) denotes d/dy.
30
K-H & R-T instabilities: solving for perturbations
Solving the second equation for δvx , inserting to the first one and solving
for δρ, and inserting the final relation to the last equation we derive
(denoting ωd = ω + kV )
ρ0ωdδvy = a2
−ρ0δv′
y ωd + ρ0kV ′
δvy
ω2
d − k2a2
′
+ g
−ρ0δv′
y ωd + ρ0kV ′
δvy
ω2
d − k2a2
.
For relatively slow flow (V ≪ a) we can assume a → ∞ (incompresible
flow) and the dispersion relation becomes
ρ0ωdδv′
y − ρ0kV ′
δvy
′
− ρ0ωdk2
δvy = 0.
31
K-H & R-T instabilities: back to the original problem
For an assumed velocity and density profile profile
V (y) =
V1, for y > 0,
V2, for y < 0,
ρ0(y) =
ρ1, for y > 0,
ρ2, for y < 0,
in each half-space. We have the dispersion relation δv′′
y − k2
δvy = 0 that
has the solution (assuming δvy → 0 for y → ∞)
δvy ∼
exp(−ky), for y > 0,
exp(ky), for y < 0.
The displacement δy at the boundary should be continuous, consequently
Dδy
Dt
=
∂
∂t
+ V
∂
∂x
δy = δvy
and therefore vy /(ω + kV ) shoud be continuous. The solution becomes
δvy ∼
(ω + kV1) exp(−ky), for y > 0,
(ω + kV2) exp(ky), for y < 0.
32
Kelvin-Helmholtz instability
We assume constant density, no gravitational field, and velocity shear
V (y) =
V1, for y > 0,
V2, for y < 0.
From the requirement that the left-hand side of the dispersion relation
ρ0ωdδv′
y − ρ0kV ′
δvy
′
= ρ0ωdk2
δvy
should be continuous at the boundary we have
ρ0ωdδv′
y 1
= ρ0ωdδv′
y 2
,
which, after substitution of ωd and δv′
y gives the dispersion relation
(ω + kV1)2
+ (ω + kV2)2
= 0.
Solving for ω gives instability for V1 = V2:
ω = −
1
2
k(V1 + V2) ±
1
2
ik(V1 − V2).
33
Kelvin-Helmholtz instability: going nonlinear
34
Kelvin-Helmholtz instability: Earth’s atmosphere
35
Kelvin-Helmholtz instability: Earth’s atmosphere
36
Kelvin-Helmholtz instability: Solar prominence
37
Kelvin-Helmholtz instability: Atmosphere of Saturn
38
Rayleigh-Taylor instability
We assume zero velocity V (y) and the density
ρ0(y) =
ρ1, for y > 0,
ρ2, for y < 0.
From the requirement that the left-hand side of the dispersion relation
(assuming ω ≫ k2
a2
)
ρ0ω2
δvy + gρ0δv′
y = a2
−ρ0δv′
y
ωd
′
should be continuous at the boundary we have
ρ0ω2
δvy + gρ0δv′
y 1
= ρ0ω2
δvy + gρ0δv′
y 2
.
Inserting the solution δvy ∼ exp(±ky) gives the dispersion relation
ω2
= gk
ρ2 − ρ1
ρ2 + ρ1
.
The flow is stable (ω2
> 0) for ρ2 > ρ1, while for ρ2 < ρ1 the
Rayleigh-Taylor instability appears.
39
Fingers of Rayleigh-Taylor instability
Figures show the developement of the finger typical for Rayleigh-Taylor
instability. The instability is stabilized by the surface tension for large
wavenumbers (Chandrasekhar). Figure shows also Kelvin-Helmholtz
instabilities on the boundary of finger. 40
Visualisation of the Rayleigh-Taylor instability
41
Crab nebula
42
Crab nebula structure due to Rayleigh-Taylor instability
While the supernova nebula becomes flat, the swept-up, accelerating
shell is subject to the Rayleigh–Taylor instability (Kulsrud et al. 1965,
Chevalier & Gull 1975, Blondin & Chevalier 2017).
43
Suggested reading
S. Chandrasekhar: Hydrodynamic and Hydromagnetic Stability
R. Courant & K. O. Friedrichs: Supersonic Flow and Shock Waves
A. Feldmeier: Theoretical Fluid Dynamics
A. Maeder: Physics, Formation and Evolution of Rotating Stars
D. Mihalas & B. W. Mihalas: Foundations of Radiation Hydrodynamics
F. H. Shu: The physics of astrophysics: II. Hydrodynamics
T. Tajima & K. Shibata: Plasma astrophysics
L. C. Woods: Physics of Plasmas
44