Astro320 Notes
Astro320 Notes
2
Course Outline
Astronomy is the study of everything in our Universe outside of the Earth’s atmo-
sphere, including the Universe as a whole (‘cosmology’). As such, it involves all
of physics: quantum mechanics, nuclear and particle physics, classical mechanics,
special and general relativity, electromagnetism, statistical physics, hydrodynamics,
plasma physics, and even solid state physics. As we will see, though, with the excep-
tion of rocky planets and asteroids, all objects in the Universe can be characterized
as some kind of fluid. In addition, almost all information we receive from these
objects reaches us in the form of radiation (photons)1 . Consequently, in this course
on physical processes in astronomy we will focus almost exclusively on fluid dy-
namics and radiative processes. Note that we focus exclusively on neutral fluids;
the physics of electrically charged fluids, known as plasmas will not be covered.
It is assumed that the student is familiar with vector calculus, with curvi-linear
coordinate systems. A brief overview of these topics is provided in Appendices A-E.
The other appendices present detailed background information that isprovided for
the interested student, but which is not considered part of the course material.
1
Other astrophysical messengers include neutrinos, cosmic rays and gravitational waves
3
CONTENTS
Part I: Fluid Dynamics
1: Introduction to Fluids and Plasmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2: Dynamical Treatments of Fluids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3: Hydrodynamic Equations for Ideal Fluids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4: Viscosity, Conductivity & The Stress Tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5: Hydrodynamic Equations for Non-Ideal Fluids . . . . . . . . . . . . . . . . . . . . . . . . . 34
6: Equations of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
7: Vorticity & Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
8; Hydrostatics and Steady Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .53
9: Viscous Flow and Accretion Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
10: Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .77
11: Sound Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
12: Shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
13: Fluid Instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4
APPENDICES
5
LITERATURE
The material covered and presented in these lecture notes has relied heavily on a
number of excellent textbooks listed below.
• Theoretical Astrophysics
by M. Bartelmann (ISBN-978-3-527-41004-0)
• Galactic Dynamics
by J. Binney & S. Tremaine (ISBN-978-0-691-13027-9)
6
7
Part I: Fluid Dynamics
Almost everything we encounter in the Universe, from gas planets to stars, and
from the interstellar medium to galaxies, can be categorized as some kind of fluid.
Hence, understanding astrophysical processes requires a solid understanding of fluid
dynamics. The following chapters present fairly detailed description of the dynamics
of fluids with an application to astrophysics.
Fluid dynamics is a rich topic, and one could easily devote an entire course to it.
The following chapters therefore only scratch the surface of this rich topic. Readers
who want to get more indepth information are referred to the following excellent
textbooks
- The Physics of Fluids and Plasmas by A. Choudhuri
- Modern Fluid Dynamics for Physics and Astrophysics by [Link] et al.
- The Physics of Astrophysics II. Gas Dynamics by F. Shu
- Principles of Astrophysical Fluid Dynamics by C. Clarke & B. Carswell
- Modern Classical Physics by [Link] & R. Blandford
8
CHAPTER 1
What is a fluid?
A fluid is a substance that can flow, has no fixed shape, and offers little resistance
to an external stress
• In a fluid the constituent particles (atoms, ions, molecules, stars) can ‘freely’
move past one another.
• A fluid changes its shape at a steady rate when acted upon by a stress force.
What is a plasma?
A plasma is a fluid in which (some of) the consistituent particles are electrically
charged, such that the interparticle force (Coulomb force) is long-range in nature.
Fluid Demographics:
All fluids are made up of large numbers of constituent particles, which can be
molecules, atoms, ions, dark matter particles or even stars. Different types of fluids
mainly differ in the nature of their interparticle forces. Examples of inter-particle
forces are the Coulomb force (among charged particles in a plasma), vanderWaals
forces (among molecules in a neutral fluid) and gravity (among the stars in a galaxy).
Fluids can be both collisional or collisionless, where we define a collision as an
interaction between constituent particles that causes the trajectory of at least one
of these particles to be deflected ‘noticeably’. Collisions among particles drive the
system towards thermodynamic equilibrium (at least locally) and the velocity
distribution towards a Maxwell-Boltzmann distribution.
In neutral fluids the particles only interact with each other on very small scales.
Typically the inter-particle force is a vanderWaals force, which drops off very rapidly.
Put differently, the typical cross section for interactions is the size of the particles
(i.e., the Bohr radius for atoms), which is very small. Hence, to good approximation
9
Figure 1: Examples of particle trajectories in (a) a collisional, neutral fluid, (b) a
plasma, and (c) a self-gravitating collisionless, neutral fluid. Note how different the
dynamics are.
In a fully ionized plasma the particles exert Coulomb forces (F~ ∝ r −2 ) on each other.
Because these are long-range forces, the velocity of a charged particle changes more
likely due to a succession of many small deflections rather than due to one large one.
As a consequence, particles trajectories in a highly ionized plasma (see Fig. 1b) are
very different from those in a neutral fluid.
10
system. We can then write that
Here hF~ ii is the time (or ensemble) averaged force at the instantaneous position of
particle i and δ F~i (t) is the instantaneous deviation due to the discrete nature of the
particles that make up the system. As N → ∞ then δ F~i → 0 and the system is
said to be collisionless; its dynamics are governed by the collective force from all
particles rather than by collisions/interactions with individual particles.
11
charges is screened beyond the Debye length:
1/2
kB T
λD = ≃ 4.9 cm n−1/2 T 1/2
8π n e2
Here n is the number density in cm−3 , T is the temperature in degrees Kelvin, and
e is the electrical charge of an electron in e.s.u. Related to the Debye length is the
Plasma parameter
1
g≡ ≃ 8.6 × 10−3 n1/2 T −3/2
n λ3D
As an example, let’s consider three different astrophysical plasmas: the ISM (inter-
stellar medium), the ICM (intra-cluster medium), and the interior of the Sun. The
warm phase of the ISM has a temperature of T ∼ 104 K and a number density of
n ∼ 1 cm−3 . This implies ND ∼ 1.2 × 108 . Hence, the warm phase of the ISM can be
treated as a collisionless plasma on sufficiently small time-scales (for example when
treating high-frequency plasma waves). The ICM has a much lower average density
of ∼ 10−4 cm−3 and a much higher temperature (∼ 107 K). This implies a much
larger number of particles per Debye volume of ND ∼ 4 × 1014 . Hence, the ICM can
typically be approximated as a collisionless plasma. The interior of stars, though,
has a similar temperature of ∼ 107 K but at much higher density (n ∼ 1023 cm−3 ),
implying ND ∼ 10. Hence, stellar interiors are highly collisional plasmas!
12
Compressibility: Fluids and plasmas can be either gaseous or liquid. A gas is
compressible and will completely fill the volume available to it. A liquid, on the
other hand, is (to good approximation) incompressible, which means that a liquid
of given mass occupies a given volume.
NOTE: Although a gas is said to be compressible, many gaseous flows (and virtually
all astrophysical flows) are incompressible. When the gas is in a container, you can
easily compress it with a piston, but if I move my hand (sub-sonically) through the
air, the gas adjust itself to the perturbation in an incompressible fashion (it moves out
of the way at the speed of sound). The small compression at my hand propagates
forward at the speed of sound (sound wave) and disperses the gas particles out of
the way. In astrophysics we rarely encounter containers, and subsonic gas flow is
often treated (to good approximation) as being incompressible.
Throughout what follows, we use ‘fluid’ to mean a neutral fluid, and ‘plasma’ to
refer to a fluid in which the particles are electrically charged.
NOTE: An ideal (or perfect) fluid should NOT be confused with an ideal or perfect
gas, which is defined as a gas in which the pressure is solely due to the kinetic motions
of the constituent particles. As we show in Chapter 6, and as you have probably
seen before, this implies that the pressure can be written as P = n kB T , with n the
particle number density, kB the Boltzmann constant, and T the temperature.
13
star cover many orders of magnitude. To good approximation, its equation
of state is that of an ideal gas.
• Giant (gaseous) planets: Similar to stars, gaseous planets are large spheres
of gas, albeit with a rocky core. Contrary to stars, though, the gas is typically
so dense and cold that it can no longer be described with the equation of
state of an ideal gas.
• White Dwarfs & Neutron stars: These objects (stellar remnants) can be
described as fluids with a degenerate equation of state.
• Proto-planetary disks: the dense disks of gas and dust surrounding newly
formed stars out of which planetary systems form.
14
• Accretion disks: Accretion disks are gaseous, viscous disks in which the
viscosity (enhanced due to turbulence) causes a net rate of radial infall towards
the center of the disk, while angular momentum is being transported outwards
(accretion)
15
CHAPTER 2
2. a (set of) equation(s) to describe how the state variables change with time
This system is described by the position vectors (~x) and momentum vectors (~p) of
all the N particles, i.e., by (~x1 , ~x2 , ..., ~xN , ~p1 , ~p2 , ..., ~pN ).
If the particles are trully classical, in that they can’t emit or absorb radiation, then
one can define a Hamiltonian
N
X
H(~xi , ~pi , t) ≡ H(~x1 , ~x2 , ..., ~xN , ~p1 , p~2 , ..., p~N , t) = ~pi · ~x˙ i − L(~xi , ~x˙ i , t)
i=1
where L(~xi , ~x˙ i , t) is the system’s Lagrangian, and ~x˙ i = d~xi /dt.
The equations that describe the time-evolution of these state-variables are the Hamil-
tonian equations of motion:
∂H ∂H
~x˙ i = ; p~˙ i = −
∂~pi ∂~xi
16
Example 2: an electromagnetic field
~ x) and
The state of this system is described by the electrical and magnetic fields, E(~
~ x), respectively, and the equations that describe their evolution with time are the
B(~
Maxwell equations, which contain the terms ∂ E/∂t ~ ~
and ∂ B/∂t.
17
h h
λ= ≃√
p mkB T
Here h is the Planck constant, p is the particle’s momentum, m is the particle mass,
kB is the Boltzman constant, and T is the temperature of the fluid. This de Broglie
wavelength indicates the ‘characteristic’ size of the wave-packet that according to
quantum mechanics describes the particle, and is typically very small. Except for
extremely dense fluids such as white dwarfs and neutron stars, or ‘exotic’ types of
dark matter (i.e., ‘fuzzy dark matter’), the de Broglie wavelength is always much
smaller than the mean particle separation, and classical, Newtonian mechanics suf-
fices. As we have seen above, a classical, Newtonian system of N particles can be
described by a Hamiltonian, and the corresponding equations of motions. We refer
to this as the level-1 description of fluid dynamics (see under ‘example 1’ above).
Clearly, when N is very large is it unfeasible to solve the 2N equations of motion for
all the positions and momenta of all particles. We need another approach.
In the level-2 approach, one introduces the distribution function f (~x, p~, t), which
describes the number density of particles in 6-dimensional ‘phase-space’ (~x, p~) (i.e.,
how many particles are there with positions in the 3D volume ~x +d~x and momenta in
the 3D volume ~p + d~p). The equation that describes how f (~x, p~, t) evolves with time
is called the Boltzmann equation for a neutral fluid. If the fluid is collisionless
this reduces to the Collisionless Boltzmann equation (CBE). If the collisionless
fluid is a plasma, the same equation is called the Vlasov equation. Often the CBE
and the Vlasov equation are used without distinction.
At the final level-3, the fluid is modelled as a continuum. This means we ignore
that fluids are made up of constituent particles, and rather describe the fluid with
continuous fields, such as the density and velocity fields ρ(~x) and ~u(~x) which assign
to each point in space a scalar quantity ρ and a vector quantity ~u, respectively. For
an ideal neutral fluid, the state in this level-3 approach is fully described by four
fields: the density ρ(~x), the velocity field ~u(~x), the pressure P (~x), and the internal,
specific energy ε(~x) (or, equivalently, the temperature T (~x)). In the MHD treat-
ment of plasmas one also needs to specify the magnetic field B(~ ~ x). The equations
that describe the time-evolution of ρ(~x), ~u(~x), and ε(~x) are called the continuity
equation, the Navier-Stokes equations, and the energy equation, respectively.
Collectively, we shall refer to these as the hydrodynamic equations or fluid equa-
tions. In MHD you have to slightly modify the Navier-Stokes equations, and add
an additional induction equation describing the time-evolution of the magnetic
18
field. For an ideal (or perfect) fluid (i.e., no viscosity and/or conductivity), the
Navier-Stokes equations reduce to what are known as the Euler equations. For a
collisionless gravitational system, the equivalent of the Euler equations are called the
Jeans equations.
Throughout this course, we mainly focus on the level-3 treatment, to which we refer
hereafter as the macroscopic approach. However, for completeness we will derive
these continuum equations starting from a completely general, microscopic level-1
treatment. Along the way we will see how subtle differences in the inter-particle forces
gives rise to a rich variety in dynamics (fluid vs. plasma, collisional vs. collisionless).
1. the FE needs to be much smaller than the characteristic scale in the problem,
which is the scale over which the hydrodynamical quantities Q change by an
order of magnitude, i.e.
Q
lFE ≪ lscale ∼
∇Q
2. the FE needs to be sufficiently large that fluctuations due to the finite number
of particles (‘discreteness noise’) can be neglected, i.e.,
3
n lFE ≫1
3. the FE needs to be sufficiently large that it ‘knows’ about the local conditions
through collisions among the constituent particles, i.e.,
lFE ≫ λ
19
The ratio of the mean-free path, λ, to the characteristic scale, lscale is known as the
Knudsen number: Kn = λ/lscale . Fluids typically have Kn ≪ 1; if not, then one
is not justified in using the continuum approach (level-3) to fluid dynamics, and one
is forced to resort to a more statistical approach (level-2).
Note that fluid elements can NOT be defined for a collisionless fluid (which has
an infinite mean-free path). This is one of the reasons why one cannot use the
macroscopic approach to derive the equations that govern a collisionless fluid.
In the case of an ideal (or perfect) fluid (i.e., with zero viscosity and conductivity),
the Navier-Stokes equations (which are the hydrodynamical momentum equations)
reduce to what are called the Euler equations. In that case, the evolution of fluid
elements is describe by the following set of hydrodynamical equations:
20
• If the EoS is barotropic, i.e., if P = P (ρ), then the energy equation is not needed
to close the set of equations. There are two barotropic EoS that are encountered
frequently in astrophysics: the isothermal EoS, which describes a fluid for which
cooling and heating always balance each other to maintain a constant temperature,
and the adiabatic EoS, in which there is no net heating or cooling (other than
adiabatic heating or cooling due to the compression or expansion of volume, i.e., the
P dV work). We will discuss these cases in more detail later in the course.
• No EoS exists for a collisionless fluid. Consequently, for a collisionless fluid one
can never close the set of fluid equations, unless one makes a number of simplifying
assumptions (i.e., one postulates various symmetries)
• If the fluid is not ideal, then the momentum equations include terms that contain
the (kinetic) viscosity, ν, and the energy equation includes a term that contains
the conductivity, K. Both ν and K depend on the mean-free path of the constituent
particles and therefore depend on the temperature and collisional cross-section of
the particles. Closure of the set of hydrodynamic equations then demands additional
constitutive equations ν(T ) and K(T ). Often, though, ν and K are simply assumed
to be constant (the T -dependence is ignored).
• If the fluid is self-gravitating (which is the case, for example, for stars and
galaxies) there is an additional unknown, the gravitational potential Φ. However,
there is also an additional equation, the Poisson equation relating Φ to ρ, so that
the set of equations remains closed.
• In the case of a plasma, the charged particles give rise to electric and magnetic
fields. Each fluid element now carries 6 additional scalars (Ex , Ey , Ez , Bx , By , Bz ),
and the set of equations has to be complemented with the Maxwell equations that
describe the time evolution of E~ and B.
~
21
Fluid Dynamics: Eulerian vs. Lagrangian Formalism:
One distinguishes two different formalisms for treating fluid dynamics:
• Eulerian Formalism: in this formalism one solves the fluid equations ‘at
fixed positions’: the evolution of a quantity Q is described by the local (or
partial, or Eulerian) derivative ∂Q/∂t. An Eulerian hydrodynamics code is a
‘grid-based code’, which solves the hydro equations on a fixed grid, or using
an adaptive grid, which refines resolution where needed. The latter is called
Adaptive Mesh Refinement (AMR).
dQ ∂Q
= + ~u · ∇Q
dt ∂t
~ x, t), it is straightforward
Using a similar derivation, but now for a vector quantity A(~
to show that
22
~
dA ~
∂A
= ~
+ (~u · ∇) A
dt ∂t
which, in index-notation, is written as
dAi ∂Ai ∂Ai
= + uj
dt ∂t ∂xj
Another way to derive the above relation between the Eulerian and Lagrangian
derivatives, is to think of dQ/dt as
dQ Q(~x + δ~x, t + δt) − Q(~x, t)
= lim
dt δt→0 δt
Using that
~x(t + δt) − ~x(t) δ~x
~u = lim =
δt→0 δt δt
and
Q(~x + δ~x, t) − Q(~x, t)
∇Q = lim
x→0
δ~ δ~x
it is straightforward to show that this results in the same expression for the substan-
tial derivative as above.
• Streaklines: the locus of points of all the fluid particles that have passed con-
tinuously through a particular spatial point in the past. Dye steadily injected
into the fluid at a fixed point extends along a streakline.
23
Figure 2: Streaklines showing laminar flow across an airfoil; made by injecting dye
at regular intervals in the flow
• Particle paths: (aka pathlines) are the trajectories that individual fluid ele-
ments follow. The direction the path takes is determined by the streamlines of
the fluid at each moment in time.
Only if the flow is steady, which means that all partial time derivatives (i.e., ∂~u/∂t =
∂ρ/∂t = ∂P/∂t) vanish, will streamlines be identical to streaklines be identical to
particle paths. For a non-steady flow, they will differ from each other.
24
CHAPTER 3
Without any formal derivation (this comes later) we now present the hydrodynamic
equations for an ideal, neutral fluid. Note that these equations adopt the level-3
continuum approach discussed in the previous chapter.
Lagrangian Eulerian
dρ ∂ρ
Continuity Eq: = −ρ ∇ · ~u + ∇ · (ρ~u) = 0
dt ∂t
d~u ∇P ∂~u ∇P
Momentum Eqs: =− − ∇Φ + (~u · ∇) ~u = − − ∇Φ
dt ρ ∂t ρ
dε P L ∂ε P L
Energy Eq: = − ∇ · ~u − + ~u · ∇ε = − ∇ · ~u −
dt ρ ρ ∂t ρ ρ
NOTE: students should become familiar with switching between the Eulerian and
Lagrangian equations, and between the vector notation shown above and the
index notation. The latter is often easier to work with. When writing down the
index versions, make sure that each term carries the same index, and make use of the
Einstein summation convention. The only somewhat tricky term is the (~u · ∇) ~u-term
in the Eulerian momentum equations, which in index form is given by uj (∂ui /∂xj ),
where i is the index carried by each term of the equation.
25
(∇ · ~u = 0), which is also called solenoidal.
Momentum Equations: these equations simply state than one can accelerate a
fluid element with either a gradient in the pressure, P , or a gradient in the gravita-
tional potential, Φ. Basically these momentum equations are nothing but Newton’s
F~ = m~a applied to a fluid element. In the above form, valid for an inviscid, ideal
fluid, the momentum equations are called the Euler equations.
Energy Equation: the energy equation states that the only way that the specific,
internal energy, ε, of a fluid element can change, in the absence of conduction, is
by adiabatic compression or expansion, which requires a non-zero divergence of the
velocity field (i.e., ∇ · ~u 6= 0), or by radiation (emission or absorption of photons).
The latter is expressed via the net volumetric cooling rate,
dQ
L=ρ =C−H
dt
Here Q is the thermodynamic heat, and C and H are the net volumetric cooling and
heating rates, respectively.
If the ideal fluid is governed by self-gravity (as opposed to, is placed in an external
gravitational field), then one needs to complement the hydrodynamical equations
with the Poisson equation: ∇2 Φ = 4πGρ. In addition, closure requires an addi-
tional constitutive relations in the form of an equation-of-state P = P (ρ, ε). If
the ideal fluid obeys the ideal gas law, then we have the following two constitutive
relations:
kB T 1 kB T
P = ρ, ε=
µ mp γ − 1 µ mp
(see Chapter 6 for details). Here µ is the mean molecular weight of the fluid in units
of the proton mass, mp , and γ is the adiabatic index, which is often taken to be 5/3
as appropriate for a mono-atomic gas.
∂ρ
+ ∇ · (ρ~u) = 0
∂t
∂ρ~u
+ ∇ · Π = −ρ∇Φ
∂t
∂E ∂Φ
+ ∇ · [(E + P ) ~u] = ρ −L
∂t ∂t
Here
Π = ρ ~u ⊗ ~u + P
is the momentum flux density tensor (of rank 2), and
1 2
E=ρ u +Φ+ε
2
NOTE: In the expression for the momentum flux density tensor A⊗ ~ B ~ is the tensor
product of A ~ and B~ defined such that (A
~ ⊗ B)
~ ij = ai bj (see Appendix A). Hence,
the index-form of the momentum flux density tensor is simply Πij = ρ ui uj + P δij ,
with δij the Kronecker delta function. Note that this expression is ONLY valid for
an ideal fluid; in the next chapter we shall derive a more general expression for the
momentum flux density tensor.
Note also that whereas there is no source or sink term for the density, gradients in the
gravitational field act as a source of momentum, while its time-variability can cause
an increase or decrease in the energy density of the fluid (if the fluid is collisionless,
we call this violent relaxation). Another source/sink term for the energy density
is radiation (emission or absorption of photons).
27
CHAPTER 4
The hydrodynamic equations presented in the previous chapter are only valid for
an ideal fluid, i.e., a fluid without viscosity and conduction. We now examine the
origin of conduction and viscosity, and link the latter to the stress tensor, which
is an important quantity in all of fluid dynamics.
In an ideal fluid, the particles effectively have a mean-free path of zero, such that they
cannot communicate with their neighboring particles. In reality, though, the mean-
free path, λmfp = (nσ)−1 is finite, and particles ”communicate” with each other
through collisions. These collisions cause an exchange of momentum and energy
among the particles involved, acting as a relaxation mechanism. Note that in a
collisionless system the mean-free path is effectively infinite, and there is no two-
body relaxation, only collective relaxation mechanisms (i.e., violent relaxation or
wave-particle interactions).
• When there are gradients in velocity (”shear”) then the collisions among neigh-
boring fluid elements give rise to a net transport of momentum. The collisions
drive the system towards equilibrium, i.e., towards no shear. Hence, the collisions
act as a resistance to shear, which is called viscosity. See Fig. 3 for an illustration.
• When there are gradients in temperature (or, in other words, in specific inter-
nal energy), then the collisions give rise to a net transport of energy. Again,
the collisions drive the system towards equilibrium, in which the gradients vanish,
and the rate at which the fluid can erase a non-zero ∇T is called the (thermal)
conductivity.
28
Figure 3: Illustration of origin of viscosity and shear stress. Three neighboring fluids
elements (1, 2 and 3) have different streaming velocities, ~u. Due to the microscopic
motions and collisions (characterized by a non-zero mean free path), there is a net
transfer of momentum from the faster moving fluid elements to the slower moving
fluid elements. This net transfer of momentum will tend to erase the shear in ~u(~x),
and therefore manifests itself as a shear-resistance, known as viscosity. Due to the
transfer of momentum, the fluid elements deform; in our figure, 1 transfers linear
momentum to the top of 2, while 3 extracts linear momentum from the bottom of 2.
Consequently, fluid element 2 is sheared as depicted in the figure at time t+∆t. From
the perspective of fluid element 2, some internal force (from within its boundaries)
has exerted a shear-stress on its bounding surface.
29
Non-uniform Gases” by S. Chapman and T. Cowling. Using the Chapman-Enskog
expansion one finds the following expressions for µ and K:
1/2
a m kB T 5
µ= , K = cV µ
σ π 2
Here a is a numerical factor that depends on the details of the interparticle forces, σ
is the collisional cross section, and cV is the specific heat (i.e., per unit mass). Thus,
for a given fluid (given σ and m) we basically have that µ = µ(T ) and K = K(T ).
Note that µ ∝ T 1/2 ; viscosity increases with temperature. This only holds for gases!
For liquids we know from experience that viscosity decreases with increasing tem-
perature (think of honey). Since in astrophysics we are mainly concerned with gas,
µ ∝ T 1/2 will be a good approximation for most of what follows.
Now that we have a rough idea of what viscosity (resistance to shear) and conduc-
tivity (resistance to temperature gradients) are, we have to ask how to incorporate
them into our hydrodynamic equations.
~v = ~u + w
~
where h~v i = ~u, hwi
~ = 0 and h.i indicates the average over a fluid element. If we
define vi as the velocity in the i-direction, we have that
hvi vj i = ui uj + hwi wj i
These different velocities allow us to define a number of different velocity tensors:
Stress Tensor: σij ≡ −ρhwi wj i σ = −ρw
~ ⊗w ~
Momentum Flux Density Tensor: Πij ≡ +ρhvi vj i Π = +ρ~v ⊗ ~v
Ram Pressure Tensor: Σij ≡ +ρui uj Σ = +ρ~u ⊗ ~u
30
which are related according to σ = Σ − Π. Note that each of these tensors is man-
ifest symmetric (i.e., σij = σji , etc.), which implies that they have 6 independent
variables.
Note that the stress tensor is related to the microscopic random motions. These
are the ones that give rise to pressure, viscosity and conductivity! The reason that
~ x, n̂) acting on a
σij is called the stress tensor is that it is related to the stress Σ(~
surface with normal vector n̂ located at ~x according to
Σi (n̂) = σij nj
Here Σi (n̂) is the i-component of the stress acting on a surface with normal n̂, whose
j-component is given by nj . Hence, in general the stress will not necessarily be along
the normal to the surface, and it is useful to decompose the stress in a normal
stress, which is the component of the stress along the normal to the surface, and a
shear stress, which is the component along the tangent to the surface.
To see that fluid elements in general are subjected to shear stress, consider the
following: Consider a flow (i.e., a river) in which we inject a small, spherical blob (a
fluid element) of dye. If the only stress to which the blob is subject is normal stress,
the only thing that can happen to the blob is an overall compression or expansion.
However, from experience we know that the blob of dye will shear into an extended,
‘spaghetti’-like feature; hence, the blob is clearly subjected to shear stress, and this
shear stress is obvisouly related to another tensor called the deformation tensor
∂ui
Tij =
∂xj
Since ∂ui /∂xj = 0 in a static fluid (~u(~x) = 0), we see that in a static fluid the stress
tensor can only depend on the normal stress, which we call the pressure.
31
The minus sign is a consequence of the sign convention of the stress.
Sign Convention: The stress Σ(~ ~ x, n̂) acting at location ~x on a surface with normal
n̂, is exerted by the fluid on the side of the surface to which the normal points, on
the fluid from which the normal points. In other words, a positive stress results in
compression. Hence, in the case of pure, normal pressure, we have that Σ = −P .
Viscous Stress Tensor: The expression for the stress tensor in the case of static
fluid motivates us to write in general
where we have introduced a new tensor, τij , which is known as the viscous stress
tensor, or the deviatoric stress tensor.
Since the deviatoric stress tensor, τij , is only non-zero in the presence of shear in the
fluid flow, this suggests that
∂uk
τij = Tijkl
∂xl
where Tijkl is a proportionality tensor of rank four. As described in Appendix F
(which is NOT part of the curriculum for this course), most (astrophysical) fluids
are Newtonian, in that they obey a number of conditions. As detailed in that
appendix, for a Newtonian fluid, the relation between the stress tensor and the
deformation tensor is given by
∂ui ∂uj 2 ∂uk ∂uk
σij = −P δij + µ + − δij + η δij
∂xj ∂xi 3 ∂xk ∂xk
Here P is the pressure, δij is the Kronecker delta function, µ is the coefficient
of shear viscosity, and η is the coefficient of bulk viscosity (aka the ‘second
viscosity’). We thus see that for a Newtonian fluid, the stress tensor, despite being
a symmetric tensor of rank two (which implies 6 independent variables), only has
three independent components: P , µ and η.
Let’s take a closer look at these three quantities, starting with the pressure P . To be
exact, P is the thermodynamic equilibrium pressure, and is normally computed
thermodynamically from some equation of state, P = P (ρ, T ). It is related to the
32
translational kinetic energy of the particles when the fluid, in equilibrium, has reached
equipartition of energy among all its degrees of freedom, including (in the case of
molecules) rotational and vibrations degrees of freedom.
Using the above expression for σij , and using that ∂uk /∂xk = ∇ · ~u (Einstein sum-
mation convention), it is easy to see that
Pm = P − η ∇ · ~u
From this expression it is clear that the bulk viscosity, η, is only non-zero if P 6=
Pm . This, in turn, can only happen if the constituent particles of the fluid have
degrees of freedom beyond position and momentum (i.e., when they are molecules
with rotational or vibrational degrees of freedom). Hence, for a fluid of monoatoms
(ideal gas), η = 0. From the fact that P = Pm + η∇ · ~u it is clear that for an
incompressible flow P = Pm and the value of η is irrelevant; bulk viscosity plays
no role in incompressible fluids or flows. The only time when Pm 6= P is when a
fluid consisting of particles with internal degrees of freedom (e.g., molecules) has just
undergone a large volumetric change (i.e., during a shock). In that case there may
be a lag between the time the translational motions reach equilibrium and the time
when the system reaches full equipartition in energy among all degrees of freedom.
In astrophysics, bulk viscosity can generally be ignored, but be aware that it may
be important in shocks. This only leaves the shear viscosity µ, which describes the
ability of the fluid to resist shear stress via momentum transport resulting from
collisions and the non-zero mean free path of the particles.
33
CHAPTER 5
As we have seen in the previous chapter, the effect of viscosity is captured by the
stress tensor, which is given by
∂ui ∂uj 2 ∂uk ∂uk
σij = −P δij + τij = −P δij + µ + − δij + η δij
∂xj ∂xi 3 ∂xk ∂xk
Note that in the limit µ → 0 and η → 0, valid for an ideal fluid, σij = −P δij . This
suggests that we can incorporate viscosity in the hydrodynamic equations by simply
replacing the pressure P with the stress tensor, i.e., P δij → −σij = P δij − τij .
It is more common, and more useful, to write out the viscous stress tensor, yielding
dui ∂P ∂ ∂ui ∂uj 2 ∂uk ∂ ∂uk ∂Φ
ρ =− + µ + − δij + η −ρ
dt ∂xi ∂xj ∂xj ∂xi 3 ∂xk ∂xi ∂xk ∂xi
34
These are the Navier-Stokes equations (in Lagragian index form) in all their glory,
containing both the shear viscosity term and the bulk viscosity term (the latter
is often ignored).
Note that µ and η are usually functions of density and temperature so that they
have spatial variations. However, it is common to assume that these are suficiently
small so that µ and η can be treated as constants, in which case they can be taken
outside the differentials. In what follows we will make this assumption as well.
where we have introduced the kinetic viscosity ν ≡ µ/ρ. Note that these equations
reduce to the Euler equations in the limit ν → 0. Also, note that the ∇(∇ · ~u)
term is only significant in the case of flows with variable compression (i.e., viscous
dissipation of accoustic waves or shocks), and can often be ignored. This leaves the
ν∇2~u term as the main addition to the Euler equations. Yet, this simple ‘diffuse’
term (describing viscous momentum diffusion) dramatically changes the charac-
ter of the equation, as it introduces a higher spatial derivative. Hence, additional
boundary conditions are required to solve the equations. When solving problems
with solid boundaries (not common in astrophysics), this condition is typically that
the tangential (or shear) velocity at the boundary vanishes. Although this may sound
ad hoc, it is supported by observation; for example, the blades of a fan collect dust.
Recall that when writing the Navier-Stokes equation in Eulerian form, we have that
d~u/dt → ∂~u/∂t + ~u · ∇~u. It is often useful to rewrite this extra term using the vector
calculus identity
~u · ~u
~u · ∇~u = ∇ + (∇ × ~u) × ~u
2
35
Hence, for an irrotational flow (i.e., a flow for which ∇ × ~u = 0), we have that
~u · ∇~u = 12 ∇u2 , where u ≡ |~u|.
Now that we have added the effect of viscosity, what remains is to add conduction.
We can make progress by realizing that, on the microscopic level, conduction arises
from collisions among the constituent particles, causing a flux in internal energy. The
internal energy density of a fluid element is h 21 ρw 2 i, where w
~ = ~v − ~u is the random
motion of the particle wrt the fluid element (see Chapter 4), and the angle brackets
indicate an ensemble average over the particles that make up the fluid element. Based
on this we see that the conductive flux in the i-direction can be written as
1
Fcond,i = h ρw 2 wi i = hρεwi i
2
From experience we also know that we can write the conductive flux as
F~cond = −K ∇T
36
with K the thermal conductivity.
Next we realize that conduction only causes a net change in the internal energy at
some fixed position if the divergence in the conductive flux (∇· F~cond ) at that position
is non-zero. This suggests that the final form of the energy equation, for a non-ideal
fluid, and in Lagrangian vector form, has to be
dε
ρ = −P ∇ · ~u − ∇ · F~cond + V − L
dt
To summarize, below we list the full set of equations of gravitational, radial hydro-
dynamics (ignoring bulk viscosity) 2 .
dρ
Continuity Eq. = −ρ ∇ · ~u
dt
d~u 2 1
Momentum Eqs. ρ = −∇P + µ ∇ ~u + ∇(∇ · ~u) − ρ ∇Φ
dt 3
dε
Energy Eq. ρ = −P ∇ · ~u − ∇ · F~cond − L + V
dt
∂ui
Diss/Cond/Rad V ≡ τik , Fcond,k = hρεwk i , L≡C −H
∂xk
2
Diss/Cond/Rad stands for Dissipation, Conduction, Radiation
37
CHAPTER 6
Equations of State
Closure: The hydrodynamic equations for an ideal fluid (continuity eq, momentum
eqs, and energy eq) are 5 equations with 7 unknowns: (ρ, ~u, P , T (or ε), and Φ.
With the addition of the Poisson equation, which relates ρ and Φ. The seventh
and final equation that ensures closure is the equation of state P = P (ρ, T ). Note
that if the EoS is barotropic, i.e., P = P (ρ), then the continuity, momentum and
Poisson equations for a closed set, and the energy equation is not required.
Ideal Gas: a hypothetical gas that consists of identical point particles (i.e. of zero
volume) that undergo perfectly elastic collisions and for which interparticle forces
can be neglected.
An ideal gas obeys the ideal gas law: P V = N kB T .
kB T
P = P (ρ, T ) = ρ
µ mp
NOTE: astrophysical gases are often well described by the ideal gas law. Even for a
fully ionized gas, the interparticle forces (Coulomb force) can typically be neglected
(i.e., the potential energies involved are typically < 10% of the kinetic energies).
Ideal gas law breaks down for dense, and cool gases, such as those present in gaseous
planets.
38
Maxwell-Boltzmann Distribution: the distribution of particle momenta, p~ =
m~v , of an ideal gas follows the Maxwell-Boltzmann distribution.
3/2
3 1 p2
P(~p) d p~ = exp − d3 p~
2πmkB T 2mkB T
where p2 = ~p · ~p. This distribution follows from maximizing entropy under the
following assumptions:
NOTE: if there are temperature gradients in the gas, then the particle momenta only
follow the Maxwell-Boltzmann distribution locally, with T = T (~x) being the local
temperature.
P = ζ n hEi
where ζ = 2/3 (ζ = 1/3) in the case of a non-relativistic (relativistic) fluid, and
Z ∞
hEi = E P(E) dE
0
is the average, translational energy of the particles. In the case of our ideal (non-
relativistic) fluid,
2 Z ∞ 2
p p 3
hEi = = P(p) dp = kB T
2m 0 2m 2
39
Hence, we find that the EoS for an ideal gas is indeed given by
2 kB T
P = n hEi = n kB T = ρ
3 µmp
Specific Internal Energy: the internal energy per unit mass for an ideal gas is
hEi 3 kB T
ε= =
µmp 2 µmp
Actually, the above derivation is only valid for a true ‘ideal gas’, in which the particles
are point particles. More generally,
1 kB T
ε=
γ − 1 µmp
where γ is the adiabatic index, which for an ideal gas is equal to γ = (q +5)/(q +3),
with q the internal degrees of freedom of the fluid particles: q = 0 for point particles
(resulting in γ = 5/3), while diatomic particles have q = 2 (at sufficiently low
temperatures, such that they only have rotational, and no vibrational degrees of
freedom). The fact that q = 2 in that case arises from the fact that a diatomic
molecule only has two relevant rotation axes; the third axis is the symmetry axis of
the molecule, along which the molecule has negligible (zero in case of point particles)
moment of inertia. Consequently, rotation around this symmetry axis carries no
energy.
Photon gas: Having discussed the EoS of an ideal gas, we now focus on a gas of
photons. Photons have energy E = hν and momentum p = E/c = hν/c, with h the
Planck constant.
Black Body: an idealized physical body that absorbs all incident radiation. A black
body (BB) in thermal equilibrium emits electro-magnetic radiation called black
body radiation.
The spectral number density distribution of BB photons is given by
8πν 2 1
nγ (ν, T ) = 3 hν/k T −1
c e B
40
which implies a spectral energy distribution
8πhν 3 1
u(ν, T ) = nγ (ν, T ) hν =
c3 ehν/kB T − 1
and thus an energy density of
Z ∞
4σSB 4
u(T ) = u(ν, T ) dν = T ≡ ar T 4
0 c
where
2π 5 kB4
σSB =
15h3 c2
is the Stefan-Boltzmann constant and ar ≃ 7.6 × 10−15 erg cm−3 K−4 is called the
radiation constant.
Radiation Pressure: when the photons are reflected off a wall, or when they
are absorbed and subsequently re-emitted by that wall, they transfer twice their
momentum in the normal direction to that wall. Since photons are relativistic, we
have that the EoS for a photon gas is given by
1 1 1 aT 4
P = n hEi = nγ hhνi = u(T ) =
3 3 3 3
where we have used that u(T ) = nγ hEi.
41
for quarks). Finally, µ is called the chemical potential, and is a form of potential
energy that is related (in a complicated way) to the number density and temperature
of the particles (see Appendix G).
Classical limit: In the limit where the mean interparticle separation is much larger
than the de Broglie wavelength of the particles, so that quantum effects (e.g., Heisen-
berg’s uncertainty principle) can be ignored, the above distribution function of mo-
menta can be accurately approximated by the Maxwell-Boltzmann distribution.
Pauli Exclusion Principle: no more than one fermion of a given spin state can
occupy a given phase-space element h3 . Hence, for electrons, which have g = 2, the
maximum phase-space density is 2/h3 .
42
fully degenerate, then
N 3
Vx Vp = h
2
Using that ne = N/Vx , we find that
1/3
3
pF = ne h
8π
43
White Dwarfs and the Chandrasekhar limit: White dwarfs are the end-states
of stars with mass low enough that they don’t form a neutron star. When the
pressure support from nuclear fusion in a star comes to a halt, the core will start
to contract until degeneracy pressure kicks in. The star consists of a fully ionized
plasma. Assume for simplicity that the plasma consists purely of hydrogen, so that
the number density of protons is equal to that of electrons: np = ne . Because of
equipartition
p2p p2
= e
2mp 2me
p
Since mp ≫ me we have also that pp ≫ pe (in fact pp /pe = mp /me ≃ 43).
Consequently, when cooling or compressing the core of a star, the electrons will
become degenerate well before the protons do. Hence, white dwarfs are held up
against collapse by the degeneracy pressure from electrons. Since the electrons
are typically non-relativistic, the EoS of the white dwarf is: P ∝ ρ5/3 . If the white
dwarf becomes more and more massive (i.e., because it is accreting mass from a
companion star), the Pauli-exclusion principle causes the Fermi momentum, pF , to
increase to relativistic values. This softens the EoS towards P ∝ ρ4/3 . Such an
equation of state is too soft to stabilize the white dwarf against gravitational collapse;
the white dwarf collapses until it becomes a neutron star, at which stage it is
supported against further collapse by the degeneracy pressure from neutrons. This
happens when the mass of the white dwarf reaches Mlim ≃ 1.44M⊙ , the so-called
Chandrasekhar limit.
Non-Relativistic Relativistic
non-degenerate P ∝ ρT P ∝ T4
degenerate P ∝ ρ5/3 P ∝ ρ4/3
Summary of equations of state for different kind of fluids
44
CHAPTER 7
Vorticity: The vorticity of a flow is defined as the curl of the velocity field:
vorticity : ~ = ∇ × ~u
w
It is a microscopic measure of rotation (vector) at a given point in the fluid, which
can be envisioned by placing a paddle wheel into the flow. If it spins about its axis
at a rate Ω, then w = |w|
~ = 2Ω.
Circulation: The circulation around a closed contour C is defined as the line integral
of the velocity along that contour:
I Z
circulation : ΓC = ~u · d~l = w ~
~ · dS
C S
Vortex line: a line that points in the direction of the vorticity vector. Hence, a
vortex line relates to w,
~ as a streamline relates to ~u (cf. Chapter 2).
In an inviscid fluid the vortex lines/tubes move with the fluid: a vortex line an-
chored to some fluid element remains anchored to that fluid element.
45
Figure 4: Evolution of a vortex tube. Solid dots correspond to fluid elements. Due
to the shear in the velocity field, the vortex tube is stretched and tilted. However, as
long as the fluid is inviscid and barotropic Kelvin’s circularity theorem assures that
the circularity is conserved with time. In addition, since vorticity is divergence-free
(‘solenoidal’), the circularity along different cross sections of the same vortex-tube is
the same.
46
∂w
~ ∇P
= ∇ × (~u × w)
~ −∇× + ν∇2 w
~
∂t ρ
~ = ∇S × A
To write this in Lagrangian form, we first use that ∇ × (S A) ~ + S (∇ × A)
~
[see Appendix A] to write
1 1 1 ρ∇(1) − 1∇ρ ∇P × ∇ρ
∇ × ( ∇P ) = ∇( ) × ∇P + (∇ × ∇P ) = × ∇P =
ρ ρ ρ ρ2 ρ2
where we have used, once more, that curl(grad S) = 0. Next, using the vector
identities from Appendix A, we write
∇ × (w
~ × ~u) = w(∇
~ · ~u) − (w
~ · ∇)~u − ~u(∇ · w)
~ + (~u · ∇)w
~
The third term vanishes because ∇ · w
~ = ∇ · (∇ × ~u) = 0. Hence, using that ∂ w/∂t
~ +
(~u · ∇)w
~ = dw/dt
~ we finally can write the vorticity equation in Lagrangian
form:
dw
~ ∇ρ × ∇P
~ · ∇)~u − w(∇
= (w ~ · ~u) + 2
+ ν∇2 w
~
dt ρ
This equation describes how the vorticity of a fluid element evolves with time. We
now describe the various terms of the rhs of this equation in turn:
• (w
~ · ∇)~u: This term represents the stretching and tilting of vortex tubes due
to velocity gradients. To see this, we pick w
~ to be pointing in the z-direction.
Then
• w(∇
~ · ~u): This term describes stretching of vortex tubes due to flow com-
pressibility. This term is zero for an incompressible fluid or flow (∇ · ~u = 0).
Note that, again under the assumption that the vorticity is pointing in the
z-direction,
47
∂ux ∂uy ∂uz
w(∇
~ · ~u) = wz + + ~ez
∂x ∂y ∂z
• ν∇2 w:
~ This term describes the diffusion of vorticity due to viscosity, and
is obviously zero for an inviscid fluid (ν = 0). Typically, viscosity gener-
ates/creates vorticity at a bounding surface: due to the no-slip boundary con-
dition shear arises giving rise to vorticity, which is subsequently diffused into
the fluid by the viscosity. In the interior of a fluid, no new vorticity is generated;
rather, viscosity diffuses and dissipates vorticity.
• ∇ × F~ : There is a fifth term that can create vorticity, which however does not
appear in the vorticity equation above. The reason is that we assumed that the
only external force is gravity, which is a conservative force and can therefore be
written as the gradient of a (gravitational) potential. More generally, though,
there may be non-conservative, external body forces present, which would give
rise to a ∇ × F~ term in the rhs of the vorticity equation. An example of a non-
conservative force creating vorticity is the Coriolis force, which is responsible
for creating hurricanes.
48
Figure 5: The baroclinic creation of vorticity in a pyroclastic flow. High density fluid
flows down a mountain and shoves itself under lower-density material, thus creating
non-zero baroclinicity.
Using the definition of circulation, it can be shown (here without proof) that
Z
dΓ ∂w~ ~
= + ∇ × (w~ × ~u) · dS
dt S ∂t
Using the vorticity equation, this can be rewritten as
Z
dΓ ∇ρ × ∇P
= 2
+ ν∇ w~ + ∇ × F~ · dS
~
dt S ρ2
49
NOTE: By comparing the equations expressing dw/dt ~ and dΓ/dt it is clear that
the stretching a tilting terms present in the equation describing dw/dt,
~ are absent
in the equation describing dΓ/dt. This implies that stretching and tilting changes
the vorticity, but keeps the circularity invariant. This is basically the first theorem
of Helmholtz described below.
Kelvin’s Circulation Theorem: The number of vortex lines that thread any
element of area that moves with the fluid (i.e., the circulation) remains unchanged
in time for an inviscid, barotropic fluid, in the absence of non-conservative forces.
We end this chapter on vorticity and circulation with the three theorems of Helmholtz,
which hold in the absence of non-conservative forces (i.e., F~ = 0).
where A1 and A2 are the areas of the cross sections that bound the volume V of the
vortex tube. Using Stokes’ curl theorem, we have that
50
Z I
~ · n̂ dA =
w ~u · d~l
A C
Hence we have that ΓC1 = ΓC2 where C1 and C2 are the curves bounding A1 and A2 ,
respectively.
Helmholtz Theorem 2: A vortex line cannot end in a fluid. Vortex lines and tubes
must appear as closed loops, extend to infinity, or start/end at solid boundaries.
51
Figure 6: A beluga whale demonstrating Kelvin’s circulation theorem and Helmholtz’
second theorem by producing a closed vortex tube under water, made out of air.
52
CHAPTER 8
Having derived all the relevant equations for hydrodynamics, we now start examining
several specific flows. Since a fully general solution of the Navier-Stokes equation is
(still) lacking (this is one of the seven Millenium Prize Problems, a solution of which
will earn you $1,000,000), we can only make progress if we make several assumptions.
We start with arguably the simplest possible flow, namely ‘no flow’. This is the area
of hydrostatics in which ~u(~x, t) = 0. And since we seek a static solution, we also
must have that all ∂/∂t-terms vanish. Finally, in what follows we shall also ignore
radiative processes (i.e., we set L = 0).
Applying these restrictions to the continuity, momentum and energy equations (see
box at the end of Chapter 5) yields the following two non-trivial equations:
∇P = −ρ ∇Φ
∇ · F~cond = 0
The first equation is the well known equation of hydrostatic equilibrium, stating
that the gravitational force is balanced by pressure gradients, while the second equa-
tion states that in a static fluid the conductive flux needs to be divergence-free.
To further simplify matters, let’s assume (i) spherical symmetry, and (ii) a barotropic
equation of state, i.e., P = P (ρ).
dP G M(r) ρ(r)
=−
dr r2
53
In addition, if the gas is self-gravitating (such as in a star) then we also have that
dM
= 4πρ(r) r 2
dr
For a barotropic EoS this is a closed set of equations, and the density profile can be
solved for (given proper boundary conditions). Of particular interest in astrophysics,
is the case of a polytropic EoS: P ∝ ρΓ , where Γ is the polytropic index. Note
that Γ = 1 and Γ = γ for isothermal and adiabatic equations of state, respectively.
A spherically symmetric, polytropic fluid in HE is called a polytropic sphere.
Here n = 1/(Γ − 1) is related to the polytropic index (in fact, confusingly, some texts
refer to n as the polytropic index),
1/2
4πGρc
ξ= r
Φ0 − Φc
is a dimensionless radius,
Φ0 − Φ(r)
θ=
Φ0 − Φc
with Φc and Φ0 the values of the gravitational potential at the center (r = 0) and
at the surface of the star (where ρ = 0), respectively. The density is related to θ
according to ρ = ρc θn with ρc the central density.
54
(see Chapter 6) and is therefore described by a polytrope of index n = 3/2. In the
relativistic case P ∝ ρ4/3 which results in a polytrope of index n = 3.
Heat transport in stars: Typically, ignoring abundance gradients, stars have the
equation of state of an ideal gas, P = P (ρ, T ). This implies that the equations of
stellar structure need to be complemented by an equation of the form
dT
= F (r)
dr
Since T is a measure of the internal energy, the rhs of this equation describes the
heat flux, F (r).
55
Recall from Chapter 4 that the thermal conductivity K ∝ (kB T )1/2 /σ where σ
is the collisional cross section. Using that kB T ∝ v 2 and that the mean-free path of
the particles is λmfp = 1/(nσ), we have that
K ∝ n λmfp v
with v the thermal, microscopic velocity of the particles (recall that ~u = 0). Since
radiative heat transport in a star is basically the conduction of photons, and since
c ≫ ve and the mean-free part of photons is much larger than that of electrons (after
all, the cross section for Thomson scattering, σT , is much smaller than the typical
cross section for Coulomb interactions), we have that in stars radiation is a far more
efficient heat transport mechanism than conduction. An exception are relativistic,
degenerate cores, for which ve ∼ c and photons and electrons have comparable mean-
free paths.
Trivia: On average it takes ∼ 200.000 years for a photon created at the core of the
Sun in nuclear burning to make its way to the Sun’s photosphere; from there it only
takes ∼ 8 minutes to travel to the Earth.
Hydrostatic Mass Estimates: Now let us consider the case of an ideal gas, for
which
kB T
P = ρ,
µmp
but this time the gas is not self-gravitating; rather, the gravitational potential may
be considered ‘external’. A good example is the ICM; the hot gas that permeates
clusters. From the EoS we have that
56
dP ∂P dρ ∂P dT P dρ P dT
= + = +
dr ∂ρ dr ∂T dr ρ dr T dr
P r dρ r dT P d ln ρ d ln T
= + = +
r ρ dr T dr r d ln r d ln r
Substitution of this equation in the equation for Hydrostatic equilibrium (HE) yields
kB T (r) r d ln ρ d ln T
M(r) = − +
µmp G d ln r d ln r
This equation is often used to measure the ‘hydrostatic’ mass of a galaxy cluster;
X-ray measurements can be used to infer ρ(r) and T (r) (after deprojection, which is
analytical in the case of spherical symmetry). Substitution of these two radial depen-
dencies in the above equation then yields an estimate for the cluster’s mass profile,
M(r). Note, though, that this mass estimate is based on three crucial assump-
tions: (i) sphericity, (ii) hydrostatic equilibrium, and (iii) an ideal-gas EoS. Clusters
typically are not spherical, often are turbulent (such that ~u 6= 0, violating the as-
sumption of HE), and can have significant contributions from non-thermal pressure
due to magnetic fields, cosmic rays and/or turbulence. Including these non-thermal
pressure sources the above equation becomes
kB T (r) r d ln ρ d ln T Pnt d ln Pnt
M(r) = − + +
µmp G d ln r d ln r Pth d ln r
were Pnt and Pth are the non-thermal and thermal contributions to the total gas
pressure. Unfortunately, it is extremely difficult to measure Pnt reliably, which is
therefore often ignored. This may result in systematic biases of the inferred cluster
mass (typically called the ‘hydrostatic mass’).
The Solar corona is a large, spherical region of hot (T ∼ 106 K) plasma extending
well beyond its photosphere. Let’s assume that the heat is somehow (magnetic
reconnection?) produced in the lower layers of the corona, and try to infer the density,
temperature and pressure profiles under the assumption of hydrostatic equilibrium.
57
We have the boundary condition of the temperature at the base, which we assume
to be T0 = 3 × 106 K, at a radius of r = r0 ∼ R⊙ ≃ 6.96 × 1010 cm. The mass of the
corona is negligble, and we therefore have that
dP G M⊙ µmp P
= − 2
dr r kB T
d dT
K r2 = 0
dr dr
where we have used the ideal gas EoS to substitute for ρ. Note that the latter of these
equations follows from ∇ · F~ = 0, which is the energy equation in HE. As we have
seen above K ∝ nλmfp T 1/2 . In a plasma one furthermore has that λmfp ∝ n−1 T 2 ,
which implies that K ∝ T 5/2 . Hence, the second equation can be written as
dT
r 2 T 5/2 = constant
dr
which implies
−2/7
r
T = T0
r0
Note that this equation satisfies our boundary condition, and that T∞ = limr→∞ T (r) =
0. Substituting this expression for T in the HE equation yields
dP G M⊙ µmp dr
=− 2/7 r 12/7
P kB T0 r0
Note that
7 G M⊙ µmp
lim P = P0 exp − 6 0
=
r→∞ 5 kB T0 r0
58
Hence, you need an external pressure to confine the corona. Well, that seems OK,
given that the Sun is embedded in an ISM, whose pressure we can compute taking
characteristic values for the warm phase (T ∼ 104 K and n ∼ 1 cm−3 ). Note that the
other phases (cold and hot) have the same pressure. Plugging in the numbers, we
find that
P∞ ρ0
∼ 10
PISM ρISM
Since ρ0 ≫ ρISM we thus infer that the ISM pressure falls short, by orders of magni-
tude, to be able to confine the corona....
As first inferred by Parker in 1958, the correct implication of this puzzling result is
that a hydrostatic corona is impossible; instead, Parker made the daring suggestion
that there should be a solar wind, which was observationally confirmed a few years
later.
————————————————-
Having addressed hydrostatics (‘no flow’), we now consider the next simplest flow;
steady flow, which is characterised by ~u(~x, t) = ~u(~x). For steady flow ∂~u/∂t = 0,
and fluid elements move along the streamlines (see Chapter 2).
The enthalpy, H, is a measure for the total energy of a thermodynamic system that
includes the internal energy, U, and the amount of energy required to make room
for it by displacing its environment and establishing its volume and pressure:
59
H = U + PV
The differential of the enthalpy can be written as
dH = dU + P dV + V dP
Using the first law of thermodynamics, according to which dU = dQ − P dV , and
the second law of thermodynamics, according to which dQ = T dS, we can rewrite
this as
dH = T dS + V dP
which, in specific form, becomes
dP
dh = T ds +
ρ
(i.e., we have s = S/m). This relation is one of the Gibbs relations frequently
encountered in thermodynamics. NOTE: for completeness, we point out that this
expression ignores changes in the chemical potential (see Appendix J).
∇P
= ∇h − T ∇s
ρ
(for a formal proof, see at the end of this chapter). Now recall from the previous
chapter on vorticity that the baroclinic term is given by
∇P ∇ρ × ∇P
∇× =
ρ ρ2
Using the above relation, and using that the curl of the gradient of a scalar vanishes,
we can rewrite this baroclinic term as ∇ × (T ∇s). Now using that ∇ × S A ~ =
∇S × A ~ + S(∇ × A)
~ for a scalar S and a vector A ~ (see Appendix A), we can write
∇×(T ∇s) = ∇T ×∇s. This is another form for the baroclinic term. It demonstrates
that another way to create baroclinicity, and thus vorticity, is by having the gradient
in entropy be misaligned with the gradient in temperature.
60
Using the momentum equation for a steady, ideal fluid, and substituting ∇P/ρ →
∇h − T ∇s, we obtain
∇B = T ∇s + ~u × w
~
u2 u2
B≡ +Φ+h= + Φ + ε + P/ρ
2 2
Let’s investigate what happens to the Bernoulli function for an ideal fluid in a
steady flow. We start by pointing out that in an ideal fluid there is no conduction
and no dissipation. As a consequence, the flow of an ideal fluid conserves entropy;
ds/dt = 0. We say that the flow is isentropic.
Using that
ds ∂s
= + u · ∇s = 0
dt ∂t
61
we see that for a steady flow of ideal fluid we always have that u · ∇s = 0. In words,
there can’t be an entropy gradient in the direction of the flow.
~u · ∇s = 0 and ~u · ∇B = 0
Furthermore, using that
dB ∂B
= + ~u · ∇B
dt ∂t
we immediately see that steady flow of ideal fluid obeys
dB
=0
dt
62
Hence, steady flow of an ideal fluid conserves both entropy and the Bernoulli function!
Using the definition of the Bernoulli function we can write this as
dB d~u dΦ ds 1 dP
= ~u · + +T + =0
dt dt dt dt ρ dt
Since ds/dt = 0 for an ideal fluid, we have that if the flow is such that the gravita-
tional potential along the flow doesn’t change significantly (such that dΦ/dt ≃ 0),
we find that
d~u 1 dP
~u · =−
dt ρ dt
This is known as Bernoulli’s theorem, and states that as the speed of a steady flow
increases, the internal pressure of the ideal fluid must decrease (this can be rather
counter-intuitive). Applications of Bernoulli’s theorem discussed in class include the
shower curtain and the pitot tube (a flow measurement device used to measure fluid
flow velocity).
∇B − T ∇s = ~u × w
~
Suppose we have an isentropic, irrotational fluid, which means that ∇s = 0 and
~ = 0 everywhere. Then we also have that ∇B = 0; the Bernoulli function is
w
everywhere the same. Now suppose such a flow encounters a shock. As we will
see in Chapter 12, if the shock is adiabatic in that no radiative losses occur, then
conservation of energy implies that the Bernoulli function remains constant across
the shock. However, a shock will increase the entropy of the gas. And if the strength
of the shock varies in the direction perpendicular to the flow (for example, when the
flow hits a curved shock), then behind the shock we have ∇B = 0 but ∇s 6= 0.
This implies that behind the shock we must have created vorticity. This is one of the
most important mechanisms for creating vorticity in astrophysics!
————————————————-
63
Potential flow: The final flow to consider in this chapter is potential flow. Consider
~ ≡ ∇ × ~u = 0 everywhere. This implies that
an irrotational flow, which satisfies w
there is a scalar function, φu (x), such that ~u = ∇φu , which is why φu (x) is called
the velocity potential. The corresponding flow ~u(~x) is called potential flow.
If the fluid is ideal (i.e., ν = K = 0), and barotropic or isentropic, such that the flow
fluid has vanishing baroclinicity, then Kelvin’s circulation theorem assures that
the flow will remain irrotational throughout (no vorticity can be created), provided
that all forces acting on the fluid are conservative.
∇ · ~u = ∇2 φu = 0
This is the well known Laplace equation, familiar from electrostatics. Mathemat-
ically, this equation is of the elliptic PDE type which requires well defined boundary
conditions in order for a solution to both exist and be unique. A classical case of
potential flow is the flow around a solid body placed in a large fluid volume. In this
case, an obvious boundary condition is the one stating that the velocity component
perpendicular to the surface of the body at the body (assumed at rest) is zero. This
is called a Neumann boundary condition and is given by
∂φu
= ~n · ∇φu = 0
∂n
with ~n the normal vector. The Laplace equation with this type of boundary condition
constitutes a well-posed problem with a unique solution. An example of potential
flow around a solid body is shown in Fig. 2 in Chapter 2. We will not examine any
specific examples of potential flow, as this means having to solve a Laplace equation,
which is purely a mathematical exersize. We end, though, by pointing out that real
fluids are never perfectly inviscid (ideal fluids don’t exist). And any flow past a
surface involves a boundary layer inside of which viscosity creates vorticity (due to
no-slip boundary condition, which states that the tangential velocity at the surface
of the body must vanish). Hence, potential flow can never fully describe the flow
around a solid body; otherwise one would run into d’Alembert’s paradox which
is that steady potential flow around a body exerts zero force on the body; in other
words, it costs no energy to move a body through the fluid at constant speed. We
know from everyday experience that this is indeed not true. The solution to the
64
paradox is that viscosity created in the boundary layer, and subsequently dissipated,
results in friction.
Although potential flow around an object can thus never be a full description of the
flow, in many cases, the boundary layer is very thin, and away from the boundary
layer the solutions of potential flow still provide an accurate description of the flow.
————————————————-
To see this, use that the natural variables of h are the specific entropy, s, and the
pressure P . Hence, h = h(s, P ), and we thus have that
∂h ∂h
dh = ds + dP
∂s ∂P
From a comparison with the previous expression for dh, we see that
∂h ∂h 1
=T, =
∂s ∂P ρ
which allows us to derive
∂h ∂h ∂h
∇h = ~ex + ~ey + ~ez
∂x ∂y ∂z
∂h ∂s ∂h ∂P ∂h ∂s ∂h ∂P ∂h ∂s ∂h ∂P
= + ~ex + + ~ey + + ~ez
∂s ∂x ∂P ∂x ∂s ∂y ∂P ∂y ∂s ∂z ∂P ∂z
∂h ∂s ∂s ∂s ∂h ∂P ∂P ∂P
= ~ex + ~ey + ~ez + ~ex + ~ey + ~ez
∂s ∂x ∂y ∂z ∂P ∂x ∂y ∂z
1
= T ∇s + ∇P
ρ
————————————————-
65
CHAPTER 9
As we have seen in our discussion on potential flow in the previous chapter, realistic
flow past an object always involves a boundary layer in which viscosity results in
vorticity. Even if the viscosity of the fluid is small, the no-slip boundary condition
typically implies a region where the shear is substantial, and viscocity thus manifests
itself.
In this chapter we examine two examples of viscous flow. We start with a well-
known example from engineering, known as Poiseuille-Hagen flow through a pipe.
Although not really an example of astrophysical flow, it is a good illustration of how
viscosity manifests itself as a consequence of the no-slip boundary condition. The
second example that we consider is viscous flow in a thin accretion disk. This flow,
which was first worked out in detail in a famous paper by Shakura & Sunyaev in
1973, is still used today to describe accretion disks in AGN and around stars.
————————————————-
Pipe Flow: Consider the steady flow of an incompressible viscous fluid through
a circular pipe of radius Rpipe and lenght L. Let ρ be the density of the fluid as
it flows through the pipe, and let ν = µ/ρ be its kinetic viscosity. Since the
flow is incompressible, we have that fluid density will be ρ throughout. If we pick a
Cartesian coordinate system with the z-axis along the symmetry axis of the cylinder,
then the velocity field of our flow is given by
~u = uz (x, y, z) ~ez
In other words, ux = uy = 0.
66
Figure 7: Poiseuille-Hagen flow of a viscous fluid through a pipe of radius Rpipe and
lenght L.
and using that all partial time-derivatives of a steady flow vanish, we obtain that
∂ρux ∂ρuy ∂ρuz ∂uz
+ + =0 ⇒ =0
∂x ∂y ∂z ∂z
where we have used that ∂ρ/∂z = 0 because of the incompressibility of the flow.
Hence, we can update our velocity field to be ~u = uz (x, y) ~ez .
Next we write down the momentum equations for a steady, incompressible flow,
which are given by
∇P
(~u · ∇)~u = − + ν∇2 ~u − ∇Φ
ρ
In what follows we assume the pipe to be perpendicular to ∇Φ, so that we may
ignore the last term in the above expression. For the x- and y- components of the
momentum equation, one obtains that ∂P/∂x = ∂P/∂y = 0. For the z-component,
we instead have
∂uz 1 ∂P
uz =− + ν∇2 uz
∂z ρ ∂z
Combining this with our result from the continuity equation, we obtain that
1 ∂P
= ν∇2 uz
ρ ∂z
Next we use that ∂P/∂z cannot depend on z; otherwise uz would depend on z, but
according to the continuity equation ∂uz /∂z = 0. This means that the pressure
67
gradient in the z-direction must be constant, which we write as −∆P/L, where ∆P
is the pressure difference between the beginning and end of the pipe, and the minus
sign us used to indicate that the fluid pressure declines as it flows throught the pipe.
∆P 2
uz (R) = Rpipe − R2
4ρ ν L
As is evident from the above expression, for a given pressure difference ∆P , the flow
speed u ∝ ν −1 (i.e., a more viscous fluid will flow slower). In addition, for a given
fluid viscosity, applying a larger pressure difference ∆P results in a larger flow speed
(u ∝ ∆P ).
Now let us compute the amount of fluid that flows through the pipe per unit time:
R
Zpipe
π ∆P 4
Ṁ = 2π ρ uz (R) R dR = R
8 ν L pipe
0
Note the strong dependence on the pipe radius; this makes it clear that a clogging of
the pipe has a drastic impact on the mass flow rate (relevant for both arteries and oil-
pipelines). The above expression also gives one a relatively easy method to measure
68
the viscosity of a fluid: take a pipe of known Rpipe and L, apply a pressure difference
∆P across the pipe, and measure the mass flow rate, Ṁ; the above expression allows
one to then compute ν.
The Poiseuille velocity flow field has been experimentally confirmed, but only for
slow flow! When |~u| gets too large (i.e., ∆P is too large), then the flows becomes
irregular in time and space; turbulence develops and |~u| drops due to the enhanced
drag from the turbulence. This will be discussed in more detail in Chapter 12.
————————————————-
Accretion Disks: We now move to a viscous flow that is more relevant for as-
trophysics; accretion flow. Consider a thin accretion disk surrounding an accreting
object of mass M• ≫ Mdisk (such that we may ignore the disk’s self-gravity). Because
of the symmetries involved, we adopt cylindrical coordinates, (R, θ, z), with the
z-axis perpendicular to the disk. We also have that ∂/∂θ is zero, and we set uz = 0
throughout.
Let’s start with the continuity equation, which in our case reads
∂ρ 1 ∂
+ (R ρ uR ) = 0
∂t R ∂R
(see Appendix D for how to express the divergence in cylindrical coordinates).
NOTE: There are several terms in the above expression that may seem ‘surprising’.
The important thing to remember in writing down the equations in curvi-linear
69
coordinates is that operators can also act on unit-direction vectors. For example,
the θ-component of ∇2~u is NOT ∇2 uθ . That is because the operator ∇2 acts on
uR~eR + uθ~eθ + uz~ez , and the directions of ~eR and ~eθ depend on position! The same
holds for the convective operator (~u · ∇) ~u. The full expressions for both cylindrical
and spherical coordinates are written out in Appendix D.
Setting all the terms containing ∂/∂θ and/or uz to zero, the Navier-Stokes equation
simplifies considerably to
2
∂uθ ∂uθ uR uθ ∂ uθ ∂ 2 uθ 1 ∂uθ uθ
ρ + uR + =µ + + − 2
∂t ∂R R ∂R2 ∂z 2 R ∂R R
where we have replaced the kinetic viscosity, ν, with µ = νρ.
Next we multiply the continuity equation by Ruθ which we can then write as
∂(Σ R uθ ) ∂(Ruθ ) ∂(Σ R uR uθ ) ∂uθ
−Σ + − R Σ uR =0
∂t ∂t ∂R ∂R
Adding this to R times the Navier-Stokes equation, and rearranging terms, yields
∂(Σ R uθ ) ∂(Σ R uR uθ )
+ + Σ uR uθ = G(µ, R)
∂t ∂R
where G(µ, R) = RF (µ). Next we introduce the angular frequency Ω ≡ uθ /R
which allows us to rewrite the above expression as
∂(Σ R2 Ω) 1 ∂
+ Σ R3 Ω uR = G(µ, R)
∂t R ∂R
70
Note that Σ R2 Ω = Σ R uθ is the angular momentum per unit surface area. Hence
the above equation describes the evolution of angular momentum in the accretion
disk. It is also clear, therefore, that G(µ, R) must describe the viscous torque on
the disk material, per unit surface area. To derive an expression for it, recall that
Z 2
∂ uθ 1 ∂uθ uθ
G(µ, R) = R dz µ + − 2
∂R2 R ∂R R
where we have ignored the ∂ 2 uθ /∂z 2 term which is assumed to be small. Using that
µ = νρ and that µ is independent of R and z (this is an assumption that underlies
the Navier-Stokes equation from which we started) we have that
2
∂ uθ 1 ∂uθ uθ
G(µ, R) = ν R Σ + − 2
∂R2 R ∂R R
Next we use that uθ = Ω R to write
∂uθ dΩ
= Ω+R
∂R dR
Substituting this in the above expression for G(µ, R) yield
2
2d Ω dΩ 1 ∂ 3 dΩ
G(µ, R) = ν Σ R + 3R = ν ΣR
dR2 dR R ∂R dR
Substituting this expression for the viscous torque in the evolution equation for the
angular momentum per unit surface density, we finally obtain the full set of equations
that govern our thin accretion disk:
∂ 1 ∂ 1 ∂ 3 dΩ
Σ R2 Ω + Σ R 3 Ω uR = ν ΣR
∂t R ∂R R ∂R dR
∂Σ 1 ∂
+ (R Σ uR ) = 0
∂t R ∂R
1/2
G M•
Ω=
R3
71
These three equations describe the dynamics of a thin, viscous accretion disk. The
third equation indicates that we assume that the fluid is in Keplerian motion around
the accreting object of mass M• . As discussed further below, this is a reasonable
assumption as long as the accretion disk is thin.
Ṁ (R) = −2πΣ R uR
Now let us consider a steady accretion disk. This implies that ∂/∂t = 0 and
that Ṁ (R) = Ṁ ≡ Ṁ• (the mass flux is constant throughout the disk, otherwise
∂Σ/∂t 6= 0). In particular, the continuity equation implies that
R Σ u R = C1
Using the above expression for the mass inflow rate, we see that
Ṁ•
C1 = −
2π
Ṁ•
C2 = R•2 Ω• C1 = − (G M• R• )1/2
2π
72
we have that
−1
Ṁ• 2 3 dΩ
νΣ = − R Ω + (G M• R• )1/2 R
2π dR
" 1/2 #
Ṁ• R•
= + 1−
3π R
This shows that the mass inflow rate and kinetic viscosity depend linearly on each
other.
The gravitational energy lost by the inspiraling material is converted into heat. This
is done through viscous dissipation: viscosity robs the disk material of angular
momentum which in turn causes it to spiral in.
(see Chapter 4). Note that the last term in the above expression vanishes because
the fluid is incompressible, such that
" 2 #
∂ui ∂uj ∂ui
V=µ +
∂xj ∂xi ∂xj
In our case, using that ∂/∂θ = ∂/∂z = 0 and that uz = 0, the only surviving terms
are
" 2 2 # " 2 2 #
∂uR ∂uθ ∂uR ∂uR ∂uR ∂uθ
V =µ + + =µ 2 +
∂R ∂R ∂R ∂R ∂R ∂R
73
If we make the reasonable assumption that uR ≪ uθ , we can ignore the first term,
such that we finally obtain
2 2
∂uθ 2 dΩ
V=µ = µR
∂R dR
which expresses the viscous dissipation per unit volume. Note that there is no viscous
dissipation if dΩ/dR = 0, i.e., in the case of solid body rotation. This makes sense
since in that case there is no shear in the disk.
Using once more that dΩ/dR = −(3/2)Ω/R, and integrating over the entire disk
yields the accretion luminosity of a thin accretion disk:
Z∞
dE G M• Ṁ•
Lacc ≡ 2π R dR =
dt 2 R•
R•
To put this in perspective, realize that the gravitation energy of mass m at radius
R• is G M• m/ R• . Thus, Lacc is exactly half of the gravitational energy lost due to
the inflow. This obviously begs the question where the other half went...The answer
is simple; it is stored in kinetic energy at the ‘boundary’ radius R• of the accreting
flow.
We end our discussion on accretion disks with a few words of caution. First of
all, our entire derivation is only valid for a thin accretion disk. In a thin disk, the
74
pressure in the disk must be small (otherwise it would puff up). This means that the
∂P/∂R term in the R-component of the Navier-Stokes equation is small compared
to ∂Φ/∂R = GM/R2 . This in turn implies that the gas will indeed be moving on
Keplerian orbits, as we have assumed. If the accretion disk is thick, the situation is
much more complicated, something that will not be covered in this course.
Finally, let us consider the time scale for accretion. As we have seen above, the
energy loss rate per unit surface area is
2
2 dΩ 9 G M•
ν ΣR = ν
dR 4 R3
We can compare this with the gravitational potential energy of disk material per
unit surface area, which is
G M• Σ
E=
R
E 4 R2 R2
tacc ≡ = ∼
dE/dt 9 ν ν
To estimate this time-scale, we first estimate the molecular viscosity. Recall that
ν ∝ λmfpv with v a typical velocity of the fluid particles. In virtually all cases
encountered in astrophysics, we have that the size of the accretion disk, R, is many,
many orders of magnitude larger than λmfp . As a consequence, the corresponding
tacc easily exceeds the Hubble time!
The conclusion is that molecular viscosity is way too small to result in any signif-
icant accretion in objects of astrophysical size. Hence, other source of viscosity are
required, which is a topic of ongoing discussion in the literature. Probably the most
promising candidates are turbulence (in different forms), and the magneto-rotational
instability (MRI). Given the uncertainties involved, it is common practive to simply
write −1
P 1 dΩ
ν=α
ρ R dR
where α is a ‘free parameter’. A thin accretion disk modelled this way is often called
an alpha-accretion disk. If you wonder what the origin is of the above expression;
75
Figure 8: Image of the central region of NGC 4261 taken with the Hubble Space
Telescope. It reveals a ∼ 100pc scale disk of dust and gas, which happens to be per-
pendicular to a radio jet that emerges from this galaxy. This is an alledged ‘accretion
disk’ supplying fuel to the central black hole in this galaxy.
it simply comes from assuming that the only non-vanishing off-diagonal term of the
stress tensor is taken to be αP (where P is the value along the diagonal of the stress
tensor).
76
CHAPTER 10
Turbulence
1
~u · ∇~u = ∇u2 − ~u × w
~
2
which describes the ”inertial acceleration” and is ultimately responsible for the origin
of the chaotic character of many flows and of turbulence. Because of this non-
linearity, we cannot say whether a solution to the Navier-Stokes equation with nice
and smooth initial conditions will remain nice and smooth for all time (at least not
in 3D).
Laminar flow: occurs when a fluid flows in parallel layers, without lateral mixing
(no cross currents perpendicular to the direction of flow). It is characterized by high
momentum diffusion and low momentum convection.
The Reynold’s number: In order to gauge the importance of viscosity for a fluid,
it is useful to compare
2 the ratio of the inertial acceleration (~u · ∇~u) to the viscous
1
acceleration (ν ∇ ~u + 3 ∇(∇ · ~u) ). This ratio is called the Reynold’s number, R,
and can be expressed in terms of the typical velocity scale U ∼ |~u| and length scale
L ∼ 1/∇ of the flow, as
~u · ∇~u U 2 /L UL
R= ∼ =
ν ∇2~u + 13 ∇(∇ · ~u) νU/L2 ν
If R ≫ 1 then viscosity can be ignored (and one can use the Euler equations to
describe the flow). However, if R ≪ 1 then viscosity is important.
77
Figure 9: Illustration of laminar vs. turbulent flow.
Similarity: Flows with the same Reynold’s number are similar. This is evident
from rewriting the Navier-Stokes equation in terms of the following dimensionless
variables
~u ~x U P Φ ˜ = L∇
ũ = x̃ = t̃ = t p̃ = Φ̃ = ∇
U L L ρ U2 U2
This yields (after multiplying the Navier-Stokes equation with L/U 2 ):
∂ ũ ˜ + ∇p̃
˜ +∇˜ Φ̃ = 1 1
˜ ũ + ∇(
2 ˜ ∇
˜ · ũ)
+ ũ · ∇ũ ∇
∂ t̃ R 3
which shows that the form of the solution depends only on R. This principle is
extremely powerful as it allows one to making scale models (i.e., when developing
airplanes, cars etc). NOTE: the above equation is only correct for an incompressible
fluid, i.e., a fluid that obeys ∇ρ = 0. If this is not the case the term P̃ (∇ρ/ρ) needs
to be added at the rhs of the equation, braking its scale-free nature.
78
Figure 10: Illustration of flows at different Reynolds number.
• R > 103 : vortices are unstable, resulting in a turbulent wake behind the
cylinder that is ‘unpredictable’.
79
Figure 11: The image shows the von Kármán Vortex street behind a 6.35 mm di-
ameter circular cylinder in water at Reynolds number of 168. The visualization was
done using hydrogen bubble technique. Credit: Sanjay Kumar & George Laughlin,
Department of Engineering, The University of Texas at Brownsville
The following movie shows a R = 250 flow past a cylinder. Initially one can witness
separation, and the creation of two counter-rotating vortices, which then suddenly
become ‘unstable’, resulting in the von Kármán vortex street:
[Link]
80
Figure 12: Typical Reynolds numbers for various biological organisms. Reynolds
numbers are estimated using the length scales indicated, the “rule-of-thumb” in the
text, and material properties of water.
81
Boundary Layers: Even when R ≫ 1, viscosity always remains important in thin
boundary layers adjacent to any solid surface. This boundary layer must exist in
order to satisfy the no-slip boundary condition. If the Reynolds number exceeds
a critical value, the boundary layer becomes turbulent. Turbulent layers and their
associated turbulent wakes exert a much bigger drag on moving bodies than their
laminar counterparts.
Momentum Diffusion & Reynolds stress: This gives rise to an interesting phe-
nomenon. Consider flow through a pipe. If you increase the viscosity (i.e., decrease
R), then it requires a larger force to achieve a certain flow rate (think of how much
harder it is to push honey through a pipe compared to water). However, this trend
is not monotonic. For sufficiently low viscosity (large R), one finds that the trend
reverses, and that is becomes harder again to push the fluid through the pipe. This
is a consequence of turbulence, which causes momentum diffusion within the flow,
which acts very much like viscosity. However, this momentum diffusion is not due
to the viscous stress tensor, τij , but rather to the Reynolds stress tensor Rij .
To understand the ‘origin’ of the Reynolds stress tensor,consider the following:
ui = ūi + u′i
This is knowns as the Reynolds decomposition. The ‘mean’ component can be a
time-average, a spatial average, or an ensemble average, depending on the detailed
characteristics of the flow. Note that this is reminiscent of how we decomposed the
microscopic velocities of the fluid particles in a ‘mean’ velocity (describing the fluid
elements) and a ‘random, microscopic’ velocity (~v = ~u + w).
~
Substituting this into the Navier-Stokes equation, and taking the average of that, we
obtain
∂ ūi ∂ ūi 1 ∂
+ ūj = σ ij − ρu′i u′j
∂t ∂xj ρ ∂xj
where, for simplicity, we have ignored gravity (the ∇Φ-term). This equation looks
identical to the Navier-Stokes equation (in absence of gravity), except for the −ρu′i u′j
term, which is what we call the Reynolds stress tensor:
82
Rij = −ρu′i u′j
Note that u′i u′j means the same averaging (time, space or ensemble) as above, but
now for the product of u′i and u′j . Note that ū′i = 0, by construction. However,
the expectation value for the product of u′i and u′j is generally not. As is evident
from the equation, the Reynolds stresses (which reflect momentum diffusion due
to turbulence) act in exactly the same way as the viscous stresses. However, they
are only present when the flow is turbulent.
Note also that the Reynolds stress tensor is related to the two-point correlation
tensor
83
• Turbulent flows have a high rate of viscous energy dissipation.
• Advected tracers are rapidly mixed by turbulent flows.
However, one further property of turbulence seems to be more fun-
damental than all of these because it largely explains why turbulence
demands a statistical treatment...turbulence is chaotic.
Turbulence kicks in at sufficiently high Reynolds number (typically R > 103 − 104 ).
Turbulent flow is characterized by irregular and seemingly random motion. Large
vortices (called eddies) are created. These contain a large amount of kinetic energy.
Due to vortex stretching these eddies are stretched thin until they ‘break up’ in
smaller eddies. This results in a cascade in which the turbulent energy is transported
from large scales to small scales. This cascade is largely inviscid, conserving the total
turbulent energy. However, once the length scale of the eddies becomes comparable
to the mean free path of the particles, the energy is dissipated; the kinetic energy
associated with the eddies is transformed into internal energy. The scale at which
this happens is called the Kolmogorov length scale. The length scales between
the scale of turbulence ‘injection’ and the Kolomogorov length scale at which it
is dissipated is called the inertial range. Over this inertial range turbulence is
believed/observed to be scale invariant. The ratio between the injection scale, L,
and the dissipation scale, l, is proportional to the Reynolds number according to
L/l ∝ R3/4 . Hence, two turbulent flows that look similar on large scales (comparable
L), will dissipate their energies on different scales, l, if their Reynolds numbers are
different.
84
CHAPTER 11
Sound Waves
If the perturbation is small, we may assume that the velocity gradients are so small
that viscous effects are negligble (i.e., we can set ν = 0). In addition, we assume that
the time scale for conductive heat transport is large, so that energy exchange due to
conduction can also safely be ignored. In the absence of these dissipative processes,
the wave-induced changes in gas properties are adiabatic.
Let (ρ0 , P0 , ~u0) be a uniform, equilibrium solution of the Euler fluid equations
(i.e., ignore viscosity). Also, in what follows we will ignore gravity (i.e., ∇Φ = 0).
Uniformity implies that ∇ρ0 = ∇P0 = ∇~u0 = 0. In addition, since the only al-
lowed motion is uniform motion of the entire system, we can always use a Galilean
coordinate transformation so that ~u0 = 0, which is what we adopt in what follows.
85
Substitution into the continuity and momentum equations, one obtains that ∂ρ0 /∂t =
∂~u0 /∂t = 0, indicative of an equilibrium solution as claimed.
Perturbation Analysis: Consider a small perturbation away from the above equi-
librium solution:
ρ0 → ρ0 + ρ1
P0 → P0 + P1
~u0 → ~u0 + ~u1 = ~u1
where |ρ1 /ρ0 | ≪ 1, |P1 /P0 | ≪ 1 and ~u1 is small (compared to the sound speed, to
be derived below).
Next we linearize these equations, which means we use that the perturbed values
are all small such that terms that contain products of two or more of these quantities
are always negligible compared to those that contain only one such quantity. Hence,
the above equations reduce to
∂ρ1
+ ρ0 ∇~u1 = 0
∂t
∂~u1 ∇P1
+ = 0
∂t ρ0
86
These equations describe the evolution of perturbations in an inviscid and uniform
fluid. As always, these equations need an additional equation for closure. As men-
tioned above, we don’t need the energy equation: instead, we can use that the
flow is adiabatic, which implies that P ∝ ργ .
where we have used (∂P/∂ρ)0 as shorthand for the partial derivative of P (ρ) at
ρ = ρ0 . And since the flow is isentropic, we have that the partial derivative is for
constant entropy. Using that P (ρ0 ) = P0 and P (ρ0 + ρ1 ) = P0 + P1 , we find that,
when linearized,
∂P
P1 = ρ1
∂ρ 0
Note that P1 6= P (ρ1 ); rather P1 is the perturbation in pressure associated with the
perturbation ρ1 in the density.
Taking the partial time derivative of the above continuity equation, and using that
∂ρ0 /∂t = 0, gives
∂ 2 ρ1 ∂~u1
2
+ ρ0 ∇ · =0
∂t ∂t
Substituting the above momentum equation, and realizing that (∂P/∂ρ)0 is a
constant, then yields
∂ 2 ρ1 ∂P
− ∇2 ρ1 = 0
∂t2 ∂ρ 0
87
with ~k the wavevector, k = |~k| = 2π/λ the wavenumber, λ the wavelength,
ω = 2πν the angular frequency, and ν the frequency.
To gain some insight, consider the 1D case: ρ1 ∝ ei(kx−ωt) ∝ eik(x−vp t) , where we have
defined the phase velocity vp ≡ ω/k. This is the velocity with which the wave
pattern propagates through space. For our perturbation of a compressible fluid, this
phase velocity is called the sound speed, cs . Substituting the solution ρ1 ∝ ei(kx−ωt)
into the wave equation, we see that
s
ω ∂P
cs = =
k ∂ρ s
where we have made it explicit that the flow is assumed to be isentropic. Note that
the partial derivative is for the unperturbed medium. This sound speed is sometimes
called the adiabatic speed of sound, to emphasize that it relies on the assumption
of an adiabatic perturbation. If the fluid is an ideal gas, then
s
kB T
cs = γ
µ mp
which shows that the adiabatic sound speed of an ideal fluid increases with temper-
ature.
We can repeat the above derivation by relaxing the assumption of isentropic flow,
and assuming instead that (more generally) the flow is polytropic. In that case,
P ∝ ρΓ , with Γ the polytropic index (Note: a polytropic EoS is an example of a
barotropic EoS). The only thing that changes is that now the sound speed becomes
s s
∂P P
cs = = Γ
∂ρ ρ
which shows that the sound speed is larger for a stiffer EoS (i.e., a larger value of Γ).
Note also that, for our barotropic fluid, the sound speed is independent of ω. This
implies that all waves move equally fast; the shape of a wave packet is preserved
88
as it moves. We say that an ideal (inviscid) fluid with a barotropic EoS is a non-
dispersive medium.
To gain further insight, let us look once more at the (1D) solution for our perturba-
tion:
• The eikx part describes a periodic, spatial oscillation with wavelength λ = 2π/k.
We will return to this in Chapter 15, when we discuss the Jeans stability criterion.
89
CHAPTER 12
Shocks
When discussing sound waves in the previous chapter, we considered small (linear)
perturbations. In this Chapter we consider the case in which the perturbations are
large (non-linear). Typically, a large disturbance results in an abrupt discontinuity
in the fluid, called a shock. Note: not all discontinuities are shocks, but all shocks
are discontinuities.
90
Figure 13: The steepening of a sound wave into a shock due to non-linearity (i.e.,
sound speed depends on density).
91
Mach Number: if v is the flow speed of the fluid, and cs is the sound speed, then
the Mach number of the flow is defined as
v
M=
cs
Note: simply accelerating a flow to supersonic speeds does not necessarily generate
a shock. Shocks only arise when an obstruction in the flow causes a deceleration of
fluid moving at supersonic speeds. The reason is that disturbances cannot propagate
upstream, so that the flow cannot ‘adjust itself’ to the obstacle because there is no
way of propagating a signal (which always goes at the sound speed) in the upstream
direction. Consequently, the flow remains undisturbed until it hits the obstacle,
resulting in a discontinuous change in flow properties; a shock.
Structure of a Shock: Fig. 14 shows the structure of a planar shock. The shock
has a finite, non-zero width (typically a few mean-free paths of the fluid particles),
and separates the ‘up-stream’, pre-shocked gas, from the ‘down-stream’, shocked gas.
For reasons that will become clear in what follows, it is useful to split the downstream
region in two sub-regions; one in which the fluid is out of thermal equilibrium, with
net cooling L > 0, and, further away from the shock, a region where the downstream
gas is (once again) in thermal equilibrium (i.e., L = 0). If the transition between
these two sub-regions falls well outside the shock (i.e., if x3 ≫ x2 ) the shock is said
to be adiabatic. In that case, we can derive a relation between the upstream (pre-
shocked) properties (ρ1 , P1 , T1 , u1) and the downstream (post-shocked) properties
(ρ2 , P2 , T2 , u2 ); these relations are called the Rankine-Hugoniot jump conditions.
Linking the properties in region three (ρ3 , P3 , T3 , u3) to those in the pre-shocked gas
is in general not possible, except in the case where T3 = T1 . In this case one may
consider the shock to be isothermal.
92
ρ1 P1 s ρ2 P2 ρ3 P3
h v2 T2
v1 T1 v3 T3
o
c
L=0 k L>0 L=0
x1 x2 x3
Figure 14: Structure of a planar shock.
ρ1 u1 = ρ2 u2
This equation describes mass conservation across a shock.
93
∂ ∂ ∂Φ
(ρ ux ) = − (ρ ux ux + P ) − ρ
∂t ∂x ∂x
Integrating this equation over V and ignoring any gradient in Φ across the shock, we
obtain
ρ1 u21 + P1 = ρ2 u22 + P2
This equation describes how the shock converts ram pressure into thermal
pressure.
Finally, applying the same to the energy equation under the assumption that the
shock is adiabatic (i.e., dQ/dt = 0), one finds that (E + P )u has to be the same on
both sides of the shock, i.e.,
1 2 P
u +Φ+ε+ ρ u = constant
2 ρ
We have already seen that ρ u is constant. Hence, if we once more ignore gradients
in Φ across the shock, we obtain that
1 2 1
u1 + ε1 + P1 /ρ1 = u22 + ε2 + P2 /ρ2
2 2
This equation describes how the shock converts kinetic energy into enthalpy.
Qualitatively, a shock converts an ordered flow upstream into a disordered (hot) flow
downstream.
The three equations in the rectangular boxes are known as the Rankine-Hugoniot
(RH) jump conditions for an adiabatic shock. Using straightforward but
tedious algebra, these RH jump conditions can be written in a more useful form
using the Mach number M1 of the upstream gas:
−1
ρ2 u1 1 γ −1 1
= = + 1−
ρ1 u2 M21 γ + 1 M21
P2 2γ γ−1
= M21 −
P1 γ+1 γ+1
T2 P2 ρ2 γ −1 2 2 1 4γ γ −1
= = γM1 − + −
T1 P1 ρ1 γ+1 γ+1 M21 γ−1 γ+1
94
Here we have used that for an ideal gas
kB T
P = (γ − 1) ρ ε = ρ
µ mp
Given that M1 > 1, we see that ρ2 > ρ1 (shocks compress), u2 < u1 (shocks
decelerate), P2 > P1 (shocks increase pressure), and T2 > T1 (shocks heat).
The latter may seem surprising, given that the shock is considered to be adiabatic:
although the process has been adiabatic, in that dQ/dt = 0, the gas has changed its
adiabat; its entropy has increased as a consequence of the shock converting kinetic
energy into thermal, internal energy. In general, in the presence of viscosity, a
change that is adiabatic does not imply that the states before and after are simply
linked by the relation P = K ργ , with K some constant. Shocks are always viscous,
which causes K to change across the shock, such that the entropy increases; it is this
aspect of the shock that causes irreversibility, thus defining an ”arrow of time”.
where we have used that γ = 5/3 for a monoatomic gas. Thus, with an adia-
batic shock you can achieve a maximum compression in density of a factor four!
Physically, the reason why there is a maximal compression is that the pressure and
temperature of the downstream fluid diverge as M21 . This huge increase in down-
stream pressure inhibits the amount of compression of the downstream gas. However,
this is only true under the assumption that the shock is adiabtic. The downstream,
post-shocked gas is out of thermal equilibrium, and in general will be cooling (i.e.,
L > 0). At a certain distance past the shock (i.e., when x = x3 in Fig. 14), the
fluid will re-establish thermal equilibrium (i.e., L = 0). In some special cases, one
can obtain the properties of the fluid in the new equilibrium state; one such case is
the example of an isothermal shock, for which the downstream gas has the same
temperature as the upstream gas (i.e., T3 = T1 ).
In the case of an isothermal shock, the first two Rankine-Hugoniot jump con-
95
ditions are still valid, i.e.,
ρ1 u1 = ρ3 u3
ρ1 u21 + P1 = ρ3 u23 + P3
However, the third condition, which derives from the energy equation, is no longer
valid. After all, in deriving that one we had assumed that the shock was adiabatic.
In the case of an isothermal shock we have to replace the third RH jump condition
with T1 = T3 . The latter implies that c2s = P3 /ρ3 = P1 /ρ1 , and allows us to rewrite
the second RH condition as
Here the second step follows from using the first RH jump condition. If we now
substitute this result back into the first RH jump condition we obtain that
2
ρ3 u1 u1
= = = M21
ρ1 u3 cs
Hence, in the case of isothermal shock (or an adiabatic shock, but sufficiently far
behind the shock in the downstream fluid), we have that there is no restriction to
how much compression the shock can achieve; depending on the Mach number of the
shock, the compression can be huge.
96
Figure 15: An actual example of a supernova blastwave. The red colors show the
optical light emitted by the supernova ejecta, while the green colors indicate X-ray
emission coming from the hot bubble of gas that has been shock-heated when the
blast-wave ran over it.
97
gas, which ‘pushes’ the shock outwards. As more and more material is swept-up, and
accelerated outwards, the mass of the shell increases, which causes the velocity of the
shell to decelerate. At the early stages, the cooling of the hot bubble is negligble, and
the blastwave is said to be in the adiabatic phase, also known as the Sedov-Taylor
phase. At some point, though, the hot bubble starts to cool, radiating away the
kinetic energy of the supernova, and lowering the interior pressure up to the point
that it no longer pushes the shell outwards. This is called the radiative phase.
From this point on, the shell expands purely by its inertia, being slowed down by the
work it does against the surrounding material. This phase is called the snow-plow
phase. Ultimately, the velocity of the shell becomes comparable to the sound speed
of the surrounding material, after which it continues to move outward as a sound
wave, slowly dissipating into the surroundings.
During the adiabatic phase, we can use a simple dimensional analysis to solve for
the evolution of the shock radius, rsh , with time. Since the only physical parameters
that can determine rsh in this adiabatic phase are time, t, the initial energy of the
SN explosion, ε0 , and the density of the surrounding medium, ρ0 , we have that
It is easy to check that there is only one set of values for η, α and β for which the
product on the right has the dimensions of length (which is the dimension of rsh .
This solution has η = 2/5, α = 1/5 and β = −1/5, such that
1/5
ε
rsh = A t2/5
ρ0
98
CHAPTER 13
Fluid Instabilities
99
time, τs , while re-establishing thermal equilibrium
proceeds much slower, on the conduction time, τc . Given that τs ≪ τc we can assume
that Pb∗ = P ′ , and treat the displacement as adiabatic. The latter implies that the
process can be described by an adiabatic EoS: P ∝ ργ . Hence, we have that
1/γ 1/γ 1/γ
Pb∗ P′ 1 dP
ρ∗b = ρb = ρb = ρb 1+ δz
Pb P P dz
In the limit of small displacements δz, we can use Taylor series expansion to show
that, to first order,
ρ dP
ρ∗b = ρ + δz
γ P dz
where we have used that initially ρb = ρ, and that the Taylor series expansion,
f (x) ≃ f (0)+f ′(0)x+ 21 f ′′ (0)x2 +..., of f (x) = [1+x]1/γ is given by f (x) ≃ 1+ γ1 x+....
Suppose we have a stratified medium in which dρ/dz < 0 and dP/dz < 0. In that
case, if ρ∗b > ρ′ the blob will be heavier than its surrounding and it will sink back to its
original position; the system is stable to convection. If, on the other hand, ρ∗b < ρ′
then the displacement has made the blob more buoyant, resulting in instability.
Hence, using that ρ′ = ρ + (dρ/dz) δz we see that stability requires that
dρ ρ dP
<
dz γ P dz
This is called the Schwarzschild criterion for convective stability.
It is often convenient to rewrite this criterion in a form that contains the temperature.
Using that
µ mp
ρ = ρ(P, T ) = P
kB T
it is straightforward to show that
dρ ρ dP ρ dT
= −
dz P dz T dz
Substitution in ρ′ = ρ + (dρ/dz) δz then yields that
∗ ′ 1 ρ dP ρ dT
ρb − ρ = −(1 − ) + δz
γ P dz T dz
100
Since stability requires that ρ∗b − ρ′ > 0, and using that δz > 0, dP/dz < 0 and
dT /dz < 0 we can rewrite the above Schwarzschild criterion for stability as
dT 1 T dP
< 1−
dz γ P dz
This shows that if the temperature gradient becomes too large the system becomes
convectively unstable: blobs will rise up until they start to loose their thermal en-
ergy to the ambient medium, resulting in convective energy transport that tries to
“overturn” the hot (high entropy) and cold (low entropy) material. In fact, without
any proof we mention that in terms of the specific entropy, s, one can also write
the Schwarzschild criterion for convective stability as ds/dz > 0.
dT 1 T dP
< 1−
dz γ P dz
dρ ρ dP
<
dz γ P dz
ds
>0
dz
It is easy to see where the RT instability comes from. Consider a fluid of density
ρ2 sitting on top of a fluid of density ρ1 < ρ2 in a gravitational field that is point-
ing in the downward direction. Consider a small perturbation in which the initially
horizontal interface takes on a small amplitude, sinusoidal deformation. Since this
101
Figure 16: Example of Rayleigh-Taylor instability in a hydro-dynamical simulation.
implies moving a certain volume of denser material down, and an equally large vol-
ume of the lighter material up, it is immediately clear that the potential energy of
this ‘perturbed’ configuration is lower than that of the initial state, and therefore
energetically favorable. Simply put, the initial configuration is unstable to small
deformations of the interface.
Stability analysis (i.e., perturbation analysis of the fluid equations) shows that the
dispersion relation corresponding to the RT instability is given by
r
g ρ2 − ρ1
ω = ±i k
k ρ2 + ρ1
where g is the gravitational acceleration, and the factor (ρ2 − ρ1 )/(ρ2 + ρ1 ) is called
the Atwood number. Since the wavenumber of the perturbation k > 0 we see
that ω is imaginary, which implies that the perturbations will grow exponentially
(i.e., the system is unstable). If ρ1 > ρ2 though, ω is real, and the system is stable
(perturbations to the interface propagate as waves).
102
Figure 17: Illustration of onset of Kelvin-Helmholtz instability
Stability analysis (i.e., perturbation analysis of the fluid equations) shows that the
the dispersion relation corresponding to the KH instability is given by
ωR (ρ1 u1 + ρ2 u2 )
=
k ρ1 + ρ2
and
ωI (ρ1 ρ2 )1/2
= (u1 − u2 )
k ρ1 + ρ2
Since the imaginary part is non-zero, except for u1 = u2 , we we have that, in principle,
any velocity difference across an interface is KH unstable. In practice, surface
tension can stabilize the short wavelength modes so that typically KH instability
kicks in above some velocity treshold.
103
The mode that will destroy the cloud has k ∼ 1/Rc , so that the time-scale for cloud
destruction is
1 Rc δ + 2
τKH ≃ ≃
ω cs,h (δ + 1)1/2
Assuming pressure equilibrium between cloud and ICM, and adopting the EoS of an
ideal gas, implies that ρh Th = ρc Tc , so that
1/2 1/2
cs,h T ρc
= h1/2 = 1/2 = (δ + 1)1/2
cs,c Tc ρh
Hence, one finds that the Kelvin-Helmholtz time for cloud destruction is
1 Rc δ + 2
τKH ≃ ≃
ω cs,c δ + 1
Note that τKH ∼ ζ(Rc/cs,c ) = ζτs , with ζ = 1(2) for δ ≫ 1(≪ 1). Hence, the Kelvin-
Helmholtz instability will typically destroy clouds falling into a hot ”atmosphere”
on a time scale between one and two sound crossing times, τs , of the cloud. Note,
though, that magnetic fields and/or radiative cooling at the interface may stabilize
the clouds.
104
and a Jeans mass 3
4 λJ π
MJ = πρ0 = ρ0 λ3J
3 2 6
From the dispersion relation one immediately sees that the system is unstable (i.e.,
ω is imaginary) if k < kJ (or, equivalently, λ > λJ or M > MJ ). This is called the
Jeans criterion for gravitational instability. It expresses when pressure forces
(which try to disperse matter) are no longer able to overcome gravity (which tries to
make matter collapse), resulting in exponential gravitational collapse on a time scale
r
3π
τff =
32 G ρ
In deriving the Jeans Stability criterion you will encounter a somewhat puzzling issue.
Consider the Poisson equation for the unperturbed medium (which has density ρ0
and gravitational potential Φ0 ):
∇2 Φ0 = 4πGρ0
105
Figure 18: The locus of ther-
mal equilibrium (L = 0) in
the (ρ, T ) plane, illustrating the
principle of thermal instability.
The dashed line indicates a line
of constant pressure.
The condition L(ρ, T ) = 0 corresponds to a curve in the (ρ, T )-plane with a shape
similar to that shown in Fig. 18. It has flat parts at T ∼ 106 K, at T ∼ 104K, at
T ∼ 10 − 100K. This can be understood from simple atomic physics (see for example
§ 8.5.1 of Mo, van den Bosch & White, 2010). Above the TE curve we have that
L > 0 (net cooling), while below it L < 0 (net heating). The dotted curve indicates
a line of constant pressure (T ∝ ρ−1 ). Consider a blob in thermal and mechanical
(pressure) equilibrium with its ambient medium, and with a pressure indicated by
the dashed line. There are five possible solutions for the density and temperature of
the blob, two of which are indicated by P1 and P2 ; here confusingly the P refers to
‘point’ rather than ‘pressure’. Suppose I have a blob located at point P2 . If I heat
the blob, displacing it from TE along the constant pressure curve (i.e., the blob is
assumed small enough that the sound crossing time, on which the blob re-established
mechanical equilibrium, is short). The blob now finds itself in the region where L > 0
(i.e, net cooling), so that it will cool back to its original location on the TE-curve;
the blob is stable. For similar reasons, it is easy to see that a blob located at point
P1 is unstable. This instability is called thermal instability, and it explains
why the ISM is a three-phase medium, with gas of three different temperatures
(T ∼ 106 K, 104 K, and ∼ 10 − 100 K) coexisting in pressure equilibrium. Gas at any
other temperature but in pressure equilibrium is thermally unstable.
106
It is easy to see that the requirement for thermal instability translates into
∂L
<0
∂T P
which is known as the Field criterion for thermal instability (after astrophysicist
George B. Field).
we thus see that a polytropic ideal gas must have that T ∝ ρΓ−1 . Substituting that
in the expression for the Jeans mass, we obtain that
3 3 4
MJ ∝ ρ 2 Γ−2 = ρ 2 (Γ− 3 )
Thus, we see that for Γ > 4/3 the Jeans mass will increase with increasing density,
while the opposite is true for Γ < 4/3. Now consider a system that is (initially)
larger than the Jeans mass. Since pressure can no longer support it against its own
gravity, the system will start to collapse, which increases the density. If Γ < 4/3,
the Jeans mass will becomes smaller as a consequence of the collapse, and now small
subregions of the system will find themselves having a mass larger than the Jeans
mass ⇒ the system will start to fragment.
If the collapse is adiabatic (i.e., we can ignore cooling), then Γ = γ = 5/3 > 4/3 and
there will be no fragmentation. However, if cooling is very efficient, such that while
the cloud collapses it maintains the same temperature, the EoS is now isothermal,
which implies that Γ = 1 < 4/3: the cloud will fragment into smaller collapsing
clouds. Fragmentation is believed to underly the formation of star clusters.
107
A very similar process operates related to the thermal instability. In the discussion
of the Field criterion we had made the assumption “the blob is assumed small
enough that the sound crossing time, on which the blob re-established mechanical
equilibrium, is short”. Here ‘short’ means compared to the cooling time of the cloud.
Let’s define the cooling length lcool ≡ cs τcool , where cs is the cloud’s sound speed and
τcool is the cooling time (the time scale on which it radiates away most of its internal
energy). The above assumption thus implies that the size of the cloud, lcloud ≪
lcool . As a consequence, whenever the cloud cools somewhat, it can immediately
re-establish pressure equilibrium with its surrounding (i.e., the sound crossing time,
τs = lcloud /cs is much smaller than the cooling time τcool = lcool /cs ).
Now consider a case in which lcloud ≫ lcool (i.e., τcool ≪ τs ). As the cloud cools,
it cannot maintain pressure equilibrium with its surroundings; it takes too long for
mechanical equilibrium to be established over the entire cloud. What happens is that
smaller subregions, of order the size lcool , will fragment. The smaller fragments will
be able to maintain pressure equilibrium with their surroundings. But as the small
cloudlets cool further, the cooling length lcool shrinks. To see this, realize that when T
drops this lowers the sound speed and decreases the cooling time; after all, we are in
the regime of thermal instability, so (∂L/∂T )P < 0. As a consequence, lcool = cs τcool
drops as well. So the small cloudlett soon finds itself larger than the cooling length,
and it in turn will fragment. This process of shattering continues until the cooling
time becomes sufficiently long and the cloudletts are no longer thermally unstable
(see McCourt et al., 2018, MNRAS, 473, 5407 for details).
108
Part II: Collisionless Dynamics
The following chapters give an elementary introduction into the rich topic of colli-
sionless dynamics. The main goal is to highlight how the lack of collisions among the
constituent particles give rise to a dynamics that differs remarkably from collisional
fluids. We also briefly discuss the theory or orbits, which are the building blocks of
collisionless systems, the Virial theorem, and the gravothermal catastrophe, which
is a consequence of the negative heat capacity of a gravitational system. Finally, we
briefly discuss interactions (‘collisions’) among collisionless systems.
Collisionless Dynamics is a rich topic, and one could easily devote an entire course
to it (for example the Yale Graduate Course ‘ASTR 518; Galactic Dynamics’). The
following chapters therefore only scratch the surface of this rich topic. Readers who
want to get more indepth information are referred to the following excellent text-
books
- Galactic Dynamics by J. Binney & S. Tremaine
- Galactic Nuclei by D. Merritt
- Galaxy Formation and Evolution by H.J. Mo, F. van den Bosch & S. White
109
CHAPTER 14
Consider a density distribution ρ(~x). What is the gravitational force F~g acting on
a particle of mass m at location ~x ? We can sum the small constributions δ F~g from
different regions ~x ′ ± d3~x ′ , given by
m δm(~x ′ ) ~x ′ − ~x ~x ′ − ~x
δ F~g (~x) = G ′ 2 ′
= Gm ′ 3
ρ(~x ′ )d3~x ′
|~x − ~x| |~x − ~x| |~x − ~x|
R
Adding up all the small contributions yields F~g (~x) = δ F~g (~x) ≡ m ~g (~x), where
Z
~x ′ − ~x
~g (~x) = G d3~x ′ ′ ρ(~x ′ )
|~x − ~x|3
is the gravitational field (i.e., the force per unit mass). Using that
~x ′ − ~x 1
= ∇x
|~x ′ − ~x|3 |~x ′ − ~x|
110
we can rewrite g(~x) as
Z Z
3 ′ 1 ′ Gρ(~x ′ )
~g (~x) = G d ~x ∇x ρ(~x ) = ∇x d3~x ′ ≡ −∇x Φ
|~x ′ − ~x| |~x ′ − ~x|
where in the last step we have defined the gravitational potential
Z
ρ(~x ′ )
Φ(~x) = −G d3~x ′ ′
|~x − ~x|
It can be shown, that the above expression is equivalent to what is known as the
Poisson equation:
∇2 Φ = 4π G ρ
For a derivation, see Section 3.2 of Astrophysical Fluid Dynamics by Clarke &
Carswell, or Section 2.1 of Galactic Dynamics by Binney & Tremaine.
In general, it is extremely complicated to solve the Poisson equation for Φ(~x) given
ρ(~x) [see Chapter 2 of Galactic Dynamics by Binney & Tremaine for a detailed
discussion]. However, under certain symmetries, solutions to the Poisson equation
are fairly straightforward. In particular, under spherical symmetry the general
solution to the Poisson equation is
Z r Z ∞
1 ′ ′2 ′ ′ ′ ′
Φ(r) = −4πG ρ(r ) r dr + ρ(r ) r dr
r 0 r
Note that the potential at r depends on the mass distribution outside of r. However,
if we now compute the gravitational force per unit mass
dΦ G M(r)
F~g (r) = − êr = − êr
dr r2
where Z r
M(r) ≡ 4π ρ(r ′ ) r ′2 dr
0
is the enclosed mass within r. This shows that the gravitational force does not
depend on the mass distribution outside of r.
111
Newton’s first theorem: a body that is inside a spherical shell of matter experi-
ences no net gravitational force from that shell. The equivalent in general relativity
is called Birkhoff’s theorem.
This is easily understood from the fact that the solid angles that extent from a
point inside a sphere to opposing directions have areas on the sphere that scale as r 2
(where r is the distance from the point to the sphere), while the gravitational force
per unit mass scales as r −2 . Hence, the gravitational forces from the two opposing
areas exactly cancel.
Escape velocity: the velocity needed for a particle or fluid element to escape to
infinity. Since E = v 2 /2 + Φ(~x), and escape requires E > 0, the escape velocity is
p
Vesc (~x) = 2 |Φ(~x)|
independent of the symmetry (or lack thereof) of the mass distribution.
Since gas cannot be on self-intersecting orbits, gas in disk galaxies generally orbits
on circular orbits. The measured rotation velocities therefore reflect the circular
velocities, which can be used to infer the enclosed mass as a function of radius. This
method is generaly used to infer the presence of dark matter halos surrounding
disk galaxies.
112
Consider a gravitational system consisting of N particles (e.g., stars, fluid elements).
The total energy of the system is E = K + W , where
P
N
1
Total Kinetic Energy: K= 2
mi vi2
i=1
P
N P
G mi mj
Total Potential Energy: W = − 21 |~
ri −~
rj |
i=1 j6=i
The latter follows from the fact that gravitational binding energy between a pair
of masses is proportional to the product of their masses, and inversely proportional
to their separation. The factor 1/2 corrects for double counting the number of pairs.
where rij = |~ri − ~rj |. In the continuum limit this simply becomes
Z
1
W = ρ(~x) Φ(~x) d3~x
2
One can show (see e.g., Binney & Tremaine 2008) that this is equal to the trace of
the Chandrasekhar Potential Energy Tensor
Z
∂Φ 3
Wij ≡ − ρ(~x) xi d ~x
∂xj
In particular,
3
X Z
W = Tr(Wij ) = Wii = − ρ(~x) ~x · ∇Φ d3~x
i=1
113
which is another, equally valid, expression for the gravitational potential energy in
the continuum limit.
2K + W = 0
Combining the virial equation with the expression for the total energy, E = K +W ,
we see that for a system that obeys the virial theorem
E = −K = W/2
If we assume, for simplicity, that all galaxies have equal mass then we can rewrite
this as
N N
1 X 2 G (Nm)2 1 X X 1
Nm vi − =0
N i=1 2 N 2 i=1 rij
j6=i
2 hv 2 i
M=
G h1/ri
with
114
XX 1 N
1
h1/ri =
N(N − 1) i=1 j6=i rij
G M2
W =−
rg
Using the relations above, it is clear that rg = 2/h1/ri. We can now rewrite the
above equation for M in the form
rg hv 2 i
M=
G
Hence, one can infer the mass of our cluster of galaxies from its velocity dispersion
and its gravitation radius. In general, though, neither of these is observable, and one
uses instead
2
Reff hvlos i
M =α
G
where vlos is the line-of-sight velocity, Reff is some measure for the ‘effective’ radius
of the system in question, and α is a parameter of order unity that depends on the
radial distribution of the galaxies. Note that, under the assumption of isotropy,
2
hvlos i = hv 2 i/3 and one can also infer the mean reciprocal pair separation from the
projected pair separations; in other words under the assumption of isotropy one can
infer α, and thus use the above equation to compute the total, gravitational mass of
the cluster. This method was applied by Fritz Zwicky in 1933, who inferred that
the total dynamical mass in the Coma cluster is much larger than the sum of the
masses of its galaxies. This was the first observational evidence for dark matter,
although it took the astronomical community until the late 70’s to generally accept
this notion.
115
For a self-gravitating fluid
N
X 1 1 3
K= mi vi2 = N m hv 2 i = N kB T
i=1
2 2 2
where the last step follows from the kinetic theory of ideal gases of monoatomic
particles. In fact, we can use the above equation for any fluid (including a collisionless
one), if we interpret T as an effective temperature that measures the rms velocity
of the constituent particles. If the system is in virial equilibrium, then
3
E = −K = − N kB T
2
which, as we show next, has some important implications...
Heat Capacity: the amount of heat required to increase the temperature by one
degree Kelvin (or Celsius). For a self-gravitating fluid this is
dE 3
= − N kB
C≡
dT 2
which is negative! This implies that by losing energy, a gravitational system
gets hotter!! This is a very counter-intuitive result, that often leads to confusion and
wrong expectations. Below we give three examples of implications of the negative
heat capacity of gravitating systems,
116
potential energy that becomes more negative). In order for the star to remain in
virial equilibrium its kinetic energy, which is proportional to temperature, has to
increase; the star’s energy loss results in an increase of its temperature.
In the Sun, hydrogen burning produces energy that replenishes the energy loss from
the surface. As a consequence, the system is in equilibrium, and will not contract.
However, once the Sun has used up all its hydrogren, it will start to contract and heat
up, because of the negative heat capacity. This continues until the temperature in
the core becomes sufficiently high that helium can start to fuse into heavier elements,
and the Sun settles in a new equilibrium.
Example 3: Core Collapse a system with negative heat capacity in contact with
a heat bath is thermodynamically unstable. Consider a self-gravitating fluid of ‘tem-
perature’ T1 , which is in contact with a heat bath of temperature T2 . Suppose the
system is in thermal equilibrium, so that T1 = T2 . If, due to some small disturbance,
a small amount of heat is tranferred from the system to the heat bath, the negative
heat capacity implies that this results in T1 > T2 . Since heat always flows from hot
to cold, more heat will now flow from the system to the heat bath, further increasing
the temperature difference, and T1 will continue to rise without limit. This run-away
instability is called the gravothermal catastrophe. An example of this instability
is the core collapse of globular clusters: Suppose the formation of a gravitational
system results in the system having a declining velocity dispersion profile, σ 2 (r) (i.e.,
σ decreases with increasing radius). This implies that the central region is (dynami-
cally) hotter than the outskirts. IF heat can flow from the center to those outskirts,
the gravothermal catastrophe kicks in, and σ in the central regions will grow with-
out limits. Since σ 2 = GM(r)/r, the central mass therefore gets compressed into
a smaller and smaller region, while the outer regions expand. This is called core
collapse. Note that this does NOT lead to the formation of a supermassive black
hole, because regions at smaller r always shrink faster than regions at somewhat
larger r. In dark matter halos, and elliptical galaxies, the velocity dispersion profile
is often declining with radius. However, in those systems the two-body relaxation
time is soo long that there is basically no heat flow (which requires two-body in-
teractions). However, globular clusters, which consist of N ∼ 104 stars, and have
a crossing time of only tcross ∼ 5 × 106 yr, have a two-body relaxation time of only
∼ 5 × 108 yr. Hence, heat flow in globular clusters is not negligible, and they can
(and do) undergo core collapse. The collapse does not proceed indefinitely, because
of binaries (see Galactic Dynamics by Binney & Tremaine for more details).
117
CHAPTER 15
In this chapter we consider collisionless fluids, such as galaxies and dark matter halos.
As discussed in previous chapters, their dynamics is governed by the Collisionless
Boltzmann equation (CBE)
df ∂f ∂f ∂Φ ∂f
= + vi − =0
dt ∂t ∂xi ∂xi ∂vi
By taking the velocity moment of the CBE (see Chapter 8), we obtain the Jeans
equations
∂ui ∂ui 1 ∂ σ̂ij ∂Φ
+ uj =− −
∂t ∂xj ρ ∂xj ∂xi
which are the equivalent of the Navier-Stokes equations (or Euler equations), but for
a collisionless fluid. The quantity σ̂ij in the above expression is the stress tensor,
defined as
σ̂ij = −ρ hwi wj i = −ρ(hvi vj i − hvi i hvj i)
In this chapter, we write a hat on top of the stress tensor, in order to distinguish it
from the velocity dispersion tensor given by
σ̂ij
σij2 = hvi vj i − hvi i hvj i = −
ρ
This notation may cause some confusion, but it is adapted here in order to be con-
sistent with the notation in standard textbooks on galactic dynamics. For the same
reason, in what follows we will write hvi i in stead of ui (also because ui was defined
as the velocity of a fluid element, but for a collisionless fluid the concept of a fluid
element is not defined).
As we have discussed in detail in Chapters 4 and 5, for a collisional fluid the stress
tensor is given by
σ̂ij = −ρσij2 = −P δij + τij
118
and therefore completely specified by two scalar quantities; the pressure P and the
shear viscosity µ (as always, we ignore bulk viscosity). Both P and µ are related to
ρ and T via constitutive equations, which allow for closure in the equations.
In the case of a collisionless fluid, though, no consistutive relations exist, and the
(symmetric) velocity dispersion tensor has 6 unknowns. As a consequence, the Jeans
equations do not form a closed set. Adding higher-order moment equations of the
CBE will yield more equations, but this also adds new, higher-order unknowns such
as hvi vj vk i, etc. As a consequence, the set of CBE moment equations never closes!
Note that σij2 is a local quantity; σij2 = σij2 (~x). At each point ~x it defines the
velocity ellipsoid; an ellipsoid whose principal axes are defined by the orthogonal
eigenvectors of σij2 with lengths that are proportional to the square roots of the
respective eigenvalues.
Since these eigenvalues are typically not the same, a collisionless fluid experiences
anisotropic pressure-like forces. In order to be able to close the set of Jeans equa-
tions, it is common to make certain assumptions about the symmetry of the fluid.
For example, a common assumption is that the fluid is isotropic, such that the
(local) velocity dispersion tensor is specified by a single quantity; the local velocity
dispersion σ 2 . Note, though, that if with this approach, a solution is found, the
solution may not correspond to a physical distribution function (DF) (i.e., in order
to be physical, f ≥ 0 everywhere). Thus, although any real DF obeys the Jeans
equations, not every solution to the Jeans equations corresponds to a physical DF!!!
As a worked out example, we now derive the Jeans equations under cylindrical
symmetry. We therefore write the Jeans equations in the cylindrical coordinate
system (R, φ, z). The first step is to write the CBE in cylindrical coordinates.
df ∂f ∂f ∂f ∂f ∂f ∂f ∂f
= + Ṙ + φ̇ + ż + v̇R + v̇φ + v̇z
dt ∂t ∂R ∂φ ∂z ∂vR ∂vφ ∂vz
vR = Ṙ ⇒ v̇R = R̈
vφ = Rφ̇ ⇒ v̇φ = Ṙφ̇ + Rφ̈
vz = ż ⇒ v̇z = z̈
∂Φ v2
v̇R = − ∂R + Rφ
v v
v̇φ = − R1 ∂R
∂Φ
+ RR φ
v̇z = − ∂Φ∂z
The Jeans equations follow from multiplication with vR , vφ , and vz and integrat-
ing over velocity space. Note that the cylindrical symmetry requires that all
derivatives with respect to φ vanish. The remaining terms are:
120
Z Z
∂f ∂ ∂(ρhvR i)
vR d3~v = vR f d3~v =
∂t ∂t ∂t
Z Z
∂f 3 ∂ ∂(ρhvR2 i)
vR2 d ~v = vR2 f d3~v =
∂R ∂R ∂R
Z Z
∂f ∂ ∂(ρhvR vz i)
vR vz d3~v = vR vz f d3~v =
∂z ∂z ∂z
Z 2 Z 2 Z
vR vφ ∂f 3 1 ∂(vR vφ f ) 3 ∂(vR vφ2 ) 3 hvφ2 i
d ~v = d ~v − f d ~v = −ρ
R ∂vR R ∂vR ∂vR R
Z Z Z
∂Φ ∂f 3 ∂Φ ∂(vR f ) 3 ∂vR 3 ∂Φ
vR d ~v = d ~v − f d ~v = −ρ
∂R ∂vR ∂R ∂vR ∂vR ∂R
Z 2 Z Z
vR vφ ∂f 3 1 ∂(vR2 vφ f ) 3 ∂(vR2 vφ ) 3 hv 2 i
d ~v = d ~v − f d ~v = −ρ R
R ∂vφ R ∂vφ ∂vφ R
Z Z Z
∂Φ ∂f 3 ∂Φ ∂(vR f ) 3 ∂vz 3
vR d ~v = d ~v − f d ~v = 0
∂z ∂vz ∂z ∂vz ∂vR
Working out the similar terms for the other Jeans equations we finally obtain the
Jeans Equations in Cylindrical Coordinates:
2
∂(ρhvR i) ∂(ρhvR2 i) ∂(ρhvR vz i) hvR i − hvφ2 i ∂Φ
+ + +ρ + =0
∂t ∂R ∂z R ∂R
These are 3 equations with 9 unknowns, which can only be solved if we make addi-
tional assumptions. In particular, one often makes the following assumptions:
∂
1 System is static ⇒ the ∂t
-terms are zero and hvR i = hvz i = 0.
2 Velocity dispersion tensor is diagonal ⇒ hvi vj i = 0 (if i 6= j).
3 Meridional isotropy ⇒ hvR2 i = hvz2i = σR2 = σz2 ≡ σ 2 .
121
Under these assumptions we have 3 unknowns left: hvφ i, hvφ2 i, and σ 2 , and the Jeans
equations reduce to
2
∂(ρσ 2 ) σ − hvφ2 i ∂Φ
+ρ + =0
∂R R ∂R
∂(ρσ 2 ) ∂Φ
+ρ =0
∂z ∂z
Since we now only have two equations left, the system is still not closed. If from
the surface brightness we can estimate the mass density, ρ(R, z), and hence (using
the Poisson equation) the potential Φ(R, z), we can solve the second of these Jeans
equations for the meridional velocity dispersion:
Z∞
2 1 ∂Φ
σ (R, z) = ρ dz
ρ ∂z
z
and the first Jeans equation then gives the mean square azimuthal velocity
hvφ2 i = hvφ i2 + σφ2 :
∂Φ R ∂(ρσ 2 )
hvφ2 i(R, z) 2
= σ (R, z) + R +
∂R ρ ∂R
Thus, although hvφ2 i is uniquely specified by the Jeans equations, we don’t know how
it splits in the actual azimuthal streaming, hvφ i, and the azimuthal dispersion,
σφ2 . Additional assumptions are needed for this.
————————————————-
122
A similar analysis, but for a spherically symmetric system, using the spherical coordi-
nate system (r, θ, φ), gives the following Jeans equations in Spherical Symmetry
∂(ρhvr i) ∂(ρhvr2 i) ρ 2 ∂Φ
+ + 2hvr i − hvθ2 i − hvφ2 i + ρ =0
∂t ∂r r ∂r
∂(ρhvθ i) ∂(ρhvr vθ i) ρ
+ + 3hvr vθ i + hvθ2 i − hvφ2 i cotθ = 0
∂t ∂r r
∂(ρhvφ i) ∂(ρhvr vφ i) ρ
+ + [3hvr vφ i + 2hvθ vφ icotθ] = 0
∂t ∂r r
If we now make the additional assumptions that the system is static and that also
the kinematic properties of the system are spherically symmetric, then there can
be no streaming motions and all mixed second-order moments vanish. Consequently,
the velocity dispersion tensor is diagonal with σθ2 = σφ2 . Under these assumptions
only one of the three Jeans equations remains:
∂(ρσr2 ) 2ρ 2 ∂Φ
+ σr − σθ2 + ρ =0
∂r r ∂r
Notice that this single equation still constains two unknown, σr2 (r) and σθ2 (r) (if we
assume that the density and potential are known), and can thus not be solved.
It is useful to define the anisotropy parameter
where the second equality only holds under the assumption that the kinematics are
spherically symmetric.
1 ∂(ρhvr2 i) βhvr2 i dΦ
+2 =−
ρ ∂r r dr
If we now use that dΦ/dr = GM(r)/r, we can write the following expression for the
123
enclosed (dynamical) mass:
rhvr2 i d ln ρ d lnhvr2 i
M(r) = − + + 2β
G d ln r d ln r
Hence, if we know ρ(r), hvr2 i(r), and β(r), we can use the spherical Jeans equation
to infer the mass profile M(r).
with Υ(r) the mass-to-light ratio. Similarly, the line-of-sight velocity disper-
sion, σp2 (R), which can be inferred from spectroscopy, is related to both hvr2 i(r) and
β(r) according to (see Fig. 19)
Z∞
ν r dr
Σ(R)σp2 (R) = 2 h(vr cos α − vθ sin α)2 i √
r 2 − R2
R
Z∞
ν r dr
= 2 hvr2 i cos2 α + hvθ2 i sin2 α √
r 2 − R2
R
Z∞
R2 ν hvr2 i r dr
= 2 1−β 2 √
r r 2 − R2
R
The 3D luminosity density is trivially obtained from the observed Σ(R) using the
Abel transform
Z∞
1 dΣ dR
ν(r) = − √
π dR R2 − r 2
r
In general, we have three unknowns: M(r) [or equivalently ρ(r) or Υ(r)], hvr2 i(r) and
β(r). With our two observables Σ(R) and σp2 (R), these can only be determined if we
make additional assumptions.
124
Figure 19: Geometry related to projection
EXAMPLE 1: Assume isotropy: β(r) = 0. In this case we can use the Abel
transform to obtain
Z∞
2 1 d(Σσp2 ) dR
ν(r)hvr i(r) = − √
π dR R2 − r 2
r
125
We can now use the spherical Jeans Equation to write β(r) in terms of M(r),
ν(r) and hvr2 i(r). Substituting this in the equation for Σ(R)σp2 (R) yields a solution
for hvr2i(r), and thus for β(r). As long as β(r) ≤ 1 the model is said to be self-
consistent within the context of the Jeans equations.
Almost always, radically different models (based on radically different assumptions)
can be constructed, that are all consistent with the data and the Jeans equations.
This is often referred to as the mass-anisotropy degeneracy. Note, however, that
none of these models need to be physical: they can still have f < 0.
126
CHAPTER 16
Orbit Theory
The ‘building-blocks’ of collisionless systems, such as galaxies and dark matter halos,
are orbits. In this Chapter, we very briefly highlight a few aspects of orbit theory.
A more detailed account of orbit theory in the context of astrophysical systems can
be found in the excellent textbook ”Galactic Dynamics” by Binney & Tremaine.
for any t1 and t2 . The value of the integral of motion can be the same for different
orbits. Note that an integral of motion can not depend on time. Orbits can have
from zero to five integrals of motion. If the Hamiltonian does not depend on time,
then energy is always an integral of motion.
127
that is one lower than that of the non-resonant regular orbits (i.e., lωi + mωj is an
extra isolating integral of motion). Orbits with fewer than n isolating integrals of
motion are called irregular or stochastic.
Every spherical potential admits at least four isolating integrals of motion, namely
energy, E, and the three components of the angular momentum vector L. ~ Orbits in
a flattened, axisymmetric potential frequently (but not always) admit three isolating
integrals of motion: E, Lz (where the z-axis is the system’s symmetry axis), and a
non-classical third integral I3 (the integral is called non-classical since there is no
analytical expression of I3 as function of the phase-space variables).
dI ∂I dxi ∂I dvi ∂I
= + = ~v · ∇I − ∇Φ · =0
dt ∂xi dt ∂vi dt ∂~v
Compare this to the CBE for a steady-state (static) system:
∂f
~v · ∇f − ∇Φ · =0
∂~v
Thus the condition for I to be an integral of motion is identical with the condition
for I to be a steady-state solution of the CBE. This implies the following:
Jeans Theorem: Any steady-state solution of the CBE depends on the phase-
space coordinates only through integrals of motion. Any function of these integrals is
a steady-state solution of the CBE.
~
Hence, the DF of any steady-state spherical system can be expressed as f = f (E, L).
2
If the system is spherically symmetric in all its properties, then f = f (E, L ), i.e.,
the DF can only depend on the magnitude of the angular momentum vector, not on
its direction.
128
An even simpler case to consider is the one in which f = f (E): Since E = Φ(~r) +
1 2
[v + vθ2 + vφ2 ] we have that
2 r
Z
2 1 2 1 2 2 2
hvr i = dvr dvθ dvφ vr f Φ + [vr + vθ + vφ ]
ρ 2
Z
2 1 2 1 2 2 2
hvθ i = dvr dvθ dvφ vθ f Φ + [vr + vθ + vφ ]
ρ 2
Z
2 1 2 1 2 2 2
hvφ i = dvr dvθ dvφ vφ f Φ + [vr + vθ + vφ ]
ρ 2
Since these equations differ only in the labelling of one of the variables of integration,
it is immediately evident that hvr2 i = hvθ2 i = hvφ2 i. Hence, assuming that f = f (E) is
identical to assuming that the system is isotropic (and thus β(r) = 0). And since
Z
1 1 2 2 2
hvi i = dvr dvθ dvφ vi f Φ + [vr + vθ + vφ ]
ρ 2
it is also immediately evident that hvr i = hvθ i = hvφ i = 0. Thus, a system with
f = f (E) has no net sense of rotation.
The more general f (E, L2 ) models typically are anisotropic. Models with 0 < β ≤ 1
are radially anisotropic. In the extreme case of β = 1 all orbits are purely radial
and f = g(E) δ(L), with g(E) some function of energy. Tangentially anisotropic
models have β < 0, with β = −∞ corresponding to a model in which all orbits are
circular. In that case f = g(E) δ[L − Lmax (E)], where Lmax (E) is the maximum
angular momentum for energy E. Another special case is the one in which β(r) = β
is constant; such models have f = g(E) L−2β .
129
2
∂(ρhvR2 i) ∂(ρhvR vz i) hvR i − hvφ2 i ∂Φ
+ +ρ + =0
∂R ∂z R ∂R
∂(ρhvR vz i) ∂(ρhvz2 i) hvR vz i ∂Φ
+ +ρ + =0
∂R ∂z R ∂z
which clearly doesn’t suffice to solve for the four unknowns (modelling three-integral
axisymmetric systems is best done using the Schwarzschild orbit superposition tech-
nique). To make progress with Jeans modeling, one has to make additional assump-
tions. A typical assumption is that the DF has the two-integral form f = f (E, Lz ).
In that case, hvR vz i = 0 [velocity ellipsoid now is aligned with (R, φ, z)] and hvR2 i =
hvz2 i (see Binney & Tremaine 2008), so that the Jeans equations reduce to
2
∂(ρhvR2 i) hvR i − hvφ2 i ∂Φ
+ρ + =0
∂R R ∂R
∂(ρhvz2 i) ∂Φ
+ρ =0
∂z ∂z
which is a closed set for the two unknowns hvR2 i (= hvz2 i) and hvφ2 i. Note, however,
that the Jeans equations provide no information about how hvφ2 i splits in streaming
and random motions. In practice one often writes that
hvφ i = k hvφ2 i − hvR2 i
with k a free parameter. When k = 1 the azimuthal dispersion is σφ2 ≡ hvφ2 i − hvφ i2 =
σR2 = σz2 everywhere. Such models are called oblate isotropic rotators.
130
CHAPTER 17
Consider an encounter between two collisionless N-body systems (i.e., dark matter
halos or galaxies): a perturber P and a system S. Let q denote a particle of S and
let b be the impact parameter, v∞ the initial speed of the encounter, and R0 the
distance of closest approach (see Fig. 20).
There are only two cases in which we can calculate the outcome of the encounter
analytically:
• high speed encounter (v∞ ≫ vcrit ). In this case the encounter is said to
be impulsive and one can use the impulsive approximation to compute its
outcome.
• large mass ratio (MP ≪ MS ). In this case one can use the treatment of
dynamical friction to describe how P loses orbital energy and angular mo-
mentum to S.
In all other cases, one basically has to resort to numerical simulations to study the
outcome of the encounter. In what follows we present treatments of first the impulse
approximation and then dynamical friction.
131
Figure 20: Schematic illustration of an encounter with impact parameter b between
a perturber P and a subject S.
In the large v∞ limit, we have that the distance of closest approach R0 → b, and the
velocity of P wrt S is vP (t) ≃ v∞~ey ≡ vP~ey . Consequently, we have that
~
R(t) = (b, vP t, 0)
132
Let ~r be the position vector of q wrt S and adopt the distant encounter approx-
imation, which means that b ≫ max[RS , RP ], where RS and RP are the sizes of S
and P , respectively. This means that we may treat P as a point mass MP , so that
GMP
ΦP (~r) = −
~
|~r − R|
~ we we have that
Using geometry, and defining φ as the angle between ~r and R,
• The first term on rhs is a constant, not yielding any force (i.e., ∇r ΦP = 0).
• The second term on the rhs describes how the center of mass of S changes its
velocity due to the encounter with P .
• The third term on the rhs corresponds to the tidal force per unit mass and is
the term of interest to us.
133
GMP 3 ′ 2 2 1 ′2
Φ3 (~r) = − 3 r cos φ − r
R 2 2
GMP ′2 1 ′2 1 ′2
= − 3 x − y − z
R 2 2
GMP 2
Fx = x(2 − 3 sin θ) − 3y sin θ cos θ
R3
GMP 2
Fy = y(2 − 3 cos θ) − 3x sin θ cos θ
R3
GMP
Fx = − 3 z
R
Using these, we have that
Z Z Z π/2
dvx dt
∆vx = dt = Fx dt = Fx dθ
dt −π/2 dθ
with similar expressions for ∆vy and ∆vz . Using that θ = tan−1 (vP t/b) one has that
dt/dθ = b/(vP cos2 θ). Substituting the above expressions for the tidal force, and
using that R = b/ cos θ, one finds, after some algebra, that
2GMP
∆~v = (∆vx , ∆vy , ∆vz ) = (x, 0, −z)
vP b2
Substitution in the expression for ∆ES yields
134
Z
1 2 G2 MP2
∆ES = |∆~v |2 ρ(r) d3~r = MS hx2 + z 2 i
2 vP2 b4
Under the assumption that S is spherically symmetric we have that hx2 + z 2 i =
2
3
hx2 + y 2 + z 2 i = 23 hr 2 i and we obtain the final expression for the energy increase of
S as a consequence of the impulsive encounter with P :
2
4 MP hr 2 i
∆ES = G2 MS
3 vP b4
This derivation, which is originally due to Spitzer (1958), is surprisingly accurate for
encounters with b > 5max[RP , RS ], even for relatively slow encounters with v∞ ∼ σS .
For smaller impact parameters one has to make a correction (see Galaxy Formation
and Evolution by Mo, van den Bosch & White 2010 for details).
The impulse approximation shows that high-speed encounters can pump energy
into the systems involved. This energy is tapped from the orbital energy of the two
systems wrt each other. Note that ∆ES ∝ b−4 , so that close encounters are far more
important than distant encounters.
After the encounter, S has gained kinetic energy (in the amount of ∆ES ), but
its potential energy has remained unchanged (recall, this is the assumption that
underlies the impulse approximation). As a consequence, after the encounter S will
no longer be in virial equilibrium; S will have to readjust itself to re-establish
virial equilibrium.
135
Let K0 and E0 be the initial (pre-encounter) kinetic and total energy of S. The
virial theorem ensures that E0 = −K0 . The encounter causes an increase of
(kinetic) energy, so that K0 → K0 + ∆ES and E0 → E0 + ∆ES . After S has re-
established virial equilibrium, we have that K1 = −E1 = −(E0 + ∆ES ) = K0 − ∆ES .
Thus, we see that virialization after the encounter changes the kinetic energy of
S from K0 + ∆ES to K0 − ∆ES ! The gravitational energy after the encounter
is W1 = 2E1 = 2E0 + 2∆ES = W0 + 2∆ES , which is less negative than before
the encounter. Using the definition of the gravitational radius (see Chapter 18),
rg = GMS2 /|W |, from which it is clear that the (gravitational) radius of S increases
due to the impulsive encounter. Note that here we have ignored the complication
coming from the fact that the injection of energy ∆ES may result in unbinding some
of the mass of S.
Although these views are similar, there are some subtle differences. For example,
according to the first two descriptions dynamical friction is a local effect. The
136
third description, on the other hand, treats dynamical friction more as a global
effect. As we will see, there are circumstances under which these views make different
predictions, and if that is the case, the third and latter view presents itself as the
better one.
Chandrasekhar derived an expression for the dynamical friction force which, although
it is based on a number of questionable assumptions, yields results in reasonable
agreement with simulations. This so-called Chandrasekhar dynamical friction
force is given by
Here ρ(< vS ) is the density of particles of mass m that have a speed vm < vS , and ln Λ
is called the Coulomb logarithm. It’s value is uncertain (typically 3 ∼ < ln Λ < 30).
∼
One often approximates it as ln Λ ∼ ln(Mh /MS ), where Mh is the total mass of
the system of particles of mass m, but this should only be considered a very rough
estimate at best. The uncertainties for the Coulomb logarithm derive from the
oversimplified assumptions made by Chandrasekhar, which include that the medium
through which the subject mass is moving is infinite, uniform and with an isotropic
velocity distribution f (vm ) for the sea of particles.
Similar to frictional drag in fluid mechanics, F~df is always pointing in the direction
opposite of vS .
Note that F~df is independent of the mass m of the constituent particles, and pro-
portional to MS2 . The latter arises, within the second or third view depicted above,
from the fact that the wake or response density has a mass that is proportional to
MS , and the gravitational force between the subject mass and the wake/response
density therefore scales as MS2 .
137
Figure 21: Examples of the response density in a host system due to a perturber
orbiting inside it. The back-reaction of this response density on the perturber causes
the latter to experience dynamical friction. The position of the perturber is indicated
by an asterisk. [Source: Weinberg, 1989, MNRAS, 239, 549]
Let us assume that the host mass is a singular isothermal sphere with density
and potential given by
Vc2
ρ(r) = 2
Φ(r) = Vc2 ln r
4πGr
2
where Vc = GMh /rh with rh the radius of the host mass. If we further assume
that this host mass has, at each point, an isotropic and Maxwellian velocity
distrubution, then
138
2
ρ(r) vm
f (vm ) = exp − 2
(2πσ 2 )3/2 2σ
√
with σ = Vc / 2.
139
that virialized dark matter halos all have the same average density). Using that
ln Λ ∼ ln(Mh /MS ) and assuming that the subject mass starts out from an initial
radius ri = rh , we obtain a dynamical friction time
Mh /MS
tdf = 0.12 tH
ln(Mh /MS )
Hence, the time tdf on which dynamical friction brings an object of mass MS moving
in a host of mass Mh from an initial radius of ri = rh to r = 0 is shorter than the
Hubble time as long as MS ∼ > M /30. Hence, dynamical friction is only effective for
h
fairly massive objects, relative to the mass of the host. In fact, if you take into
account that the subject mass experiences mass stripping as well (due to the tidal
interactions with the host), the dynamical friction time increases by a factor 2 to 3,
and tdf < tH actually requires that MS ∼> M /10.
h
140
Part III: Radiative Processes
With the exception of meteorites, neutrinos, gravitational waves, and cosmic rays, all
information about the Universe reaches us in the form of radiation. Understanding
how radiation is produced, and how it interacts with matter on its way from the
source to our telescopes is therefore of crucial importance for astrophysics. The
following chapters give an elementary introduction into these topics.
Radiative processes is a rich topic, and one could easily devote an entire course to it.
The following chapters therefore only scratch the surface of this rich topic. Readers
who want to get more indepth information are referred to the following excellent
textbooks
- Radiative Processes in Astrophysics by G. Rybicki & A. Lightman
- Astrophysics: Decoding the Cosmos by J. Irwin
- The Physics of Astrophysics I. Radiation by F. Shu
- Theoretical Astrophysics by M. Bartelmann
141
CHAPTER 18
Radiation Essentials
Note that [Lν ] = erg s−1 Hz−1 , while [L] = erg s−1 .
Flux: The flux, f , of a source is the radiation energy per unit time passing through
a unit area
dL = f dA [f ] = erg s−1 cm−2
where A is the area. Similarly, we can also define the spectral flux density (or
simply ‘flux density’), as the flux per unit spectral bandwidth:
dLν = fν dA [fν ] = erg s−1 cm−2 Hz−1
In radio astronomy, one typically expresses fν in Jansky, where 1Jy = 10−23 erg s−1 cm−2 Hz−1 .
As with the SEDs, one may also express spectral flux densities as fλ . Using that
λ = c/ν, and using that fν dν = fλ dλ one has that
λ2 ν2
fν = fλ , fλ = fν
c c
142
Figure 22: Diagrams showing intensity and its dependence on direction and solid
angle. Fig. (a) depicts the ‘observational view’, where dA represents an element of
a detector. The arrows show incoming rays from the center of the source. Fig. (b)
depicts the ‘emission view’, where dA represents the surface of a star. At each point
on the surface, photons leave in all directions away from the surface.
where θ is the angle between the normal of the surface area through which the flux
is measured and the direction of the solid angle. The unit of intensity is [I] =
erg s−1 cm−2 sr−1 . Here ‘sr’ is a steradian, which is the unit of solid angle measure
(there are 4π steradians in a complete sphere). As with the flux and luminosity,
one can also define a specific intensity, Iν , which is the intensity per unit spectral
bandwidth ([Iν ] = erg s−1 cm−2 Hz−1 sr−1 ).
The flux emerging from the surface of a star with luminosity L and radius R∗ is
Z Z 2π Z π/2
L
F ≡ = I cos θ dΩ = dφ dθ I cos θ sin θ = π I
4πR∗2 half sphere 0 0
where we have used that dΩ = sin θ dθ dφ, and the fact that the integration over the
solid angle Ω is only to be performed over half a sphere. Note that an observer can
only measure the surface brightness of resolved objects; if unresolved, the observer
can only measure the objects flux.
143
is given by Z Z
f= I(Ω) cos θ dΩ ≃ I(Ω) dΩ ≡ hIi ΩS
ΩS
where we have assumed that ΩS is small, so that variations of cos θ across the object
can be neglected. Since both f ∝ r −2 and ΩS ∝ r −2 , where r is the object’s distance,
we see that the average surface brightness hIi is independent of distance.
Energy density: the energy density, u, is a measure of the radiative energy per unit
volume (i.e., [u] = erg cm−3 ). If the radiation intensity as seen from some specific
location in space is given by I(Ω), then the energy density at that location is
Z
1 4π
u= I dΩ ≡ J
c c
where Z
1
J≡ I dΩ
4π
is the mean intensity (i.e., average over 4π sterradian). If the radiation is isotropic
(i.e., the center of a star, or, to good approximation, a random location in the early
Universe), then J = I. If the radiation intensity
P is due to the summed intensity from
1
a number of individual sources, then u = c i fi , where fi is the flux due to source
i.
Recall from Chapter 6 that the number density of photons emerging from a Black
Body of temperature T is given by
8π ν 2 dν
nγ (ν, T ) dν = 3 hν/k
c e BT − 1
8π h ν 3 dν
u(ν, T ) dν = nγ (ν, T ) hν dν =
c3 ehν/kB T − 1
Using that u(ν, T ) = (4π/c)Jν (T ) we have that the mean specific intensity from
a black body [for which one typically uses the symbol Bν (T )] is given by
2 h ν3 dν
Bν (T ) dν = 2 hν/k
c e BT − 1
144
Figure 23: Various Planck curves for different temperatures, illustrating Wien’s dis-
placement law. Note how the Planck curve for a black body with the temperature of
the Sun peaks at the visible wavelengths, where the sensitivity of our eyes is maximal
which is called the Planck curve (or ‘formula’). Integrating over frequency yields
the total, mean intensity emitted from the surface of a Black Body
Z ∞
σSB 4
J = J(T ) = Bν (T ) dν = T
0 π
where σSB is the Stefan-Boltzmann constant. This implies an energy density
4π 4σSB 4
u = u(T ) = J= T ≡ ar T 4
c c
where ar ≃ 7.6 × 10−15 erg cm−3 K−4 is called the radiation constant (see also
Chapter 6).
Wien’s Displacement Law: When the temperature of a Black Body emitter in-
creases, the overall radiated energy increases and the peak of the radiation curve
moves to shorter wavelengths. It is straightforward to show that the product of the
temperature and the wavelength at which the Planck curve peaks is a constant, given
by
λmax T = 0.29
145
where T is the absolute temperature, expressed in degrees Kelvin, and λmax is ex-
pressed in cm. This relation is called Wien’s Displacement Law.
which is known as the Stefan-Boltzmann law. This law is used to define the
effective temperature of an emitter.
L = 4 π R2 σSB Teff
4
where R is the emitter’s radius. The effective temperature is sometimes also called
the radiation temperature, as a measure for the temperature associated with the
radiation field.
Here FX (λ) describes the transmission of the filter that defines waveband X, R(λ)
is the transmission efficiency of the telescope + instrument, and T (λ) describes the
transmission of the atmosphere. The combined effect of FX , R, and T is typically
‘calibrated’ using standard stars with known fλ .
146
Magnitudes: For historical reasons, the flux of an astronomical object in waveband
X is usually quoted in terms of apparent magnitude:
fX
mX = −2.5 log
fX,0
where the flux zero-point fX,0 has traditionally been taken as the flux in the X
band of the bright star Vega. In recent years it has become more common to use
‘AB-magnitudes’, for which
Z
−20 −1 −2 −1
fX,0 = 3.6308 × 10 erg s cm Hz FX (c/ν) dν
Similarly, the luminosities of objects (in waveband X) are often quoted as an abso-
lute magnitude:
MX = −2.5 log(LX ) + CX
where CX is a zero point. It is usually convenient to write LX in units of the solar
luminosity in the same band, L⊙X , so that
LX
MX = −2.5 log + M⊙X ,
L⊙X
where M⊙X is the absolute magnitude of the Sun in the waveband in consideration.
Using the relation between luminosity and flux we have that
mX − MX = 5 log(r/r0 )
where r0 is a fiducial distance at which mX and MX are defined to have the same
value. Conventionally, r0 is chosen to be 10 pc.
147
CHAPTER 19
Most of the baryonic matter in the Universe is in a gaseous state, made up of ∼ 75%
Hydrogen (H), ∼ 25% Helium (He) and only small amounts of other elements (called
‘metals’). Gases can be neutral, ionized, or partially ionized. The degree of ionization
of any given element is specified by a Roman numeral after the element name. For
example HI and HII refer to neutral and ionized hydrogen, respectively, while CI is
neutral carbon, CII is singly-ionized carbon, and CIV is triply-ionized carbon. A gas
that is highly (largely) ionized is called a plasma.
In a system in TE the energy in the radiation field is in equilibrium with the kinetic
energy of the particles. If the system is isolated (to both matter and radiation), and
in mechanical equilibrium, then over time TE will be established. For a gas in TE,
the radiation temperature, TR , is equal to the kinetic temperature, T , is equal
to the excitation temperature, Tex (see below for definitions). Since no photons
are allowed to escape from a system in TE (this would correspond to energy loss,
and thus violate TE), the photons that are produced in the gas (represented by TR )
therefore must be tightly coupled to the random motions of the particles (represented
by T ). This coupling implies that, for bound states, the excitation and de-excitation
of the atoms and ions must be dominated by collisions. In other words, the collision
timescales must be shorter than the timescales associated with photon interactions
148
or spontaneous de-excitations. If that were not the case, T could not remain
equal to TR .
Local Thermodynamic Equilibrium (LTE): True TE is rare (in almost all cases
energy does escape the system in the form of radiation, i.e., the system cools), and
often temperature gradients are present. A good, albeit imperfect, example of TE
is the Universe as a whole prior to decoupling. Although true TE is rare, in many
systems (stars, gaseous spheres, ISM), we may apply local TE (LTE), which implies
that the gas is in TE, but only locally. In a system in LTE, there will typically be
gradients in temperature, density, pressure, etc, but they are sufficiently small over
the mean-free path of a gas particle. For stellar interiors, the fact that radiation is
‘locally trapped’ explains why it takes so long for photons to diffuse from the center
(where they are created in nuclear reactions) to the surface (where they are emitted
into space). In the case of the Sun, this timescale is of the order of 200,000 years.
where m is the mass of a gas particle. The temperature T is called the kinetic
temperature, and is related to the mean-square particle speed, hv 2i, according to
1 3
m hv 2i = kB T
2 2
p
The most probable speed of the Maxwell-Boltzmann
p distribution is vmp = 2kB T /m,
while the mean speed is vmean = 8kB T /m.
149
Setting up a Maxwellian velocity distribution requires many elastic collisions. In
the limit where most collisions are elastic, the system will typically very quickly equi-
librate to thermal equilibrium. However, depending on the temperature of the gas,
collisions can also be inelastic: examples of the latter are collisional excitations,
in which part of the kinetic energy of the particle is used to excite its target to an
excited state (i.e., the kinetic energy is now temporarily stored as a potential energy).
If the de-excitation is collisional, the energy is given back to the kinetic energy of
the gas (i.e., no photon is emitted). However, if the de-excitation is spontaneous
or via stimulated emission, a photon is emitted. If the gas is optically thin to the
emitted photon, the energy will escape the gas. The net outcome of the collisional
excitation is then one of dissipation, i.e., cooling (kinetic energy of the gas is being
radiated away). In the optically thick limit, the photon will be absorbed by another
atom (or ion). If the system is in (L)TE, the radiation temperature of these
photons being emitted and absorbed is equal to the kinetic temperature of the
gas. Although in LTE a subset of the collisions are inelastic, this subset is typically
small, and one may still use the Maxwell-Boltzmann distribution to characterize the
velocities of the gas particles.
150
Statistical Equilibrium: A system is said to be in statistical equilibrium if the
level populations of its constituent atoms and ions do not change with time (i.e., if
the transition rate into any given level equals the rate out).
Here Ni is the number of atoms in which electrons are in energy level i (here i reflects
the principal quantum number, with i = 1 refering to the ground state), νij is the
frequency corresponding to the energy difference ∆Eij = h νij of the energy levels, h
is the Planck constant, and gi is the statistical weight of population i.
Statistical weight: the statistical weight (aka ‘degeneracy’) indicates the number
of states at a given principal quantum number, n. In quantum mechanics, a total
of four quantum numbers is needed to describe the state of an electron: the prin-
ciple quantum number, n, the orbital angular momentum quantum number, l, the
magnetic quantum number ml , and the electron spin quantum number, ms . For
hydrogen l can take on the values 0, 1, ..., n − 1, the quantum number ml can take on
the values −l, −l + 1, −1, 0, 1, ..., l, and ms can take on the values +1/2 (‘up’) and
−1/2 (‘down’). Hence, gn = 2 n2 .
Excitation Temperature: The above Boltzmann law defines the excitation tem-
perature, Tex , as the temperature which, when put into the Boltzmann law for level
populations, results in the observed ratio of Nj /Ni . For a gas in (local) TE, all levels
in an atom can be described by the same Tex , which is also equal to the kinetic
temperature, T , and the radiation temperature, TR . Under non-LTE conditions,
each pair of energy levels can have a different Tex .
151
Figure 24: Illustration of the energy levels of the hydrogen atom and some of the
most important transitions.
is called the partition function. Note that at low temperatures, below that needed
to put a significant fraction of atoms in the first excited state, the partition function
becomes equal to the statistical weight of the ground state: U = g1 exp −0/kB T =
g1 = 2. After all, the exponential factors for the excited states are all extremely
small. (Note: the summation is over all principal quantum numbers up to some
maximum nmax , which is required to prevent divergence; see Irwin 2007 or Rybicki
& Lightmann 1979 for details).
Saha equation: the Saha equation expresses the ratios of atoms/ions in different
ionization states. In particular, the number of atoms/ions in the (K + 1)th ionization
152
Figure 25: The ionization fraction of hydrogen as a function of temperature, computed
using the Saha equation.
where Un is the partition function of the nth ionization state, ne is the electron
number density, and χK is the energy required to remove an electron from the ground
state of the K th ionization state. Note that unlike the Boltzmann equation, the
Saha equation for the ionization fractions has a dependence on electron density. This
reflects that if there is a higher density of free electrons, there is a greater probability
that an electron will recombine with the ion, lowering the ionization state of the gas.
Hydrogen Since Hydrogen is the most common element in the Universe, it is impor-
tant to have some understanding of its structure. Fig. 18 shows the energy levels of
a hydrogen atom and some if its most important transitions. The Balmer transitions
have wavelengths that fall in the optical, and are therefore well known to (optical)
astronomers. The Lyman lines typically fall in the UV, and can only be observed
from space (or for high-redshift objects, for which the rest-frame UV is redshifted
into the optical).
Using the Saha equation, we can compute the ionization fraction for hydrogen
as a function of temperature and electron density. The ionization fraction is x =
153
NHII /[NHI + NHII ], which can be computed using the Saha equation:
3/2
NHII 15 T 1.58 × 105
= 2.41 × 10 exp −
NHI ne T
where we have used that χHI = 13.6eV, UHI = 2 and UHII = 1 (i.e., the ionized
hydrogen atom is just a free proton and only exists in a single state). Figure . 24
shows the ionization fraction x as function of temperature. Note that the transition
from almost neutral to almost completely ionized is extremely rapid! We can also
compute, using the Boltzmann law, the ratio of hydrogen atoms in the first excited
state compared to those in the ground state. The latter shows that the temperature
must exceed 30,000K in order for there to be an appreciable (10 percent) number of
hydrogen atoms in the first excited state. However, at such high temperatures, one
typically has that most of the hydrogen will be ionized (unless the electron density
is unrealistically high). This somewhat unintuitive result arises from the fact that
there are many more possible states available for a free electron than for a bound
electron in the first excited state. In conclusion, neutral hydrogen in LTE will have
virtually all its atoms in the ground state.
If densities are sufficiently low, which is typically the case in the ISM, the collisional
excitation rate (which scales with n2e ) is lower than the spontaneous de-excitation
rate. If that is the case, the hydrogen gas is no longer in LTE, and once again,
virtually all neutral hydrogen will find itself in the ground state. Neutral hydrogen
in the ISM is in the ground-state and is typically NOT in LTE. One important con-
sequence of the fact that neutral hydrogen is basically always observed in the ground
state, is that observations of hydrogen emission lines (Lyman, Balmer, Paschen, etc)
indicates that the hydrogen must be ionized; the lines arise from recombinations,
and are therefore called recombination lines. Note that Balmer lines are often
observed in absorption (for example, Balmer lines are evident in a spectrum of the
Sun). This implies that their must be hydrogen present in the first excited state,
which seems at odds with the conclusions reached above. A small fraction of excited
hydrogen atoms, though, can still produce strong absorption lines, simply because
hydrogen is so abundant.
21cm line emission: The ground state of hydrogen is split into two hyperfine
states due to the two possible orientations of the proton and electron spins: The
state in which the spins of proton and electron are aligned has slightly higher energy
than the one in which they are anti-aligned. The energy difference between these
154
two hyperfine states corresponds to a photon with a wavelength of 21cm (which
falls in the radio). The excitation temperature of this spin-flip transition is called
the spin temperature. Since, for typical interstellar medium conditions, the spin-
flip transition is collisionally induced (i.e., the rate for spontaneous de-excitation is
extremely low), the spin-temperature is typically equal to the kinetic temperate of the
hydrogen gas. The 21cm line is an important emission line to probe the distribution
(and temperature) of neutral hydrogen gas in the Universe.
155
CHAPTER 20
• light echos
• polarization
Scattering interactions are categorized as either elastic (coherent), where the pho-
ton energy is unchanged by the scattering event, or inelastic (incoherent), where
the photon energy changes.
156
Elastic scattering comes in three forms:
• Thomson scattering γ + e → γ + e
• Resonant scattering γ + X → X + → γ + X
• Rayleigh scattering γ + X → γ + X
• Compton scattering γ + e → γ ′ + e′
• Fluorescence γ + X → X ++ → γ ′ + X + → γ ′ + γ ′′ + X
Here accents indicate that the particle has a different energy (i.e., γ ′ is a photon with
a different energy than γ), and X ++ indicates a higher-excited state of X than X + .
8π 2 8πe4
σs = σT = re = 2 4
≃ 6.65 × 10−25 cm2
3 3me c
157
Figure 26: Illustration of how Thomson scattering causes polarization in the direc-
tions perpendicular to that of the incoming EM radiation. The incoming EM wave
~
causes the electron to oscillate in the direction of the oscillation of the E-field. This
acceleration of the electrical charge results in the emission of dipolar EM radiation.
158
Figure 27: The Klein-Nishina cross section for Compton scattering. As long as
hν ≪ me c2 one is in the Thomson scattering regime, and σs = σT . However, once
the photon energy becomes comparable to the rest-mass energy of the electron, Comp-
ton scattering takes over, and the cross-section (now called the Klein-Nishina cross-
section), starts to drop as ν −1 .
where θ is the angle between in incident and outgoing photon. This can also be
written as
λ′ − λ = λC (1 − cos θ)
which expresses that Compton scattering increases the wavelength of the photon by
of order the Compton wavelength λC = h/(me c) ∼ 2.43 × 10−10 cm. If λ ≫ λC
such a shift is negligble, and we are in the regime that is well described by Thomson
scattering.
159
So far we have considered the scattering of photons off of electrons at rest. A more
realistic treatment takes into account that electrons are also moving, and may do so
relativistically. This adds the possibility of the electron giving some of its kinetic
energy to the photon, which results in Inverse Compton (IC) scattering.
Whether the photon loses (Compton scattering) or gains (IC scattering) energy de-
pends on the energies of the photon and electron. Without derivation, the average
energy change of the photon per Compton scattering against electrons of temperature
Te = me hv 2 i/(3 kB) is
∆Eγ 4 kB Te − h ν
=
Eγ me c2
is the Lorentz factor. Hence, for ultra-relativistic electrons, which have a large
Lorentz factor, the frequency boost of a single IC scattering event can be enormous.
It is believed that this process, upscattering of low energy photons by the IC effect,
is at work in Active Galactic Nuclei.
160
terms of the Compton-y parameter) is measure for the electron pressure Pe ∝ ne Te
along the line-of-sight through the cluster. Observations of the SZ effect provide a
nearly redshift-independent means of detecting galaxy clusters.
Now consider the case of an EM wave of angular frequency, ω, interacting with the
atom/ion. The result is a forced, damped, harmonic oscillator, whose effective
cross section is given by
ω4
σs (ω) = σT
(ω 2 − ω02 )2 + (ω03 τe )2
(see Rybicki & Lightmann 1979 for a derivation). We can distinguish three regimes:
ω ≃ ω0 In this case
σT (Γcl /2)
σs (ω) ≃
2τe (ω − ω0 )2 + (Γcl /2)2
which corresponds to resonant scattering, in which the cross section is hugely
boosted wrt the Thomson case. NOTE: for resonant scattering to be important, it
is crucial that spontaneous de-excitation occurs before collisional excitation or de-
excitation (otherwise the photon energy is lost, and we are in the realm of absorption,
rather than scattering). Typically, this requires sufficiently low densities.
161
Figure 28: Illustration of how Rayleigh scattering causes the sky to be blue. Because
of its strong (λ−4 ) wave-length dependence, blue light is much more scattered than
red light. This causes the Sun light to appear redder than it really is, an effect that
strengthens when the path length through the atmosphere is larger (i.e., at sunrise
and sunset). The blue light is typically scattered multiple times before hitting the
observer, so that it appears to come from random directions on the sky.
ω ≪ ω0 In this case
4
ω
σs (ω) ≃ σT
ω0
which corresponds to Rayleigh scattering, which is characterized by a strong wave-
length dependence for the effective cross section of the form σs ∝ σT λ−4 .
Rayleigh scattering results from the electric polarizability of the particles. The
oscillating electric field of a light wave acts on the charges within a particle, causing
them to move at the same frequency (recall, the forcing frequency in this case is
much smaller than the natural frequency). The particle therefore becomes a small
radiating dipole whose radiation we see as scattered light.
162
We now turn our attention to a quantum-mechanical view of resonant scatter-
ing. The main difference between the classical view (above) and the quantum view
(below), is that in the latter there is not one, but many ‘natural frequencies’, νij , cor-
responding to all the possible energy-level-transitions ∆Eij = hνij that correspond
to the atom/ion in question.
Here φL (ν) is the Lorentz profile, which describes the natural line broadening
associated with the transition in question. The non-zero width of this Lorentz profile
implies that resonant scattering is not perfectly coherent; typically the energy of the
outgoing photon will be slightly different from that of the incident photon. The
probability distribution for this energy shift is described by φL (ν), and originates
from the Heisenberg Uncertainty Principle, according to which ∆E ∆t ≥ h̄/2;
hence, the uncertainty related to the time it takes for the electron to spontaneously
de-excite results in a related ‘uncertainty’ in energy.
163
Figure 29: Illustration of the scattering cross section of an atom or ion with at
least one bound electron. At high (low) frequency, scattering is in the Thomson
(Rayleigh) regime; at specific, intermediate frequencies, set by the transition energies
of the atom/ion, resonant scattering dominates; the profiles are Lorentz profiles, and
reflect the natural line broadening. The relative heights of the peaks are set by their
oscillator strengths. NOTE: figure is not to scale; typically the cross section for
resonant scattering is orders of magnitude larger than the Thomson cross section.
164
Figure 30: Example of a quasar spectrum revealing the Ly-α forest due to resonant
scattering of Ly-α photons by neutral hydrogen along the line-of-sight from quasar to
observer.
165
CHAPTER 21
HI + γ → p + e ,
where HI denotes a neutral hydrogen atom. Using that the rate, Γ, of an interaction
is always given by Γ = t−1
coll = nσv, the photoionization rate, Γγ,H , is proportional
to the number density of ionizing photons and to the photoionization cross section,
σpi (ν), according to:
Z ∞
Γγ,H = c σpi (ν) nγ (ν) dν
νt
166
4 π J(ν)
nγ (ν) = .
chν
The photoionization cross sections can be obtained from quantum electrodynamics
by calculating the bound-free transition probability of an atom in a radiation field
(see e.g., Rybicki & Lightman 1979).
p + e → HI + γ .
For hydrogen (or a hydrogenic ion, i.e., an ion with a single electron), the recom-
bination cross section to form an atom (or ion) at level n, σrec (v, n), is related to
the corresponding photoionization cross section by the Milne relation:
2
gn hν
σrec (v, n) = σpi (ν, n) ,
gn+1 me c v
where gn = 2n2 is the statistical weight of energy level n and ν and v are related by
me v 2 /2 = h(ν − νn ), with hνn the threshold energy required to ionize an atom whose
electron sits in energy state n. The recombination coefficient for a given level n is
the product of the capture cross section and velocity, σrec (v, n) v, averaged over the
velocity distribution f (v). For an optically thin gas where all photons produced by
recombination can escape without being absorbed, the total recombination coefficient
is the sum over all n:
∞
X ∞ Z
X
αA = αn = σrec (v, n) v f (v) dv
n=1 n=1
167
Strömgren sphere: A sphere of ionized hydrogen (H II) around an ionizing source
(e.g., AGN, O or B star, etc.). Ionization of hydrogen (from the ground state) requires
a photon energy of at least 13.6eV, which implies UV photons. In a (partially)
ionized medium, electrons and nuclei recombine to produce neutral atoms. The
region around an ionizing source will ultimately establish ionization equilibrium
in which the number of ionizations is equal to the number of recombinations.
Consider an ionizing source in a uniform medium of pure hydrogen. Let Ṅion be the
number of ionizing photons produced per second. The corresponding recombination
rate is given by
4
Ṅrec = ne np αrec V = n2e αB πRs3
3
where we have used that, for a pure hydrogen gas, ne = np , and Rs is the radius of
the Strömgren sphere (i.e., the radius of the sphere that is going to be ionized),
which can be written as !1/3
3 Ṅion
Rs =
4 π αB n2e
Using that the luminosity of the ionizing source, L∗ , is related to its surface intensity,
I∗ , according to
L∗ = 4 π R∗2 F∗ = 4 π 2 R∗2 I∗
where R∗ is the radius of the ionizing source (i.e., an O-star) and we have used that
F∗ = π I∗ (see Chapter 18). Hence, we have that
Z ∞ Z ∞
2 2 Bν (T ) π L∗ Bν (T )
Ṅion = 4 π R∗ dν = 4
dν
νt hν σSB Teff νt hν
where we have assumed that the ionizing source is a Black Body of temperature T ,
and, in the second part, that L∗ = 4πR∗2 σSB Teff
4
.
Thus, by measuring the luminosity and effective temperature of a star, and the radius
of its Strömgren sphere, one can infer the (electron) density of its surroundings.
168
CHAPTER 22
As we have seen, there are numerous processes by which a photon can interact with
matter. It is useful to define the mean-free path, l, for a photon and the related
opacity and optical depth.
dτν = σν n dl = κν ρ dl = αν dl
Here σν is the effective cross section ([σν ] = cm2 ), κν is the mass absorption
coefficient ([κν ] = cm2 g−1 ), αν is the absorption coefficient ([αν ] = cm−1 ), and
n and ρ are the number and mass densities, respectively. The optical depth to a
source at distance d is therefore
Z d Z d
τν = dτν = κν (l) ρ(l) dl
0 0
The ISM/IGM between source and observer is said to be optically thick (thin) if
τν > 1 (τν < 1).
169
Rosseland Mean Opacities: In the case of stars opacity is crucially important
for understanding stellar structure. Opacities within stars are typically expressed
in terms of the Rosseland mean opacities, κ, which is a weighted average of κν
over frequency. Typically, one finds that κ ∝ ρ T −3.5 , which is known as Kramer’s
opacity law, and is a consequence of the fact that the opacity is dominated by
bound-free and/or free-free absorption. A larger opacity implies stronger radiation
pressure, which gives rise to the concept of the Eddington luminosity.
Eddington Luminosity: the maximum luminosity a star (or, more general, emit-
ter) can achieve before the star’s radiation pressure starts to exceed the force of
gravity.
4 π G M∗ c
LEdd =
κν
which is called the Eddington luminosity. Stars with L > LEdd cannot exist, as
they would blow themselves appart (Frad > Fgrav ). The most massive stars known
have luminosities that are very close to their Eddington luminosity.
Since the luminosities of AGN (supermassive black holes with accretion disks) are set
by their accretion rate, the same argument implies an upper limit to the accretion
rate of AGN, known as the Eddington limit.
170
Figure 31: Empirical extinction laws, defined in terms of the ratio of color excesses,
for the Milky Way (MW), the Large Magellanic Cloud (LMC) and the Small Mag-
ellanic Cloud (SMC). Note the strong feature around 2100 Å in the MW extinction
law, believed to be due to graphite dust grains.
Extinction by Dust: Dust grains can scatter and absorb photons. Their ability
to do so depends on (i) grain size, (ii) grain composition, and (iii) the presence of a
magnetic field, which can cause grain alignment. Observationally, the extinction in
the V -band is defined by
fV IV
AV ≡ −2.5 log = −2.5 log
fV,0 IV,0
where the subscript zero refers to the unextincted flux/intensity. Using that IV =
IV,0 e−τV we have that
AV = 1.086 τV
More generally, Aλ = 1.086 τλ ; hence, an optical depth of unity roughly corresponds
to an extinction of one magnitude.
Reddening: In addition to extinction, dust also causes reddening, due to the fact
that dust extinction is more effective at shorter (bluer) wavelengths.
171
Color Excess: E(B − V ) ≡ AB − AV , which can also be defined for any other
wavebands.
Theoretical attempts to model the extinction curve have shown that dust comes in
two varieties, graphites and silicates, while the grain-size distribution is well fit by
dN/da ∝ a−3.5 and covers the range from ∼ 0.005µm to ∼ 0.25µm. Note that for
radiation with λ > amax ≃ 2500 Å dust mainly causes Rayleigh scattering.
172
CHAPTER 23
Radiative Transfer
Consider an incoming signal of specific intensity Iν,0 passing through a cloud (i.e.,
any gaseous region). As the radiation transits a small path length dr through the
cloud, its specific intensity changes by dIν = dIν,loss +dIν,gain . The loss-term describes
the combined effect of scattering and absorption, which remove photons from the line-
of-sight, while the gain-term describes all processes that add photons to the line-of-
sight; these include all emission processes from the gas itself, as well as scattering of
photons from any direction into the line-of-sight.
In what follows we ignore the contribution of scattering to dIν,gain , as this term makes
solving the equation of radiative transfer much more complicated. We will briefly
comments on that below, but for now the only process that is assumed to contribute
to dIν,gain are emission processes from the gas.
• Emission coefficient, jν , defined as the energy emitted per unit time, per unit
volume, per unit frequency, per unit solid angle (i.e., dE = jν dt dV dν dΩ, and thus
[jν ] = erg s−1 cm−3 Hz−1 sr−1 ).
dIν
= −αν Iν + jν (form I)
dr
dIν
= −Iν + Sν (form II)
dτν
173
Here Sν ≡ jν /αν is called the source function, and has units of specific intensity
(i.e., [Sν ] = erg s−1 cm−2 Hz−1 sr−1 ). In order to derive form II from form I, recall
that dτν = αν dr (see Chapter 22).
NOTE: we use the convention of τν increasing from the source towards the observer.
Some textbooks adopt the opposite convention, which results in some sign differences.
Case A No Cloud
In this case, there is no absorption (αν = 0) or emission (jν = 0), other than the
emission from the background source. Hence, we have that
dIν
=0 ⇒ Iν = Iν,0
dr
which expresses that intensity is a conserved quantity in vacuum.
where l is the size of the cloud along the line-of-sight. This equation simply expresses
that the increase of intensity is equal to the emission coefficient integrated along the
line-of-sight.
174
Case D Cloud in Thermodynamic Equilibrium w/o Background Source
Iν = Sν = Bν (T )
jν = αν Bν (T )
The latter of these equivalent relations is sometimes called Kirchoff’s law, and
simply expresses that a black body needs to establish a balance between emission
and absorption (i.e., Bν (T ) = jν /αν ).
Using that I˜ν,0 = Iν,0 e0 = Iν,0 the solution to this simple integral equation is
Z τν
−τν ′
Iν = Iν,0 e + Sν (τν′ ) e−(τν −τν ) dτν′
0
175
where τν is the total optical depth along the line of sight (i.e., through the cloud).
The above is the formal solution, which, under the simplifying assumption that the
source function is constant along the line of sight reduces to
Iν = Iν,0 e−τν + Sν 1 − e−τν
The first term expresses the attenuation of the background signal, the second
term expresses the added signal due to the emission from the cloud, while
the third term describes the cloud’s self-absorption.
Using the above formal solution to the equation of radiative transfer, we have the
following two extremes:
τν ≫ 1 ⇒ I ν = Sν
τν ≪ 1 ⇒ Iν = Iν,0 (1 − τν ) + Sν τν
where, for the latter case, we have used the Taylor series expansion for the exponen-
tial. In the high optical depth case, the observer just ‘sees’ the outer layers of the
cloud, and therefore the observed intensity is simply the source function of the cloud
(the observed signal contains no contribution from the background source). In the
small optical depth limit, the contribution from the cloud is suppressed by a factor
τν , while that from the background source is attenuated by a factor (1 − τν ). Using
that Sν = jν /αν and τν = αν l (if the absorption coefficient is constant throughout the
cloud), we see that τν Sν = jν l; in other words, the contribution from the cloud itself
is simply its emission coefficient (assumed constant throughout the cloud) multiplied
with the pathlength through the cloud.
To get some further insight into the source function and radiative transfer in general,
consider form II of the radiate transfer equation. If Iν > Sν then dIν /dτν < 0, so that
the specific intensity decreases along the line of sight. If, on the other hand, Iν < Sν
then dIν /dτν > 0, indicating that the specific intensity increases along the line of
sight. Hence, Iν tends towards Sν . If the optical depth of the cloud is sufficiently
large than this ‘tendency’ will succeed, and Iν = Sν .
An important special case of the general solution derived above is if the cloud is in
local thermal equilibrium (LTE). This is very often the case, since over the mean
free path of the photons, every system will tend to be in LTE, unless it was recently
176
disturbed and has not yet been able to equilibrate. In the case of LTE, we have that,
over a patch smaller than or equal to the mean free path of the photons, we have that
Sν ≡ jν /αν = Bν (T ), where T is the kinetic temperature (= radiation temperature)
of the patch.
Note that Iν is not constant throughout the cloud, as was the case for a cloud in
TE. In the case of LTE, however, there can be a non-zero gradient dIν /dr.
Keeping this difference in mind, we now look at the solution to our equation of
radiative transfer for a cloud in LTE at its two extremes:
τν ≫ 1 ⇒ Iν = Bν (T )
τν ≪ 1 ⇒ Iν = Iν,0 (1 − τν ) + Bν (T ) τν
The former expresses that an optically thick cloud in LTE emits black body radiation.
This is characterized by the fact that (i) if there is a background source, you can’t see
it, (ii) you can look into the source only for about one mean free path of the photons
(which is much smaller than the size of the source), and (iii) the only information
available to an observer is the temperature of the cloud (the observed intensity is a
Planck curve of temperature T ).
In the optically thin limit, the observed intensity depends on the background source
(if present), and depends on both the temperature (sets source function) and density
(sets optical depth) of the cloud (recall that τν ∝ κν ρ l).
177
In the case of LTE, however, there are radial gradients, which are responsible for
diminishing the intensity by the optical depth in the case where τν ≪ 1. This
may seem somewhat ‘counter-intuitive’, as it indicates that a cloud of larger optical
depth is more intense!!! To understand this, consider the limit τν → 0. In this case
all photons pass through the cloud with zero probability to be absorbed/scattered.
In this situation, there is simply no way to establish an equilibrium between emission
and absorption required for the establishment of a black body; or, put differently,
if there is no absorption, there is no emission either (after all, we are in LTE), and
thus, Iν = 0.
Based on the above, we have that, in the case of a cloud in LTE without background
source, Iν ≤ Bν (T ), where T is the temperature of the cloud. If we express the
intensity in terms of the brightness temperature we have that TB,ν ≤ T . Hence,
for a cloud in LTE without background source the observed brightness temperature
is a lower limit on the kinetic temperature of the cloud.
What about scattering? In the most general case, any element in the cloud re-
ceives radiation coming from all 4π sterradian, and a certain fraction of that radiation
will be scattered into the line-of-sight of an observer.
In general, the scattering can (will) be non-isotropic (e.g., Thomson scattering) and
incoherent (e.g., Compton scattering or resonant scattering), and the final equation
of radiative transfer can only be solved numerically.
jν,scat = αν,scat Jν
178
The source function due to scattering is then simply
Z
jν,scat 1
Sν ≡ = Jν = Iν dΩ
αν,scat 4π
Hence, the source function due to isotropic, coherent scattering is simply the
mean intensity.
The radiative transfer equation for pure scattering (no background source, and no
emission) is
dIν
= −αν,scat (Iν − Jν )
dr
Even this oversimplified case of pure isotropic, coherent scattering is not easily solved.
Since Jν involves an integration (over all 4π sterradian), the above equation is an
integro-differential equation, which are extremely difficult to solve in general; one
typically has to resort to numerical methods (see Rybicki & Lightmann 1979 for
more details). NOTE: although the scattering may be isotropic, the incoming radia-
tion is typically not.
179
Absorption Line: Bν (T ) < Iν,0 T < TB,ν
where T is the (kinetic) temperature of the cloud, and TB,ν is the brightness tem-
perature of the background source, at the frequency of the line. Hence, if the cloud
is colder (hotter) than the source, an absorption (emission) line will arise. In the
case of no background source, we effectively have that TB,ν = 0, and the cloud will
thus reveal an emission line. In the case where T = TB,ν no line will be visible,
independent of the optical depth of the cloud!
180
CHAPTER 24
Continuum radiation is any radiation that forms a continuous spectrum and is not
restricted to a narrow frequency range. In what follows we briefly describe five
continuum emission mechanisms:
• Two-Photon emission
• Synchrotron emission
In general, the way to proceed is to ‘derive’ the emission coeffient, jν , the absorption
coefficient, αν , and then use the equation of radiative transfer to compute the
specific intensity, Iν , (i.e., the ‘spectrum’), for a cloud of gas emitting continuum
radiation using any one of those mechanisms.
First some general remarks: when talking about continuum processes it is important
to distinguish thermal emission, in which the radiation is generated by the thermal
motion of charged particles and in which the intensity therefore depends (at least) on
temperature, i.e., Iν = Iν (T, ..), from non-thermal emission, which is everything
else.
Examples of thermal continuum emission are black body radiation and (thermal)
bremsstrahlung, while synchrotron radiation is an example of non-thermal emission.
Another non-thermal continuum mechanism is inverse compton radiation. However,
since this is basically an incoherent photon-scattering mechanism, rather than a
photon-production mechanism, we will not discuss IC scattering any further here.
181
Characteristics of Thermal Continuum Emission:
• Low Brightness Temperatures: Since one rarely encouters gases with kinetic
temperatures T > 107 − 108 K, and since TB ≤ T (see Chapter 28), if the brightness
temperature of the radiation exceeds ∼ 108 K it is most likely non-thermal in origin
(or has experienced IC scattering).
Thermal Radiation & Black Body Radiation: Thermal radiation is the con-
tinuum emission arising from particles colliding, which causes acceleration of charges
(atoms typically have electric or magnetic dipole moments, and colliding those results
in the emission of photons). This thermal radiation tries to establish thermal equi-
librium with the matter that produces it via photon-matter interactions. If thermal
equilibrium is established (locally), then the source function Sν ≡ jν /αν = Bν (T )
(Kirchoff’s law).
As we have seen in the previous Chapter;
Bν (T ) if τν ≫ 1
Iν =
τν Bν (T ) if τν ≪ 1
where τν = αν l is the optical depth through the cloud, which has a dimension l
along the line-of-sight.
182
electron moving with velocity v when experiencing a Coulomb interaction with a
charge Ze over an impact parameter b (see Rybicki & Lightmann 1979 for a de-
tailed derivation).
The next step is to integrate over all possible impact parameters. This are all impact
parameters b > bmin , where from a classical perspective bmin is set by the requirement
that the kinetic energy of the electron, Ek = 12 me v 2 , is larger than the binding
energy, Eb = Ze2 /b (otherwise we are in the regime of recombination; see below).
However, there are some quantum mechanical corrections one needs to make to
this bmin which arise from Heisenberg’s Uncertainty Principle (∆x ∆p ≥ h̄/2). This
correction factor is called the free-free Gaunt factor, gff (ν, Te ), which is close to
unity, and has only a weak frequency dependence. The final step in obtaining the
emission coeffient is the integration over the Maxwellian velocity distribution of
the electrons, characterized by Te . The result (in erg s−1 cm−3 Hz−1 sr−1 ) is:
−39 Z2
jν = 5.44 × 10 1/2
ne ni gff (ν, Te ) e−hν/kB Te
Te
In the case of a pure (ionized) hydrogen gas, Z = 1 and ni = ne . Upon inspection, it
is clear that free-free emission has a flat spectrum jν ∝ ν α with α ∼ 0 (controlled by
the weak frequency dependence of the Gaunt factor) with an exponential cut-off for
h ν > kB Te (the maximum photon energy is set by the temperature of the electrons).
This reveals that a measurement of the exponential cut-off is a direct measure of the
electron temperature.
The above emission coefficient tells us the emissive behavior of a pocket of (ion-
ized) gas without allowance for the internal absorption. Accounting for the latter
requires radiative transfer. Since Bremsstrahlung arises from collisions, we may use
the LTE approximation. Hence, Kirchoff’s law tells us that αν = jν /Bν (T ), which
allows us to compute the absorption coefficent, and thus the optical depth τν = αν l.
Substitution of Bν (T ), with T = Te , yields
τν ≃ 3.7 × 108 Z 2 Te−1/2 ν −3 [1 − e−hν/kB Te ] gff (ν, Te ) E
where Z
E≡ n2e dl ≃ n2e l
is called the emission measure, and we have assumed that ne = ni . Upon in-
spection, one notices that τν ∝ ν −2 (for hν ≪ kB Te ), indicating that the opacity
183
of the cloud increases with decreasing frequency. This opacity arises from free-free
absorption, which is simply the inverse process of free-free emission; a photon is
absorbed by an electron that is experiencing a Coulomb interaction.
184
Figure 32: Specific intensity of free-free emission (Bremssstrahlung), including the
effect of free-free self absorption at low frequencies, where the optical depth exceeds
unity. At low frequencies, one probes the Rayleigh-Jeans part of the Planck curve
corresponding to the electron temperature. At intermediate frequencies, where the
cloud is optically thin, the spectrum is flat, followed by an exponential cut-off at the
high-frequency end.
185
Two-Photon Emission: two photon emission occurs between bound states in an
atom, but it produces continuum emission rather than line emission.
Two photon emission occurs when an electron finds itself in a quantum level for which
any downward transition would violate quantum mechanical selection rules. Each
transition is therefore highly forbidden. However, there is a non-zero chance that
the electron decays under the emission of two, rather than one, photons. Energy
conservation guarantees that ν1 + ν2 = νtr = ∆Etr /h, where ∆Etr is the energy
difference associated with the transition. The most probable configuration is the
one in which ν1 = ν2 = νtr /2, but all configurations that satisfy the above energy
conservation are possible; they become less likely the larger |ν1 − νtr /2|, resulting in
a ‘continuum’ emission that appears as an extremely broad ‘emission line’. In fact,
whereas the number of photons with 0 < ν < νtr/2 is equal to that with νtr/2 < ν < νtr ,
the latter have more energy (i.e., Eγ = hν). Consequently, the spectral energy
distribution, Lν (erg s−1 Hz−1 ) is skewed towards higher frequency.
For two photon emission to occur, we require that spontaneous emission happens
before collisional de-excitation has a chance. Consequently, two-photon emission oc-
curs in low density ionized gas. The strength of the two photon emission depends on
the number of particles in the excited states. This in turn depends on the recombi-
nation rate; although two-photon emission is quantum-mechanical in nature, it can
still be throught of as ‘thermal emission’, and the density dependence is the same as
for free-free and free-bound emission (i.e., jν ∝ n2e ).
186
Figure 33: Emission spectra of plasmas with solar abundances. The histogram indi-
cates the total spectrum, including line radiation. The spectrum has been binned in
order to highlight the relative importance of line radiation. The thick solid line is the
total continuum emission, the thin solid line the contribution due to Bremsstrahlung,
the dashed line free-bound emission and the dotted line two-photon emission. Note
how recombination becomes less and less important when the gas gets hotter. [From
Kaastra et al. 2008, Space Science Reviews, 134, 155]
187
experiences a Lorentz force:
~
v ~ = e v B sin φ = e v B⊥ = e v⊥ B
F~e = e ×B
c c c c
where φ is the pitch angle between ~v and B.~ If φ = 0 the particle moves along the
magnetic field, and the Lorentz force is zero. If φ = 90o the particle will move in a
circle around the magnetic field line, while for 0o < φ < 90o the electron will spiral
(‘cork-screw’) around the magnetic field line. In the latter two cases, the electron is
being accelerated, which causes the emission of photons. Note that this applies to
both electrons and ions. However, since the cyclotron (synchrotron) emission from
ions is negligble compared to that from electrons, we will focus on the latter.
If the particle is non-relativistic, then the emission is called cyclotron emission. If,
on the other hand, the particles are relativistic, the emission is called synchrotron
emission. We will first focus on the former.
Cyclotron emission: the gyrating electron emits dipolar emission that (i) has the
frequency of gyration, and (ii) is highly polarized. Depending on the viewing angle
~ linear
the observer can see circular polarization (if line-of-sight is alined with B),
~ or elliptical polarization (for any
polarization, if line of sight is perpendicular to B,
other orientation).
The gyration frequency can be obtained by equating the Lorentz force with the
centripetal force:
2
e v⊥ me v⊥
Fe = B=
c r0
where v⊥ = v sin φ, which results in
me v⊥ c
r0 =
eB
~ This is called the gyration radius (or gyro-radius). The period of
where B = |B|.
gyration is T = 2πr0 /v⊥ , which implies a gyration frequency (i.e., the frequency
of the emitted photons) of
1 eB
ν0 = =
T 2π me c
188
Note that this frequency is independent of the velocity of the electron! It only depends
on the magnetic field strength B;
ν0 ~
|B|
= 2.8
MHz Gauss
We thus see that cyclotron emission really is line emission, rather than continuum
emission. The nature of this line emission is very different though, from ‘normal’
spectral lines which result from quantum transitions within atoms or molecules.
Note, though, that if the ‘source’ has a varying magnetic field, then the variance in
B will result in a ‘broadening’ of the line, which, if sufficiently large, may appear as
‘continuum emission’.
this reason, cyclotron emission is rarely observed, with the exception of the Sun,
some of the planets in our Solar System, and an occasional pulsar.
Synchrotron Emission: this is the same as cyclotron emission, but in the limit in
which the electrons are relativistic. As we demonstrate below, this has two important
effects: it makes the gyration frequency dependent on the energy (velocity) of the
electron, and it causes strong beaming of the electron’s dipole emission.
189
Figure 34: Illustration of how the Lorentz transformation from the electron rest frame
to the lab frame introduce relativistic beaming with an opening angle θ = 1/γ. Note
that in the electron rest frame, the synchrotron emission is dipole emission.
Note that now the gyration frequency does depend on the velocity (energy) of the
(relativistic) electrons, which in principle implies that because the electrons will
have a distribution in energies, the synchrotron emission is going to be continuum
emission. However, you can also see that the gyration frequency is even lower than
in the case of cyclotron emission, by a factor 1/γ. For the record, Lorentz factors of
up to ∼ 1011 have been measured, indicating that γ can be extremely large! Hence,
if the photon emission were to be at the gyration frequency, we would never be able
to see it, because of the plasme-frequency-shielding.
However, the gyration frequency is not the only frequency in this problem. Because
of the relativistic motion, the dipole emission from the electron, as seen from the
observer’s frame, is highly beamed (see Fig. 34), with an opening angle ∼ 1/γ (which
can thus be tiny). Consequently, the observer does not have a continuous view of the
electron, but only sees EM radiation when the beam sweeps over the line-of-sight.
The width of these ‘pulses’ are a factor 1/γ 3 shorter than the gyration period. The
corresponding frequency, called the critical frequency, is given by
3e
νcrit = γ 2 B⊥
4 π me c
which translates into
νcrit B⊥
= 4.2 γ 2
MHz Gauss
190
Figure 35: Specific intensity of synchrotron emission, including the effect of syn-
chrotron self absorption at low frequencies, where the optical depth exceeds unity.
So although the gyration frequency will be small, the critical frequency can be ex-
tremely large. This critical frequency corresponds to the shortest time period (the
pulse duration), and therefore represents the largest frequency, above which the emis-
sion is negligble. The longest time period, which is related to the gyration period,
determines the fundamental frequency
νf 2.8 ~
|B|
=
MHz γ sin2 φ Gauss
The emission spectrum due to synchrotron radiation will contain this fundamen-
tal frequency plus all its harmonics up to νcrit . Since these harmonics are ex-
tremely closely spaced (after all, the gyration frequency is extremely small), the
synchrotron spectrum for one value of γ looks essentially continuum. When taking
the γ-distribution into account (which is related to the energy distribution of the
relativistic electrons), the distribution becomes trully continuum, and the critical
and fundamental frequencies no longer can be discerned (because they depend on γ).
After integrating over the energy distribution of the relativistic electrons, which typ-
ically has a power-law distribution N(E) ∝ E −Γ one obtains the following emission
191
and absorption coefficients:
(Γ+1)/2
jν ∝ B⊥ ν −(Γ−1)/2
(Γ+2)/2
αν ∝ B⊥ ν −(Γ+4)/2
Using that typically Γ > 0, we have that τν ∝ ν a with a < 0; synchrotron self-
absorption becomes more important at lower frequencies.
where α ≡ − Γ−1 2
. Fig. 35 shown an illustration of a typical synchrotron spec-
trum: at low frequencies, where τν ≫ 1, we have that Iν ∝ ν 5/2 , which transits
to Iν ∝ ν −(Γ−1)/2 once the emitting medium becomes optically thin for synchroton
self-absorption. Note that there is no cut-off related to the critical frequency, since
νcrit = νcrit (E).
192
Supplemental Material
Appendices
Appendices A-E present background material on calculus relevant for this course.
The other Appendices present supplemental material that is NOT considered part of
this course’s curriculum. They are included to provide background information for
those readers that want to know a bit more.
193
Appendix A
Vector Calculus
~ = (a1 , a2 , a3 ) = a1 î + a2 ĵ + a3 k̂
Vector: A
p
~ =
Amplitude of vector: |A| a21 + a22 + a23
~ =1
Unit vector: |A|
Basis: In the above example, the unit vectors î, ĵ and k̂ form a vector basis.
~ B
Any 3 vectors A, ~ and C
~ can form a vector basis
~ B,
as long as det(A, ~ C)
~ 6= 0.
~ B)
~ a1 a2
Determinant: det(A, = = a1 b2 − a2 b1
b1 b2
a1 a2 a3
~ B,
~ C)
~ b2 b3 b3 b1 b1 b2
det(A, = b1 b2 b3 = a1 + a2 + a3
c2 c3 c3 c1 c1 c2
c1 c2 c3
Geometrically: ~ B)
det(A, ~ = ± area of parallelogram
~ B,
det(A, ~ C)
~ = ± volume of parallelepiped
~+B
Summation of vectors: A ~ =B
~ +A
~ = (a1 + b1 , a2 + b2 , a3 + b3 )
194
P
Einstein Summation Convention: ai bi = i ai bi = a1 b1 + a2 b2 + a3 b3 = ~a · ~b
∂Ai /∂xi = ∂A1 /∂x1 + ∂A2 /∂x2 + ∂A3 /∂x3 = ∇ · A ~
Aii = A11 + A22 + A33 = Tr A ~ (trace of A)~
~ ·B
• check orthogonality: two vectors are orthogonal if A ~ =0
~ in direction of A,
• compute projection of B ~ which is given by A
~ · B/|
~ A|~
î ĵ k̂
~×B
Cross Product (aka vector product): A ~ = a1 a2 a3 = εijk ai bj êk
b1 b2 b3
~ ~ ~ |B|
|A × B| = |A| ~ sin θ = det(A,
~ B)
~
NOTE: εijk is called the Levi-Civita tensor, which is described in Appendix E,.
In addition to the dot product and cross product, there is a third vector product
that one occasionally encounters in dynamics;
The tensor product AB is a tensor of rank two and is called a dyad. The sum of two
~ B,
or more dyads is called a dyadic. For example, let A, ~ C~ and D~ be four vectors,
from which we can form the dyads AB and CD. Their sum AB + CD is then a
dyadic. Note that in general a dyadic it not a dyad because it cannot be written
as a vector multiplied with a vector. Hence, dyadics differ from vectors in that the
sum of two vectors is a vector whereas the sum of two dyads is not necessarily a
dyad.
195
~·B
A ~ =B
~ ·A
~ ~×B
A ~ = −B
~ ×A
~
~ ·B
(αA) ~ = α(A
~ · B)
~ =A
~ · (αB)
~ ~ ×B
(αA) ~ = α(A
~ × B)
~ =A
~ × (αB)
~
~ · (B
A ~ + C)
~ =A
~·B
~ +A
~·C
~ ~ × (B
A ~ + C)
~ =A
~ ×B
~ +A
~×C
~
~·B
A ~ =0 → ~⊥B
A ~ ~ ×B
A ~ =0 → ~kB
A ~
~·A
A ~ = |A|
~2 ~ ×A
A ~=0
~ · (B
Triple Scalar Product: A ~ × C)
~ = det(A,
~ B,
~ C)
~ = εijk ai bj ck
~ · (B
A ~ × C)
~ = 0 → A, ~ B,
~ C
~ are coplanar
~ · (B
A ~ × C)
~ =B~ · (C
~ × A)
~ =C~ · (A~ × B)~
~ × (B
Triple Vector Product: A ~ × C)
~ = (A~ · C)
~ B~ − (A
~ · B)
~ C~
as is clear from above, A~ × (B
~ × C)
~ lies in plane of B
~ and C.
~
~ × B)
~ · (C
~ × D)
~ = (A
~ ~ ~ ~ ~ ~ ~ ~
Useful to remember: (A h · C) (B · D)i− (A ·hD) (B · C) i
~ × B)
(A ~ × (C~ × D)
~ = A ~ · (B
~ × D)
~ C ~− A ~ · (B
~ × C)
~ D ~
Gradient Operator: ∇ = ∇ ~ = ∂, ∂, ∂
∂x ∂y ∂z
This vector operator is sometimes called the nabla or del operator.
∂ 2 ∂ 2∂ 2
Laplacian operator: ∇2 = ∇ · ∇ = ∂x 2 + ∂y 2 + ∂z 2
∂f ∂f ∂f
Differential: f = f (x, y, z) → df = ∂x
dx + ∂y
dy + ∂z
dz
196
∂f ∂f ∂f
Gradient Vector: ∇f = gradf = , ,
∂x ∂y ∂z
the gradient vector at (x, y, z) is normal to the level surface
through the point (x, y, z).
î ĵ k̂
Curl of Vector Field: curlF~ = ∇ × F~ = ∂/∂x ∂/∂y ∂/∂z
Fx Fy Fz
~
A vector field for which ∇ × F = 0 is called irrotational or curl-free.
~ x) and B(~
Let S(~x) and T (~x) be scalar fields, and let A(~ ~ x) be vector fields:
~ = divA
∇·A ~ = scalar ~ = (∇ · ∇) A
∇2 A ~ = vector
~ = curlA
∇×A ~ = vector
197
∇ × (∇S) = 0 curl grad S = 0
~ =0
∇ · (∇ × A) ~=0
div curl A
∇(ST ) = S ∇T + T ∇S
~ = S(∇ · A)
∇ · (S A) ~ +A
~ · ∇S
~ = (∇S) × A
∇ × (S A) ~ + S(∇ × A)
~
~ × B)
∇ · (A ~ =B
~ · (∇ × A)
~ −A
~ · (∇ × B)
~
~ × B)
∇ × (A ~ = A(∇
~ · B)
~ − B(∇
~ ~ + (B
· A) ~ · ∇)A
~ − (A
~ · ∇)B
~
~ · B)
∇(A ~ = (A
~ · ∇)B
~ + (B
~ · ∇)A
~ +A
~ × (∇ × B)
~ +B
~ × (∇ × A)
~
~ × (∇ × A)
A ~ = 1 ∇(A
~ · A)
~ − (A
~ · ∇)A
~
2
~ = ∇2 (∇ × A)
∇ × (∇2 A) ~
198
Appendix B
• F~ (~x) is a gradient field, which means that there is a scalar field Φ(~x) so that
F~ = ∇Φ
H
• Path independence: c F~ · d~l = 0
• Irrotational = curl-free: ∇ × F~ = 0
199
Appendix C
Integral Theorems
NOTE: in the first equation we have used that ∇ × F~ is always pointing in the
direction of the normal n̂.
NOTE: The curve of the line intergral must have positive orientation, meaning that
d~l points counterclockwise when the normal of the surface points towards the viewer.
200
Appendix D
where ~e1 , ~e2 and ~e3 are the unit directional vectors in the new (q1 , q2 , q3 )-coordinate
system. In what follows we show how to properly treat such generalized coordinate
systems.
In general, one expresses the distance between (q1 , q2 , q3 ) and (q1 + dq1 , q2 + dq2 , q3 +
dq3 ) in an arbitrary coordinate system as
p
ds = hij dqi dqj
Here hij is called the metric tensor. In what follows, we will only consider orthog-
onal coordinate systems for which√hij = 0 if i 6= j, so that ds2 = h2i dqi2 (Einstein
summation convention) with hi = hii .
An example of an orthogonal coordinate system are the Cartesian coordinates, for
which hij = δij . After all, the distance between two points separated by the infinites-
imal displacement vector d~x = (dx, dy, dz) is ds2 = |d~x|2 = dx2 + dy 2 + dz 2 .
201
The coordinates (x, y, z) and (q1 , q2 , q3 ) are related to each other via the transfor-
mation relations
x = x(q1 , q2 , q3 )
y = y(q1 , q2 , q3 )
z = z(q1 , q2 , q3 )
∂~x
hi =
∂qi
202
Using this expression for the metric allows us to write the unit directional vectors
as
1 ∂~x
~ei =
hi ∂qi
and the differential vector in the compact form as
From the latter we also have that the infinitesimal volume element for a general
coordinate system is given by
Note that the absolute values are needed to assure that d3~x is positive.
203
~ In the Cartesian basis C = {~ex , ~ey , ~ez } we have that
Now consider a vector A.
~ C = Ax ~ex + Ay ~ey + Az ~ez
[A]
In the basis B = {~e1 , ~e2 , ~e3 }, corresponding to our generalized coordinate system,
we instead have that
~ B = A1 ~e1 + A2 ~e2 + A3 ~e3
[A]
We can rewrite the above as
e11 e21 e31 A1 e11 + A2 e21 + A3 e31
~ B = A1 e12 + A2 e22 + A3 e32 = A2 e12 + A2 e22 + A3 e32
[A]
e13 e23 e33 A3 e13 + A2 e23 + A3 e33
and thus
e11 e21 e31 A1 A1
~ B = e12 e22 e32 A2 ≡ T A2
[A]
e13 e23 e33 A3 A3
Using similar logic, one can write
ex1 ey1 ez1 Ax 1 0 0 Ax Ax
~ C = ex2 ey2 ez2 Ay = 0 1 0 Ay = I Ay
[A]
ex3 ey3 ez3 Az 0 0 1 Az Az
~ is the same object independent of its basis we have that
and since A
Ax A1
I Ay = T A2
Az A3
~ B and [A]
and thus, we see that the relation between [A] ~ C is given by
~ C = T [A]
[A] ~ B, ~ B = T−1 [A]
[A] ~C
For this reason, T is called the transformation of basis matrix. Note that the
columns of T are the unit-direction vectors ~ei , i.e., Tij = eij . Since these are or-
thogonal to each other, the matric T is said to be orthogonal, which implies that
T−1 = T T (the inverse is equal to the transpose), and det(T ) = ±1.
Now we are finally ready to determine how to write our position vector ~x in the new
basis B of our generalized coordinate system. Let’s write ~x = ai ~ei , i.e.
a1
[~x]B = a2
a3
204
We started this appendix by pointing out that it is tempting, but
p wrong, to set ai = qi
(as for the Cartesian basis). To see this, recall that |~x| = (a1 )2 + (a2 )2 + (a3 )2 ,
from which it is immediately clear that each ai needs to have the dimension of length.
Hence, when qi is an angle, clearly ai 6= qi . To compute the actual ai you need to
use the transformation of basis matrix as follows:
e11 e12 e13 x e11 x + e12 y + e13 z
[~x]B = T−1 [~x]C = e21 e22 e23 y = e21 x + e22 y + e23 z
e31 e32 e33 z e31 x + e32 y + e33 z
Hence, using our expression for the unit direction vectors, we see that
1 ∂xj 1 ∂~x
ai = xj = · ~x
hi ∂qi hi ∂qi
and by operating d/dt on [~x]B we find that the corresponding velocity vector in the
B basis is given by X
[~v ]B = hi q̇i ~ei
i
with q̇i = dqi /dt. Note that the latter can also be inferred more directly by simply
dividing the expression for the differential vector (d~x = hi qi ~ei ) by dt.
205
Next we write out the gradient, the divergence, the curl and the Laplacian for our
generalized coordinate system:
The gradient:
1 ∂ψ
∇ψ = ~ei
hi ∂qi
The divergence:
~= 1 ∂ ∂ ∂
∇·A (h2 h3 A1 ) + (h3 h1 A2 ) + (h1 h2 A3 )
h1 h2 h3 ∂q1 ∂q2 ∂q3
The Laplacian:
2 1 ∂ h2 h3 ∂ψ ∂ h3 h1 ∂ψ ∂ h1 h2 ∂ψ
∇ ψ= + +
h1 h2 h3 ∂q1 h1 ∂q1 ∂q2 h2 ∂q2 ∂q3 h3 ∂q3
206
Vector Calculus in Cylindrical Coordinates:
hR = 1 hφ = R hz = 1
AR = Ax cos φ − Ay sin φ
Aφ = −Ax sin φ + Ay cos φ
Az = Az
The Gradient:
~ = 1 ∂ (RAR ) + 1 ∂Aφ + ∂Az
∇·A
R ∂R R ∂φ ∂z
207
The Laplacian:
2 1 ∂ ∂ψ 1 ∂2ψ ∂2ψ
scalar : ∇ψ = R + + 2
R ∂R ∂R R2 ∂φ2 ∂z
2~ 2 FR 2 ∂Fθ
vector : ∇ F = ∇ FR − 2 − 2 ~eR
R R ∂θ
2 2 ∂FR Fθ
+ ∇ Fθ + 2 − 2 ~eθ
R ∂θ R
2
+ ∇ Fz ~ez
208
Vector Calculus in Spherical Coordinates:
hr = 1 hθ = r hφ = r sin θ
~v = ṙ ~er + r ~e˙ r
= ṙ ~er + r θ̇ ~eθ + r sin θ φ̇ ~eφ
The Gradient:
~ = 1 ∂ (r 2 Ar ) + 1
∇·A
∂
(sin θAθ ) +
1 ∂Aφ
2
r ∂r r sin θ ∂θ r sin θ ∂φ
The Convective Operator:
~ ~ ∂Br Aθ ∂Br Aφ ∂Br Aθ Bθ + Aφ Bφ
(A · ∇) B = Ar + + − ~er
∂r r ∂θ r sin θ ∂φ r
∂Bθ Aθ ∂Bθ Aφ ∂Bθ Aθ Br Aφ Bφ cotθ
+ Ar + + + − ~eθ
∂r r ∂θ r sin θ ∂φ r r
∂Bφ Aθ ∂Bφ Aφ ∂Bφ Aφ Br Aφ Bθ cotθ
+ Ar + + + + ~eφ
∂r r ∂θ r sin θ ∂φ r r
209
The Laplacian:
2 1 ∂ 2 ∂ψ 1 ∂ ∂φ 1 ∂2ψ
scalar : ∇ψ = 2 r + 2 sin θ + 2 2
r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂ψ 2
2~ 2 2Fr 2 ∂(Fθ sin θ) 2 ∂Fφ
vector : ∇F = ∇ Fr − 2 − 2 − 2 ~er
r r sin θ ∂θ r sin θ ∂φ
2 2 ∂Fr Fθ 2 cos θ ∂Fφ
+ ∇ Fθ + 2 − 2 − ~eθ
r ∂θ r sin θ r 2 sin2 θ ∂φ
2 2 ∂Fr 2 cos θ ∂Fθ Fφ
+ ∇ Fφ + 2 + 2 2 − 2 2 ~eφ
r sin θ ∂φ r sin θ ∂φ r sin θ
210
Appendix E
The Levi-Civita symbol, also known as the permutation symbol or the anti-
symmetric symbol,is a collection of numbers, defined from the sign of a permu-
tation of the natural numbers 1, 2, 3, ..., n. It is often encountered in linear algebra,
vector and tensor calculus, and differential geometry.
The n-dimensional Levi-Civita symbol is indicated by εi1 i2 ...in , where each index
i1 , i2 , ..., in takes values 1, 2, ..., n, and has the defining property that the symbol is
total antisymmetric in all its indices: when any two indices are interchanged, the
symbol is negated:
ε...ip...iq ... = −ε...iq ...ip ...
If any two indices are equal, the symbol is zero, and when all indices are unequal,
we have that
εi1 i2 ...in = (−1)p ε1,2,...n
where p is called the parity of the permutation. It is the number of pairwise inter-
changes necessary to unscramble i1 , i2 , ..., in into the order 1, 2, ..., n. A permutation
is said to be even (odd) if its parity is an even (odd) number.
211
Appendix F
1 ∂ui ∂uj
eij = +
2 ∂xj ∂xi
1 ∂ui ∂uj
ξij = −
2 ∂xj ∂xi
The symmetric part of the deformation tensor, eij , is called the rate of strain
tensor, while the anti-symmetric part, ξij , expresses the vorticity w ~ ≡ ∇ × ~u in
1
the velocity field, i.e., ξij = − 2 εijk wk . Note that one can always find a coordinate
system for which eij is diagonal. The axes of that coordinate frame indicate the
eigendirections of the strain (compression or stretching) on the fluid element.
In terms of the relation between the viscous stress tensor, τij , and the deformation
tensor, Tkl , there are a number of properties that are important.
212
• Locality: the τij − Tkl -relation is said to be local if the stress tensor is only
a function of the deformation tensor and thermodynamic state functions like
temperature.
• Linearity: the τij − Tkl -relation is said to be linear if the relation between
the stress and rate-of-strain is linear. This is equivalent to saying that τij does
not depend on ∇2~u or higher-order derivatives.
Note that (in a Newtonian fluid) the viscous stress tensor depends only on the sym-
metric component of the deformation tensor (the rate-of-strain tensor eij ), but not
on the antisymmetric component which describes vorticity. You can understand
the fact that viscosity and vorticity are unrelated by considering a fluid disk in solid
body rotation (i.e., ∇ · ~u = 0 and ∇ × ~u = w ~ 6= 0). In such a fluid there is no
”slippage”, hence no shear, and therefore no manifestation of viscosity.
213
Thus far we have derived that the stress tensor, σij , which in principle has 6 un-
knowns, can be reduced to a function of three unknowns only (P , µ, λ) as long as
the fluid is Newtonian. Note that these three scalars, in general, are functions of
temperature and density. We now focus on these three scalars in more detail, starting
with the pressure P . To be exact, P is the thermodynamic equilibrium pres-
sure, and is normally computed thermodynamically from some equation of state,
P = P (ρ, T ). It is related to the translational kinetic energy of the particles when
the fluid, in equilibrium, has reached equipartition of energy among all its degrees
of freedom, including (in the case of molecules) rotational and vibrations degrees of
freedom.
Pm = P − η ∇ · ~u
where
2 P − Pm
η = µ+λ=
3 ∇ · ~u
is the coefficient of bulk viscosity. We can now write the stress tensor as
∂ui ∂uj 2 ∂uk ∂uk
σij = −P δij + µ + − δij + η δij
∂xj ∂xi 3 ∂xk ∂xk
This is the full expression for the stress tensor in terms of the coefficients of shear
viscosity, µ, and bulk viscosity, η.
214
Appendix G
Consider a system which can exchange energy and particles with a reservoir, and
the volume of which can change. There are three ways for this system to increase
its internal energy; heating, changing the system’s volume (i.e., doing work on the
system), or adding particles. Hence,
dU = T dS − P dV + µ dN
Note that this is the first law of thermodynamics, but now with the added possibility
of changing the number of particles of the system. The scalar quantity µ is called
the chemical potential, and is defined by
∂U
µ=
∂N S,V
This is not to be confused with the µ used to denote the mean weight per particle,
which ALWAYS appears in combination with the proton mass, mp . As is evident
from the above expression, the chemical potential quantifies how the internal energy
of the system changes if particles are added or removed, while keeping the entropy
and volume of the system fixed. The chemical potential appears in the Fermi-Dirac
distribution describing the momentum distribution of a gas of fermions or bosons.
Consider an ideal gas, of volume V , entropy S and with internal energy U. Now
imagine adding a particle of zero energy (ǫ = 0), while keeping the volume fixed.
Since ǫ = 0, we also have that dU = 0. But what about the entropy? Well, we
have increased the number of ways in which we can redistribute the energy U (a
macrostate quantity) over the different particles (different microstates). Hence, by
adding this particle we have increased the system’s entropy. If we want to add a
particle while keeping S fixed, we need to decrease U to offset the increase in the
number of ‘degrees of freedom’ over which to distribute this energy. Hence, keeping
S (and V ) fixed, requires that the particle has negative energy, and we thus see that
µ < 0.
215
For a fully degenerate Fermi gas, we have that T = 0, and thus S = 0 (i.e., there is
only one micro-state associated with this macrostate, and that is the fully degenerate
one). If we now add a particle, and demand that we keep S = 0, then that particle
must have the Fermi energy (see Chapter 6); ǫ = Ef . Hence, for a fully degenerate
gas, µ = Ef .
To end this discussion of the chemical potential, we address the origin of its name,
which may, at first, seem weird. Let’s start with the ‘potential’ part. The origin of
this name is clear from the following. According to its definition (see above), the
chemical potential is the ‘internal energy’ per unit amount (moles). Now consider
the following correspondences:
These examples make it clear why µ is considered a ‘potential’. Finally, the word
chemical arises from the fact that the µ plays an important role in chemistry (i.e.,
when considering systems in which chemical reactions take place, which change the
particles). In this respect, it is important to be aware of the fact that µ is an
additive quantity that is conserved in a chemical reaction. Hence, for a chemical
216
reaction i + j → k + l one has that µi + µj = µk + µl . As an example, consider the
annihilation of an electron and a positron into two photons. Using that µ = 0 for
photons, we see that the chemical potential of elementary particles (i.e., electrons)
must be opposite to that of their anti-particles (i.e., positrons).
Because of the additive nature of the chemical potential, we also have that the
above equation for dU changes slightly whenever the gas consists of different particle
species; it becomes X
dU = T dS − P dV + µi dNi
i
where the summation is over all species i. If the gas consists of equal numbers
of elementary particles and anti-particles, then the total chemical potential of the
system will be P equal to zero. In fact, in many treatments of fluid dynamics it may be
assumed that i µi dNi = 0; in particular when the relevant reactions are ‘frozen’
(i.e., occur on a timescales τreact that are much longer than the dynamical timescales
τdyn of interest), so that dNi = 0, or if the reactions go so fast (τreact ≪ τdyn )
that each
P reaction and its inverse are in local thermodynamic equilibrium, in which
case i µi dNi = 0 for those species involved in the reaction. Only in the rare,
intermediate case when τreact ∼ τdyn is it important to keep track of the relative
abundances of the various chemical and/or nuclear species.
217