Astrophysical fluid dynamics

Abstract

These lecture notes and example problems are based on a course given at the University of Cambridge in Part III of the Mathematical Tripos. Fluid dynamics is involved in a very wide range of astrophysical phenomena, such as the formation and internal dynamics of stars and giant planets, the workings of jets and accretion discs around stars and black holes and the dynamics of the expanding Universe. Effects that can be important in astrophysical fluids include compressibility, self-gravitation and the dynamical influence of the magnetic field that is ‘frozen in’ to a highly conducting plasma. The basic models introduced and applied in this course are Newtonian gas dynamics and magnetohydrodynamics (MHD) for an ideal compressible fluid. The mathematical structure of the governing equations and the associated conservation laws are explored in some detail because of their importance for both analytical and numerical methods of solution, as well as for physical interpretation. Linear and nonlinear waves, including shocks and other discontinuities, are discussed. The spherical blast wave resulting from a supernova, and involving a strong shock, is a classic problem that can be solved analytically. Steady solutions with spherical or axial symmetry reveal the physics of winds and jets from stars and discs. The linearized equations determine the oscillation modes of astrophysical bodies, as well as their stability and their response to tidal forcing.

1 Introduction

1.1 Areas of application

Astrophysical fluid dynamics (AFD) is a theory relevant to the description of the interiors of stars and planets, exterior phenomena such as discs, winds and jets and also the interstellar medium, the intergalactic medium and cosmology itself. A fluid description is not applicable (i) in regions that are solidified, such as the rocky or icy cores of giant planets (under certain conditions) and the crusts of neutron stars, and (ii) in very tenuous regions where the medium is not sufficiently collisional (see § 2.9.3).

Important areas of application include:

1. (i) Instabilities in astrophysical fluids

2. (ii) Convection

3. (iii) Differential rotation and meridional flows in stars

4. (iv) Stellar oscillations driven by convection, instabilities or tidal forcing

5. (v) Astrophysical dynamos

6. (vi) Magnetospheres of stars, planets and black holes

7. (vii) Interacting binary stars and Roche-lobe overflow

8. (viii) Tidal disruption and stellar collisions

9. (ix) Supernovae

10. (x) Planetary nebulae

11. (xi) Jets and winds from stars and discs

12. (xii) Star formation and the physics of the interstellar medium

13. (xiii) Astrophysical discs, including protoplanetary discs, accretion discs in interacting binary stars and galactic nuclei, planetary rings, etc.

14. (xiv) Other accretion flows (Bondi, Bondi–Hoyle, etc.)

15. (xv) Processes related to planet formation and planet–disc interactions

16. (xvi) Planetary atmospheric dynamics

17. (xvii) Galaxy clusters and the physics of the intergalactic medium

18. (xviii) Cosmology and structure formation

1.2 Theoretical varieties

There are various flavours of AFD in common use. The basic model involves a compressible, inviscid fluid and is Newtonian (i.e. non-relativistic). This is known as hydrodynamics (HD) or gas dynamics (to distinguish it from incompressible hydrodynamics). The thermal physics of the fluid may be treated in different ways, either by assuming it to be isothermal or adiabatic, or by including radiative processes in varying levels of detail.

Magnetohydrodynamics (MHD) generalizes this theory by including the dynamical effects of a magnetic field. Often the fluid is assumed to be perfectly electrically conducting (ideal MHD). One can also include the dynamical (rather than thermal) effects of radiation, resulting in a theory of radiation (magneto)hydrodynamics. Dissipative effects such as viscosity and resistivity can be included. All these theories can also be formulated in a relativistic framework.

1. (i) HD: hydrodynamics

2. (ii) MHD: magnetohydrodynamics

3. (iii) RHD: radiation hydrodynamics

4. (iv) RMHD: radiation magnetohydrodynamics

5. (v) GRHD: general relativistic hydrodynamics

6. (vi) GRRMHD: general relativistic radiation magnetohydrodynamics, etc.

1.3 Characteristic features

AFD typically differs from ‘laboratory’ or ‘engineering’ fluid dynamics in the relative importance of certain effects. Compressibility and gravitation are often important in AFD, while magnetic fields, radiation forces and relativistic phenomena are important in some applications. Effects that are often unimportant in AFD include viscosity, surface tension and the presence of solid boundaries.

2 Ideal gas dynamics

2.1 Fluid variables

A fluid is characterized by a velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ and two independent thermodynamic properties. Most useful are the dynamical variables: the pressure $p(\boldsymbol{x},t)$ and the mass density ${\it\rho}(\boldsymbol{x},t)$ . Other properties, e.g. temperature $T$ , can be regarded as functions of $p$ and ${\it\rho}$ . The specific volume (volume per unit mass) is $v=1/{\it\rho}$ .

We neglect the possible complications of variable chemical composition associated with chemical and nuclear reactions, ionization and recombination.

2.2 Eulerian and Lagrangian viewpoints

In the Eulerian viewpoint we consider how fluid properties vary in time at a point that is fixed in space, i.e. attached to the (usually inertial) coordinate system. The Eulerian time-derivative is simply the partial differential operator

(2.1) $$\begin{eqnarray}\frac{\partial }{\partial t}.\end{eqnarray}$$

In the Lagrangian viewpoint we consider how fluid properties vary in time at a point that moves with the fluid at velocity $\boldsymbol{u}(\boldsymbol{x},t)$ . The Lagrangian time derivative is then

(2.2) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}=\frac{\partial }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}.\end{eqnarray}$$

2.3 Material points and structures

A material point is an idealized fluid element, a point that moves with the bulk velocity $\boldsymbol{u}(\boldsymbol{x},t)$ of the fluid. (Note that the true particles of which the fluid is composed have in addition a random thermal motion.) Material curves, surfaces and volumes are geometrical structures composed of fluid elements; they move with the fluid flow and are distorted by it.

An infinitesimal material line element ${\it\delta}\boldsymbol{x}$ (figure 1) evolves according to

(2.3) $$\begin{eqnarray}\frac{\text{D}{\it\delta}\boldsymbol{x}}{\text{D}t}={\it\delta}\boldsymbol{u}={\it\delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}.\end{eqnarray}$$

It changes its length and/or orientation in the presence of a velocity gradient. (Since ${\it\delta}\boldsymbol{x}$ is only a time-dependent vector rather than a vector field, the time derivative could be written as an ordinary derivative $\text{d}/\text{d}t$ . The notation $\text{D}/\text{D}t$ is used here to remind us that ${\it\delta}\boldsymbol{x}$ is a material structure that moves with the fluid.)

Figure 1. Examples of material line, surface and volume elements.

Infinitesimal material surface and volume elements can be defined from two or three material line elements according to the vector product and the triple scalar product (figure 1)

(2.4a,b ) $$\begin{eqnarray}{\it\delta}\boldsymbol{S}={\it\delta}\boldsymbol{x}^{(1)}\times {\it\delta}\boldsymbol{x}^{(2)},\quad {\it\delta}V={\it\delta}\boldsymbol{x}^{(1)}\boldsymbol{\cdot }{\it\delta}\boldsymbol{x}^{(2)}\times {\it\delta}\boldsymbol{x}^{(3)}.\end{eqnarray}$$

They therefore evolve according to

(2.5a,b ) $$\begin{eqnarray}\frac{\text{D}{\it\delta}\boldsymbol{S}}{\text{D}t}=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})\,{\it\delta}\boldsymbol{S}-(\boldsymbol{{\rm\nabla}}\boldsymbol{u})\boldsymbol{\cdot }{\it\delta}\boldsymbol{S},\quad \frac{\text{D}{\it\delta}V}{\text{D}t}=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})\,{\it\delta}V,\end{eqnarray}$$

as follows from the above equations (exercise). The second result is easier to understand: the volume element increases when the flow is divergent. These equations are most easily derived using Cartesian tensor notation. In this notation the equation for ${\it\delta}\boldsymbol{S}$ reads

(2.6) $$\begin{eqnarray}\frac{\text{D}{\it\delta}S_{i}}{\text{D}t}=\frac{\partial u_{j}}{\partial x_{j}}\,{\it\delta}S_{i}-\frac{\partial u_{j}}{\partial x_{i}}\,{\it\delta}S_{j}.\end{eqnarray}$$

2.4 Equation of mass conservation

The equation of mass conservation,

(2.7) $$\begin{eqnarray}\frac{\partial {\it\rho}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u})=0,\end{eqnarray}$$

has the typical form of a conservation law: ${\it\rho}$ is the mass density (mass per unit volume) and ${\it\rho}\boldsymbol{u}$ is the mass flux density (mass flux per unit area). An alternative form of the same equation is

(2.8) $$\begin{eqnarray}\frac{\text{D}{\it\rho}}{\text{D}t}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}.\end{eqnarray}$$

If ${\it\delta}m={\it\rho}{\it\delta}V$ is a material mass element, it can be seen that mass is conserved in the form

(2.9) $$\begin{eqnarray}\frac{\text{D}{\it\delta}m}{\text{D}t}=0.\end{eqnarray}$$

2.5 Equation of motion

The equation of motion,

(2.10) $$\begin{eqnarray}{\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}-\boldsymbol{{\rm\nabla}}p,\end{eqnarray}$$

derives from Newton’s second law per unit volume with gravitational and pressure forces. ${\it\Phi}(\boldsymbol{x},t)$ is the gravitational potential and $\boldsymbol{g}=-\boldsymbol{{\rm\nabla}}{\it\Phi}$ is the gravitational field. The force due to pressure acting on a volume $V$ with bounding surface $S$ is

(2.11) $$\begin{eqnarray}-\int _{S}p\,\text{d}\boldsymbol{S}=\int _{V}(-\boldsymbol{{\rm\nabla}}p)\,\text{d}V.\end{eqnarray}$$

Viscous forces are neglected in ideal gas dynamics.

2.6 Poisson’s equation

The gravitational potential is related to the mass density by Poisson’s equation,

(2.12) $$\begin{eqnarray}{\rm\nabla}^{2}{\it\Phi}=4{\rm\pi}G{\it\rho},\end{eqnarray}$$

where $G$ is Newton’s constant. The solution

(2.13) $$\begin{eqnarray}{\it\Phi}(\boldsymbol{x},t)={\it\Phi}_{\mathit{int}}+{\it\Phi}_{\mathit{ext}}=-G\int _{V}\frac{{\it\rho}(\boldsymbol{x}^{\prime },t)}{|\boldsymbol{x}^{\prime }-\boldsymbol{x}|}\,\text{d}^{3}\boldsymbol{x}^{\prime }-G\int _{\hat{V}}\frac{{\it\rho}(\boldsymbol{x}^{\prime },t)}{|\boldsymbol{x}^{\prime }-\boldsymbol{x}|}\,\text{d}^{3}\boldsymbol{x}^{\prime }\end{eqnarray}$$

generally involves contributions from both the fluid region $V$ under consideration and the exterior region $\hat{V}$ .

A non-self-gravitating fluid is one of negligible mass for which ${\it\Phi}_{\mathit{int}}$ can be neglected. More generally, the Cowling approximation ¹ consists of treating ${\it\Phi}$ as being specified in advance, so that Poisson’s equation is not coupled to the other equations.

2.7 Thermal energy equation

In the absence of non-adiabatic heating (e.g. by viscous dissipation or nuclear reactions) and cooling (e.g. by radiation or conduction),

(2.14) $$\begin{eqnarray}\frac{\text{D}s}{\text{D}t}=0,\end{eqnarray}$$

where $s$ is the specific entropy (entropy per unit mass). Fluid elements undergo reversible thermodynamic changes and preserve their entropy.

This condition is violated in shocks (see § 6.3).

The thermal variables $(T,s)$ can be related to the dynamical variables $(p,{\it\rho})$ via an equation of state and standard thermodynamic identities. The most important case is that of an ideal gas together with black-body radiation,

(2.15) $$\begin{eqnarray}p=p_{g}+p_{r}=\frac{k{\it\rho}T}{{\it\mu}m_{H}}+\frac{4{\it\sigma}T^{4}}{3c},\end{eqnarray}$$

where $k$ is Boltzmann’s constant, $m_{H}$ is the mass of the hydrogen atom, ${\it\sigma}$ is Stefan’s constant and $c$ is the speed of light. ${\it\mu}$ is the mean molecular weight (the average mass of the particles in units of $m_{H}$ ), equal to $2.0$ for molecular hydrogen, $1.0$ for atomic hydrogen, $0.5$ for fully ionized hydrogen and approximately $0.6$ for ionized matter of typical cosmic abundances. Radiation pressure is usually negligible except in the centres of high-mass stars and in the immediate environments of neutron stars and black holes. The pressure of an ideal gas is often written in the form ${\mathcal{R}}{\it\rho}T/{\it\mu}$ , where ${\mathcal{R}}=k/m_{H}$ is a version of the universal gas constant.

We define the first adiabatic exponent

(2.16) $$\begin{eqnarray}{\it\Gamma}_{1}=\left(\frac{\partial \ln p}{\partial \ln {\it\rho}}\right)_{s},\end{eqnarray}$$

which is related to the ratio of specific heat capacities

(2.17) $$\begin{eqnarray}{\it\gamma}=\frac{c_{p}}{c_{v}}=\frac{T\displaystyle \left(\frac{\partial s}{\partial T}\right)_{p}}{T\displaystyle \left(\frac{\partial s}{\partial T}\right)_{v}}\end{eqnarray}$$

by (exercise)

(2.18) $$\begin{eqnarray}{\it\Gamma}_{1}={\it\chi}_{{\it\rho}}{\it\gamma},\end{eqnarray}$$

where

(2.19) $$\begin{eqnarray}{\it\chi}_{{\it\rho}}=\left(\frac{\partial \ln p}{\partial \ln {\it\rho}}\right)_{T}\end{eqnarray}$$

can be found from the equation of state. We can then rewrite the thermal energy equation as

(2.20) $$\begin{eqnarray}\frac{\text{D}p}{\text{D}t}=\frac{{\it\Gamma}_{1}p}{{\it\rho}}\frac{\text{D}{\it\rho}}{\text{D}t}=-{\it\Gamma}_{1}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}.\end{eqnarray}$$

For an ideal gas with negligible radiation pressure, ${\it\chi}_{{\it\rho}}=1$ and so ${\it\Gamma}_{1}={\it\gamma}$ . Adopting this very common assumption, we write

(2.21) $$\begin{eqnarray}\frac{\text{D}p}{\text{D}t}=-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}.\end{eqnarray}$$

2.8 Simplified models

A perfect gas may be defined as an ideal gas with constant $c_{v}$ , $c_{p}$ , ${\it\gamma}$ and ${\it\mu}$ . Equipartition of energy for a classical gas with $n$ degrees of freedom per particle gives ${\it\gamma}=1+2/n$ . For a classical monatomic gas with $n=3$ translational degrees of freedom, ${\it\gamma}=5/3$ . This is relevant for fully ionized matter. For a classical diatomic gas with two additional rotational degrees of freedom, $n=5$ and ${\it\gamma}=7/5$ . This is relevant for molecular hydrogen. In reality ${\it\Gamma}_{1}$ is variable when the gas undergoes ionization or when the gas and radiation pressures are comparable. The specific internal energy (or thermal energy) of a perfect gas is

(2.22) $$\begin{eqnarray}e=\frac{p}{({\it\gamma}-1){\it\rho}}\left[=\frac{n}{{\it\mu}m_{H}}\frac{1}{2}kT\right].\end{eqnarray}$$

(Note that each particle has an internal energy of $kT/2$ per degree of freedom, and the number of particles per unit mass is $1/{\it\mu}m_{H}$ .)

A barotropic fluid is an idealized situation in which the relation $p({\it\rho})$ is known in advance. We can then dispense with the thermal energy equation. e.g. if the gas is strictly isothermal and perfect, then $p=c_{s}^{2}{\it\rho}$ with $c_{s}=\text{const.}$ being the isothermal sound speed. Alternatively, if the gas is strictly homentropic and perfect, then $p=K{\it\rho}^{{\it\gamma}}$ with $K=\text{const.}$

An incompressible fluid is an idealized situation in which $\text{D}{\it\rho}/\text{D}t=0$ , implying $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$ . This can be achieved formally by taking the limit ${\it\gamma}\rightarrow \infty$ . The approximation of incompressibility eliminates acoustic phenomena from the dynamics.

The ideal gas law itself is not valid at very high densities or where quantum degeneracy is important.

2.9 Microphysical basis

It is useful to understand the way in which the fluid dynamical equations are derived from microphysical considerations. The simplest model involves identical neutral particles of mass $m$ and negligible size with no internal degrees of freedom.

2.9.1 The Boltzmann equation

Between collisions, particles follow Hamiltonian trajectories in their six-dimensional $(\boldsymbol{x},\boldsymbol{v})$ phase space:

(2.23a,b ) $$\begin{eqnarray}{\dot{x}}_{i}=v_{i},\quad \dot{v}_{i}=a_{i}=-\frac{\partial {\it\Phi}}{\partial x_{i}}.\end{eqnarray}$$

The distribution function $f(\boldsymbol{x},\boldsymbol{v},t)$ specifies the number density of particles in phase space. The velocity moments of $f$ define the number density $n(\boldsymbol{x},t)$ in real space, the bulk velocity $\boldsymbol{u}(\boldsymbol{x},t)$ and the velocity dispersion $c(\boldsymbol{x},t)$ according to

(2.24a-c ) $$\begin{eqnarray}\int f\,\text{d}^{3}\boldsymbol{v}=n,\quad \int \boldsymbol{v}f\,\text{d}^{3}\boldsymbol{v}=n\boldsymbol{u},\quad \int |\boldsymbol{v}-\boldsymbol{u}|^{2}f\,\text{d}^{3}\boldsymbol{v}=3nc^{2}.\end{eqnarray}$$

Equivalently,

(2.25) $$\begin{eqnarray}\int v^{2}f\,\text{d}^{3}\boldsymbol{v}=n(u^{2}+3c^{2}).\end{eqnarray}$$

The relation between velocity dispersion and temperature is $kT=mc^{2}$ .

In the absence of collisions, $f$ is conserved following the Hamiltonian flow in phase space. This is because particles are conserved and the flow in phase space is incompressible (Liouville’s theorem). More generally, $f$ evolves according to Boltzmann’s equation,

(2.26) $$\begin{eqnarray}\frac{\partial f}{\partial t}+v_{j}\frac{\partial f}{\partial x_{j}}+a_{j}\frac{\partial f}{\partial v_{j}}=\left(\frac{\partial f}{\partial t}\right)_{c}.\end{eqnarray}$$

The collision term on the right-hand side is a complicated integral operator but has three simple properties corresponding to the conservation of mass, momentum and energy in collisions:

(2.27a-c ) $$\begin{eqnarray}\int m\left(\frac{\partial f}{\partial t}\right)_{c}\,\text{d}^{3}\boldsymbol{v}=0,\quad \int m\boldsymbol{v}\left(\frac{\partial f}{\partial t}\right)_{c}\,\text{d}^{3}\boldsymbol{v}=\mathbf{0},\quad \int \frac{1}{2}mv^{2}\left(\frac{\partial f}{\partial t}\right)_{c}\,\text{d}^{3}\boldsymbol{v}=0.\end{eqnarray}$$

The collision term is local in $\boldsymbol{x}$ (not even involving derivatives) although it does involve integrals over $\boldsymbol{v}$ . The equation $(\partial f/\partial t)_{c}=0$ has the general solution

(2.28) $$\begin{eqnarray}f=f_{M}=(2{\rm\pi}c^{2})^{-3/2}n\,\exp \left(-\frac{|\boldsymbol{v}-\boldsymbol{u}|^{2}}{2c^{2}}\right),\end{eqnarray}$$

with parameters $n$ , $\boldsymbol{u}$ and $c$ that may depend on $\boldsymbol{x}$ . This is the Maxwellian distribution.

2.9.2 Derivation of fluid equations

A crude but illuminating model of the collision operator is the Bhatnagar–Gross–Krook (BGK) approximation

(2.29) $$\begin{eqnarray}\left(\frac{\partial f}{\partial t}\right)_{c}\approx -\frac{1}{{\it\tau}}(f-f_{M}),\end{eqnarray}$$

where $f_{M}$ is a Maxwellian distribution with the same $n$ , $\boldsymbol{u}$ and $c$ as $f$ and ${\it\tau}$ is the relaxation time. This can be identified approximately with the mean free flight time of particles between collisions. In other words the collisions attempt to restore a Maxwellian distribution on a characteristic time scale ${\it\tau}$ . They do this by randomizing the particle velocities in a way consistent with the conservation of momentum and energy.

If the characteristic time scale of the fluid flow is much greater than ${\it\tau}$ , then the collision term dominates the Boltzmann equation and $f$ must be very close to $f_{M}$ . This is the hydrodynamic limit.

The velocity moments of $f_{M}$ can be determined from standard Gaussian integrals, in particular (exercise)

(2.30a,b ) $$\begin{eqnarray}\int f_{M}\,\text{d}^{3}\boldsymbol{v}=n,\quad \int v_{i}\,f_{M}\,\text{d}^{3}\boldsymbol{v}=nu_{i},\end{eqnarray}$$

(2.31a,b ) $$\begin{eqnarray}\int v_{i}v_{j}f_{M}\,\text{d}^{3}\boldsymbol{v}=n(u_{i}u_{j}+c^{2}{\it\delta}_{ij}),\quad \int v^{2}v_{i}\,f_{M}\,\text{d}^{3}\boldsymbol{v}=n(u^{2}+5c^{2})u_{i}.\end{eqnarray}$$

We obtain equations for mass, momentum and energy by taking moments of the Boltzmann equation weighted by $(m,mv_{i},mv^{2}/2)$ . In each case the collision term integrates to zero because of its conservative properties, and the $\partial /\partial v_{j}$ term can be integrated by parts. We replace $f$ with $f_{M}$ when evaluating the left-hand sides and note that $mn={\it\rho}$ :

(2.32) $$\begin{eqnarray}\frac{\partial {\it\rho}}{\partial t}+\frac{\partial }{\partial x_{i}}({\it\rho}u_{i})=0,\end{eqnarray}$$

(2.33) $$\begin{eqnarray}\frac{\partial }{\partial t}({\it\rho}u_{i})+\frac{\partial }{\partial x_{j}}\left[{\it\rho}(u_{i}u_{j}+c^{2}{\it\delta}_{ij})\right]-{\it\rho}a_{i}=0,\end{eqnarray}$$

(2.34) $$\begin{eqnarray}\frac{\partial }{\partial t}\left(\frac{1}{2}{\it\rho}u^{2}+\frac{3}{2}{\it\rho}c^{2}\right)+\frac{\partial }{\partial x_{i}}\left[\left(\frac{1}{2}{\it\rho}u^{2}+\frac{5}{2}{\it\rho}c^{2}\right)u_{i}\right]-{\it\rho}u_{i}a_{i}=0.\end{eqnarray}$$

These are equivalent to the equations of ideal gas dynamics in conservative form (see § 4) for a monatomic ideal gas ( ${\it\gamma}=5/3$ ). The specific internal energy is $e=(3/2)c^{2}=(3/2)kT/m$ .

This approach can be generalized to deal with molecules with internal degrees of freedom and also to plasmas or partially ionized gases where there are various species of particle with different charges and masses. The equations of MHD can be derived by including the electromagnetic forces in Boltzmann’s equation.

2.9.3 Validity of a fluid approach

The essential idea here is that deviations from the Maxwellian distribution are small when collisions are frequent compared to the characteristic time scale of the flow. In higher-order approximations these deviations can be estimated, leading to the equations of dissipative gas dynamics including transport effects (viscosity and heat conduction).

The fluid approach breaks down if the mean flight time ${\it\tau}$ is not much less than the characteristic time scale of the flow, or if the mean free path ${\it\lambda}\approx c{\it\tau}$ between collisions is not much less than the characteristic length scale of the flow. ${\it\lambda}$ can be very long (measured in astronomical units or parsecs) in very tenuous gases such as the interstellar medium, but may still be smaller than the size of the system.

Some typical order-of-magnitude estimates:

Solar-type star: centre ${\it\rho}\sim 10^{2}~\text{g}~\text{cm}^{-3}$ , $T\sim 10^{7}~\text{K}$ ; photosphere ${\it\rho}\sim 10^{-7}~\text{g}~\text{cm}^{-3}$ , $T\sim 10^{4}~\text{K}$ ; corona ${\it\rho}\sim 10^{-15}~\text{g}~\text{cm}^{-3}$ , $T\sim 10^{6}~\text{K}$ .

Interstellar medium: molecular clouds $n\sim 10^{3}~\text{cm}^{-3}$ , $T\sim 10~\text{K}$ ; cold medium (neutral) $n\sim 10-100~\text{cm}^{-3}$ , $T\sim 10^{2}~K$ ; warm medium (neutral/ionized) $n\sim 0.1-1~\text{cm}^{-3}$ , $T\sim 10^{4}~K$ ; hot medium (ionized) $n\sim 10^{-3}-10^{-2}~\text{cm}^{-3}$ , $T\sim 10^{6}~K$ .

The Coulomb cross-section for ‘collisions’ (i.e. large-angle scatterings) between charged particles (electrons or ions) is ${\it\sigma}\approx 1\times 10^{-4}(T/\text{K})^{-2}~\text{cm}^{2}$ . The mean free path is ${\it\lambda}=1/(n{\it\sigma})$ .

Related examples (see appendix A): A.1–A.4.

3 Ideal magnetohydrodynamics

3.1 Elementary derivation of the MHD equations

Magnetohydrodynamics (MHD) is the dynamics of an electrically conducting fluid (a fully or partially ionized gas or a liquid metal) containing a magnetic field. It is a fusion of fluid dynamics and electromagnetism.

3.1.1 Galilean electromagnetism

The equations of Newtonian gas dynamics are invariant under the Galilean transformation to a frame of reference moving with uniform velocity $\boldsymbol{v}$ ,

(3.1a,b ) $$\begin{eqnarray}\boldsymbol{x}^{\prime }=\boldsymbol{x}-\boldsymbol{v}t,\quad t^{\prime }=t.\end{eqnarray}$$

Under this change of frame, the fluid velocity transforms according to

(3.2) $$\begin{eqnarray}\boldsymbol{u}^{\prime }=\boldsymbol{u}-\boldsymbol{v},\end{eqnarray}$$

while scalar variables such as $p$ , ${\it\rho}$ and ${\it\Phi}$ are invariant. The Lagrangian time derivative $\text{D}/\text{D}t$ is also invariant, because the partial derivatives transform according to

(3.3a,b ) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}^{\prime }=\boldsymbol{{\rm\nabla}},\quad \frac{\partial }{\partial t^{\prime }}=\frac{\partial }{\partial t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}.\end{eqnarray}$$

In Maxwell’s electromagnetic theory the electric and magnetic fields $\boldsymbol{E}$ and $\boldsymbol{B}$ are governed by the equations

(3.4a-d ) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}}{\partial t}=-\boldsymbol{{\rm\nabla}}\times \boldsymbol{E},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0,\quad \boldsymbol{{\rm\nabla}}\times \boldsymbol{B}={\it\mu}_{0}\left(\boldsymbol{J}+{\it\epsilon}_{0}\frac{\partial \boldsymbol{E}}{\partial t}\right),\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{E}=\frac{{\it\rho}_{e}}{{\it\epsilon}_{0}},\end{eqnarray}$$

where ${\it\mu}_{0}$ and ${\it\epsilon}_{0}$ are the vacuum permeability and permittivity, $\boldsymbol{J}$ is the electric current density and ${\it\rho}_{e}$ is the electric charge density. (In these notes we use rationalized (e.g. SI) units for electromagnetism. In astrophysics it is also common to use Gaussian units, which are discussed in appendix B.)

It is well known that Maxwell’s equations are invariant under the Lorentz transformation of special relativity, with $c=({\it\mu}_{0}{\it\epsilon}_{0})^{-1/2}$ being the speed of light. These equations cannot be consistently coupled with those of Newtonian gas dynamics, which are invariant under the Galilean transformation. To derive a consistent Newtonian theory of MHD, valid for situations in which the fluid motions are slow compared to the speed of light, we must use Maxwell’s equations without the displacement current ${\it\epsilon}_{0}~\partial \boldsymbol{E}/\partial t$ ,

(3.5a-c ) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}}{\partial t}=-\boldsymbol{{\rm\nabla}}\times \boldsymbol{E},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0,\quad \boldsymbol{{\rm\nabla}}\times \boldsymbol{B}={\it\mu}_{0}\boldsymbol{J}.\end{eqnarray}$$

(We will not require the fourth Maxwell equation, involving $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{E}$ , because the charge density will be found to be unimportant.) It is easily verified (exercise) that these pre-Maxwell equations ² are indeed invariant under the Galilean transformation, provided that the fields transform according to

(3.6a-c ) $$\begin{eqnarray}\boldsymbol{E}^{\prime }=\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B},\quad \boldsymbol{B}^{\prime }=\boldsymbol{B},\quad \boldsymbol{J}^{\prime }=\boldsymbol{J}.\end{eqnarray}$$

These relations correspond to the limit of the Lorentz transformation for electromagnetic fields ³ when $|\boldsymbol{v}|\ll c$ and $|\boldsymbol{E}|\ll c|\boldsymbol{B}|$ .

Under the pre-Maxwell theory, the equation of charge conservation takes the simplified form $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{J}=0$ ; this is analogous to the use of $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$ as the equation of mass conservation in the incompressible (highly subsonic) limit of gas dynamics. The equation of energy conservation takes the simplified form

(3.7) $$\begin{eqnarray}\frac{\partial }{\partial t}\left(\frac{B^{2}}{2{\it\mu}_{0}}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{\boldsymbol{E}\times \boldsymbol{B}}{{\it\mu}_{0}}\right)=0,\end{eqnarray}$$

in which the energy density, $B^{2}/2{\it\mu}_{0}$ , is purely magnetic (because $|\boldsymbol{E}|\ll c|\boldsymbol{B}|$ ), while the energy flux density has the usual form of the Poynting vector $\boldsymbol{E}\times \boldsymbol{B}/{\it\mu}_{0}$ . We will verify the self-consistency of the approximations made in Newtonian MHD in § 3.1.4.

3.1.2 Induction equation

In the ideal MHD approximation we regard the fluid as a perfect electrical conductor. The electric field in the rest frame of the fluid therefore vanishes, implying that

(3.8) $$\begin{eqnarray}\boldsymbol{E}=-\boldsymbol{u}\times \boldsymbol{B}\end{eqnarray}$$

in a frame in which the fluid velocity is $\boldsymbol{u}(\boldsymbol{x},t)$ . This condition can be regarded as the limit of a constitutive relation such as Ohm’s law, in which the effects of resistivity (i.e. finite conductivity) are neglected.

From Maxwell’s equations, we then obtain the ideal induction equation,

(3.9) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B}).\end{eqnarray}$$

This is an evolutionary equation for $\boldsymbol{B}$ alone, $\boldsymbol{E}$ and $\boldsymbol{J}$ having been eliminated. The divergence of the induction equation,

(3.10) $$\begin{eqnarray}\frac{\partial }{\partial t}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B})=0,\end{eqnarray}$$

ensures that the solenoidal character of $\boldsymbol{B}$ is preserved.

3.1.3 The Lorentz force

A fluid carrying a current density $\boldsymbol{J}$ in a magnetic field $\boldsymbol{B}$ experiences a bulk Lorentz force

(3.11) $$\begin{eqnarray}\boldsymbol{F}_{m}=\boldsymbol{J}\times \boldsymbol{B}=\frac{1}{{\it\mu}_{0}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}\end{eqnarray}$$

per unit volume. This can be understood as the sum of the Lorentz forces on individual particles of charge $q$ and velocity $\boldsymbol{v}$ ,

(3.12) $$\begin{eqnarray}\sum q\boldsymbol{v}\times \boldsymbol{B}=\left(\sum q\boldsymbol{v}\right)\times \boldsymbol{B}.\end{eqnarray}$$

(The electrostatic force can be shown to be negligible in the limit relevant to Newtonian MHD; see § 3.1.4.)

In Cartesian coordinates

(3.13) $$\begin{eqnarray}\displaystyle ({\it\mu}_{0}\boldsymbol{F}_{m})_{i} & = & \displaystyle {\it\epsilon}_{ijk}\left({\it\epsilon}_{jlm}\frac{\partial B_{m}}{\partial x_{l}}\right)B_{k}\nonumber\\ \displaystyle & = & \displaystyle \left(\frac{\partial B_{i}}{\partial x_{k}}-\frac{\partial B_{k}}{\partial x_{i}}\right)B_{k}\nonumber\\ \displaystyle & = & \displaystyle B_{k}\frac{\partial B_{i}}{\partial x_{k}}-\frac{\partial }{\partial x_{i}}\left(\frac{B^{2}}{2}\right).\end{eqnarray}$$

Thus

(3.14) $$\begin{eqnarray}\boldsymbol{F}_{m}=\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}-\boldsymbol{{\rm\nabla}}\left(\frac{B^{2}}{2{\it\mu}_{0}}\right).\end{eqnarray}$$

The first term can be interpreted as a curvature force due to a magnetic tension $T_{m}=B^{2}/{\it\mu}_{0}$ per unit area in the field lines; if the field is of constant magnitude then this term is equal to $T_{m}$ times the curvature of the field lines, and is directed towards the centre of curvature. The second term is the gradient of an isotropic magnetic pressure

(3.15) $$\begin{eqnarray}p_{m}=\frac{B^{2}}{2{\it\mu}_{0}},\end{eqnarray}$$

which is also equal to the energy density of the magnetic field.

The magnetic tension gives rise to Alfvén waves ⁴ (see later), which travel parallel to the magnetic field with characteristic speed

(3.16) $$\begin{eqnarray}v_{a}=\left(\frac{T_{m}}{{\it\rho}}\right)^{1/2}=\frac{B}{({\it\mu}_{0}{\it\rho})^{1/2}},\end{eqnarray}$$

the Alfvén speed. This is often considered as a vector Alfvén velocity,

(3.17) $$\begin{eqnarray}\boldsymbol{v}_{a}=\frac{\boldsymbol{B}}{({\it\mu}_{0}{\it\rho})^{1/2}}.\end{eqnarray}$$

The magnetic pressure also affects the propagation of sound waves, which become magnetoacoustic waves (or magnetosonic waves; see later).

The combination

(3.18) $$\begin{eqnarray}{\it\Pi}=p+\frac{B^{2}}{2{\it\mu}_{0}}\end{eqnarray}$$

is often referred to as the total pressure, while the ratio

(3.19) $$\begin{eqnarray}{\it\beta}=\frac{p}{B^{2}/2{\it\mu}_{0}}\end{eqnarray}$$

is known as the plasma beta.

3.1.4 Self-consistency of approximations

Three effects neglected in a Newtonian theory of MHD are (i) the displacement current in Maxwell’s equations (compared to the electric current), (ii) the bulk electrostatic force on the fluid (compared to the magnetic Lorentz force) and (iii) the electrostatic energy (compared to the magnetic energy). We can verify the self-consistency of these approximations by using order-of-magnitude estimates or scaling relations. If the fluid flow has a characteristic length scale $L$ , time scale $T$ , velocity $U\sim L/T$ and magnetic field $B$ , then the electric field can be estimated from (3.8) as $E\sim UB$ . The electric current density and charge density can be estimated from Maxwell’s equations as $J\sim {\it\mu}_{0}^{-1}B/L$ and ${\it\rho}_{e}\sim {\it\epsilon}_{0}E/L$ . Hence the ratios of the three neglected effects to the terms that are retained in Newtonian MHD can be estimated as follows:

(3.20) $$\begin{eqnarray}\frac{{\it\epsilon}_{0}|\partial \boldsymbol{E}/\partial t|}{|\boldsymbol{J}|}\sim \frac{{\it\epsilon}_{0}UB/T}{{\it\mu}_{0}^{-1}B/L}\sim \frac{U^{2}}{c^{2}},\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\frac{|{\it\rho}_{e}\boldsymbol{E}|}{|\boldsymbol{J}\times \boldsymbol{B}|}\sim \frac{{\it\epsilon}_{0}E^{2}/L}{{\it\mu}_{0}^{-1}B^{2}/L}\sim \frac{U^{2}}{c^{2}},\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\frac{{\it\epsilon}_{0}|\boldsymbol{E}|^{2}/2}{|\boldsymbol{B}|^{2}/2{\it\mu}_{0}}\sim \frac{U^{2}}{c^{2}}.\end{eqnarray}$$

Therefore Newtonian MHD corresponds to a consistent approximation of relativistic MHD for highly subluminal flows that is correct to the leading order in the small parameter $U^{2}/c^{2}$ .

3.1.5 Summary of the MHD equations

The full set of ideal MHD equations is

(3.23) $$\begin{eqnarray}\frac{\partial {\it\rho}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u})=0,\end{eqnarray}$$

(3.24) $$\begin{eqnarray}{\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}-\boldsymbol{{\rm\nabla}}p+\frac{1}{{\it\mu}_{0}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B},\end{eqnarray}$$

(3.25) $$\begin{eqnarray}\frac{\text{D}s}{\text{D}t}=0,\end{eqnarray}$$

(3.26) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B}),\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0,\end{eqnarray}$$

together with the equation of state, Poisson’s equation, etc., as required. Most of these equations can be written in at least one other way that may be useful in different circumstances.

These equations display the essential nonlinearity of MHD. When the velocity field is prescribed, an artifice known as the kinematic approximation, the induction equation is a relatively straightforward linear evolutionary equation for the magnetic field. However, a sufficiently strong magnetic field will modify the velocity field through its dynamical effect, the Lorentz force. This nonlinear coupling leads to a rich variety of behaviour. Of course, the purely hydrodynamic nonlinearity of the $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ term, which is responsible for much of the complexity of fluid dynamics, is still present.

3.2 Physical interpretation of MHD

There are two aspects to MHD: the advection of $\boldsymbol{B}$ by $\boldsymbol{u}$ (induction equation) and the dynamical back-reaction of $\boldsymbol{B}$ on $\boldsymbol{u}$ (Lorentz force).

3.2.1 Kinematics of the magnetic field

The ideal induction equation

(3.28) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B})\end{eqnarray}$$

has a beautiful geometrical interpretation: the magnetic field lines are ‘frozen in’ to the fluid and can be identified with material curves. This is sometimes known as Alfvén’s theorem.

One way to show this result is to use the identity

(3.29) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B})=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}-\boldsymbol{B}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}+\boldsymbol{u}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B})\end{eqnarray}$$

to write the induction equation in the form

(3.30) $$\begin{eqnarray}\frac{\text{D}\boldsymbol{B}}{\text{D}t}=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}-\boldsymbol{B}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}),\end{eqnarray}$$

and use the equation of mass conservation,

(3.31) $$\begin{eqnarray}\frac{\text{D}{\it\rho}}{\text{D}t}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u},\end{eqnarray}$$

to obtain

(3.32) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}\left(\frac{\boldsymbol{B}}{{\it\rho}}\right)=\left(\frac{\boldsymbol{B}}{{\it\rho}}\right)\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}.\end{eqnarray}$$

This is exactly the same equation satisfied by a material line element ${\it\delta}\boldsymbol{x}$ (2.3). Therefore a magnetic field line (an integral curve of $\boldsymbol{B}/{\it\rho}$ ) is advected and distorted by the fluid in the same way as a material curve.

A complementary property is that the magnetic flux ${\it\delta}{\it\Phi}=\boldsymbol{B}\boldsymbol{\cdot }{\it\delta}\boldsymbol{S}$ through a material surface element is conserved:

(3.33) $$\begin{eqnarray}\displaystyle \frac{\text{D}{\it\delta}{\it\Phi}}{\text{D}t} & = & \displaystyle \frac{\text{D}\boldsymbol{B}}{\text{D}t}\boldsymbol{\cdot }{\it\delta}\boldsymbol{S}+\boldsymbol{B}\boldsymbol{\cdot }\frac{\text{D}{\it\delta}\boldsymbol{S}}{\text{D}t}\nonumber\\ \displaystyle & = & \displaystyle \left(B_{j}\frac{\partial u_{i}}{\partial x_{j}}-B_{i}\frac{\partial u_{j}}{\partial x_{j}}\right){\it\delta}S_{i}+B_{i}\left(\frac{\partial u_{j}}{\partial x_{j}}\,{\it\delta}S_{i}-\frac{\partial u_{j}}{\partial x_{i}}\,{\it\delta}S_{j}\right)\nonumber\\ \displaystyle & = & \displaystyle 0.\end{eqnarray}$$

By extension, we have conservation of the magnetic flux passing through any material surface.

Precisely the same equation as the ideal induction equation,

(3.34) $$\begin{eqnarray}\frac{\partial {\bf\omega}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times {\bf\omega}),\end{eqnarray}$$

is satisfied by the vorticity ${\bf\omega}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}$ in homentropic or barotropic ideal fluid dynamics in the absence of a magnetic field, in which case the vortex lines are ‘frozen in’ to the fluid (see Example A.2). The conserved quantity that is analogous to the magnetic flux through a material surface is the flux of vorticity through that surface, which, by Stokes’s theorem, is equivalent to the circulation $\oint \boldsymbol{u}\boldsymbol{\cdot }\text{d}\boldsymbol{x}$ around the bounding curve. However, the fact that ${\bf\omega}$ and $\boldsymbol{u}$ are directly related by the curl operation, whereas in MHD $\boldsymbol{B}$ and $\boldsymbol{u}$ are indirectly related through the equation of motion and the Lorentz force, means that the analogy between vorticity dynamics and MHD is limited in scope.

Related examples: A.5, A.6.

3.2.2 The Lorentz force

The Lorentz force per unit volume,

(3.35) $$\begin{eqnarray}\boldsymbol{F}_{m}=\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}-\boldsymbol{{\rm\nabla}}\left(\frac{B^{2}}{2{\it\mu}_{0}}\right),\end{eqnarray}$$

can also be written as the divergence of the Maxwell stress tensor:

(3.36) $$\begin{eqnarray}\boldsymbol{F}_{m}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D648},\quad \unicode[STIX]{x1D648}=\frac{1}{{\it\mu}_{0}}\left(\boldsymbol{B}\boldsymbol{B}-\frac{B^{2}}{2}\unicode[STIX]{x1D644}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity tensor. (The electric part of the electromagnetic stress tensor is negligible in the limit relevant for Newtonian MHD, for the same reason that the electrostatic energy is negligible.) In Cartesian coordinates

(3.37a,b ) $$\begin{eqnarray}(\boldsymbol{F}_{m})_{i}=\frac{\partial M_{ji}}{\partial x_{j}},\quad \unicode[STIX]{x1D614}_{ij}=\frac{1}{{\it\mu}_{0}}\left(B_{i}B_{j}-\frac{B^{2}}{2}{\it\delta}_{ij}\right).\end{eqnarray}$$

If the magnetic field is locally aligned with the $x$ -axis, then

(3.38) $$\begin{eqnarray}\unicode[STIX]{x1D648}=\left[\begin{array}{@{}ccc@{}}T_{m} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0\end{array}\right]-\left[\begin{array}{@{}ccc@{}}p_{m} & 0 & 0\\ 0 & p_{m} & 0\\ 0 & 0 & p_{m}\end{array}\right],\end{eqnarray}$$

showing the magnetic tension and pressure.

Combining the ideas of magnetic tension and a frozen-in field leads to the picture of field lines as elastic strings embedded in the fluid. Indeed there is a close analogy between MHD and the dynamics of dilute solutions of long-chain polymer molecules. The magnetic field imparts elasticity to the fluid.

3.2.3 Differential rotation and torsional Alfvén waves

We first consider the kinematic behaviour of a magnetic field in the presence of a prescribed velocity field involving differential rotation. In cylindrical polar coordinates $(r,{\it\phi},z)$ , let

(3.39) $$\begin{eqnarray}\boldsymbol{u}=r{\it\Omega}(r,z)\,\boldsymbol{e}_{{\it\phi}}.\end{eqnarray}$$

Consider an axisymmetric magnetic field, which we separate into poloidal (meridional: $r$ and $z$ ) and toroidal (azimuthal: ${\it\phi}$ ) parts:

(3.40) $$\begin{eqnarray}\boldsymbol{B}=\boldsymbol{B}_{p}(r,z,t)+B_{{\it\phi}}(r,z,t)\,\boldsymbol{e}_{{\it\phi}}.\end{eqnarray}$$

The ideal induction equation reduces to (exercise)

(3.41) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}_{p}}{\partial t}=0,\quad \frac{\partial B_{{\it\phi}}}{\partial t}=r\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Omega}.\end{eqnarray}$$

Differential rotation winds the poloidal field to generate a toroidal field. To obtain a steady state without winding, we require the angular velocity to be constant along each magnetic field line:

(3.42) $$\begin{eqnarray}\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Omega}=0,\end{eqnarray}$$

a result known as Ferraro’s law of isorotation ⁵ .

There is an energetic cost to winding the field, as work is done against magnetic tension. In a dynamical situation a strong magnetic field tends to enforce isorotation along its length.

We now generalize the analysis to allow for axisymmetric torsional oscillations:

(3.43) $$\begin{eqnarray}\boldsymbol{u}=r{\it\Omega}(r,z,t)\,\boldsymbol{e}_{{\it\phi}}.\end{eqnarray}$$

The azimuthal component of the equation of motion is (exercise)

(3.44) $$\begin{eqnarray}{\it\rho}r\frac{\partial {\it\Omega}}{\partial t}=\frac{1}{{\it\mu}_{0}r}\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(rB_{{\it\phi}}).\end{eqnarray}$$

This combines with the induction equation to give

(3.45) $$\begin{eqnarray}\frac{\partial ^{2}{\it\Omega}}{\partial t^{2}}=\frac{1}{{\it\mu}_{0}{\it\rho}r^{2}}\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(r^{2}\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Omega}).\end{eqnarray}$$

This equation describes torsional Alfvén waves. For example, if $\boldsymbol{B}_{p}=B_{z}\,\boldsymbol{e}_{z}$ is vertical and uniform, then

(3.46) $$\begin{eqnarray}\frac{\partial ^{2}{\it\Omega}}{\partial t^{2}}=v_{a}^{2}\frac{\partial ^{2}{\it\Omega}}{\partial z^{2}}.\end{eqnarray}$$

This is not strictly an exact nonlinear analysis because we have neglected the force balance (and indeed motion) in the meridional plane.

3.2.4 Force-free fields

In regions of low density, such as the solar corona, the magnetic field may be dynamically dominant over the effects of inertia, gravity and gas pressure. Under these circumstances we have (approximately) a force-free magnetic field such that

(3.47) $$\begin{eqnarray}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}=\mathbf{0}.\end{eqnarray}$$

Vector fields $\boldsymbol{B}$ satisfying this equation are known in a wider mathematical context as Beltrami fields. Since $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}$ must be parallel to $\boldsymbol{B}$ , we may write

(3.48) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}={\it\lambda}\boldsymbol{B},\end{eqnarray}$$

for some scalar field ${\it\lambda}(\boldsymbol{x})$ . The divergence of this equation is

(3.49) $$\begin{eqnarray}0=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\lambda},\end{eqnarray}$$

so that ${\it\lambda}$ is constant along each magnetic field line. In the special case ${\it\lambda}=\text{const.}$ , known as a linear force-free magnetic field, the curl of (3.48) results in the Helmholtz equation

(3.50) $$\begin{eqnarray}-{\rm\nabla}^{2}\boldsymbol{B}={\it\lambda}^{2}\boldsymbol{B},\end{eqnarray}$$

which admits a wide variety of non-trivial solutions.

A subset of force-free magnetic fields consists of potential or current-free magnetic fields for which

(3.51) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}=\mathbf{0}.\end{eqnarray}$$

In a true vacuum, the magnetic field must be potential. However, only an extremely low density of charge carriers (i.e. electrons) is needed to make the force-free description more relevant.

An example of a force-free field in cylindrical polar coordinates $(r,{\it\phi},z)$ is

(3.52) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{B}=B_{{\it\phi}}(r)\,\boldsymbol{e}_{{\it\phi}}+B_{z}(r)\,\boldsymbol{e}_{z},\\ \displaystyle \boldsymbol{{\rm\nabla}}\times \boldsymbol{B}=-\frac{\text{d}B_{z}}{\text{d}r}\,\boldsymbol{e}_{{\it\phi}}+\frac{1}{r}\frac{\text{d}}{\text{d}r}(rB_{{\it\phi}})\,\boldsymbol{e}_{z}.\end{array}\right\}\end{eqnarray}$$

Now $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}={\it\lambda}\boldsymbol{B}$ implies

(3.53) $$\begin{eqnarray}-\frac{1}{r}\frac{\text{d}}{\text{d}r}\left(r\frac{\text{d}B_{z}}{\text{d}r}\right)={\it\lambda}^{2}B_{z},\end{eqnarray}$$

which is the $z$ component of the Helmholtz equation. The solution regular at $r=0$ is

(3.54a,b ) $$\begin{eqnarray}B_{z}=B_{0}\text{J}_{0}({\it\lambda}r),\quad B_{{\it\phi}}=B_{0}\text{J}_{1}({\it\lambda}r),\end{eqnarray}$$

where $\text{J}_{n}$ is the Bessel function of order $n$ (figure 2). (Note that $\text{J}_{0}(x)$ satisfies $(x\text{J}_{0}^{\prime })^{\prime }+x\text{J}_{0}=0$ and $\text{J}_{1}(x)=-\text{J}_{0}^{\prime }(x)$ .) The helical nature of this field is typical of force-free fields with ${\it\lambda}\neq 0$ .

Figure 2. The Bessel functions $\text{J}_{0}(x)$ and $\text{J}_{1}(x)$ from the origin to the first zero of $\text{J}_{1}$ .

When applied to a infinite cylinder (e.g. as a simplified model of a magnetized astrophysical jet), the solution could be extended from the axis to the first zero of $\text{J}_{1}$ and then matched to a uniform external axial field $B_{z}$ . In this case the net axial current is zero. Alternatively the solution could be extended from the axis to the first zero of $\text{J}_{0}$ and matched to an external azimuthal field $B_{{\it\phi}}\propto r^{-1}$ generated by the net axial current.

3.2.5 Magnetostatic equilibrium and magnetic buoyancy

A magnetostatic equilibrium is a static solution $(\boldsymbol{u}=\mathbf{0})$ of the equation of motion, i.e. one satisfying

(3.55) $$\begin{eqnarray}\mathbf{0}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}-\boldsymbol{{\rm\nabla}}p+\frac{1}{{\it\mu}_{0}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B},\end{eqnarray}$$

together with $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$ .

While it is possible to find solutions in which the forces balance in this way, inhomogeneities in the magnetic field typically result in a lack of equilibrium. A magnetic flux tube (figure 3) is an idealized situation in which the magnetic field is localized to the interior of a tube and vanishes outside. To balance the total pressure at the interface, the gas pressure must be lower inside. Unless the temperatures are different, the density is lower inside. In a gravitational field the tube therefore experiences an upward buoyancy force and tends to rise.

Figure 3. A buoyant magnetic flux tube.

Related examples: A.7–A.9.

4 Conservation laws, symmetries and hyperbolic structure

4.1 Introduction

There are various ways in which a quantity can be said to be ‘conserved’ in fluid dynamics or MHD. If a quantity has a density (amount per unit volume) $q(\boldsymbol{x},t)$ that satisfies an equation of the conservative form

(4.1) $$\begin{eqnarray}\frac{\partial q}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{F}=0,\end{eqnarray}$$

then the vector field $\boldsymbol{F}(\boldsymbol{x},t)$ can be identified as the flux density (flux per unit area) of the quantity. The rate of change of the total amount of the quantity in a time-independent volume $V$ ,

(4.2) $$\begin{eqnarray}Q=\int _{V}q\,\text{d}V,\end{eqnarray}$$

is then equal to minus the flux of $\boldsymbol{F}$ through the bounding surface $S$ :

(4.3) $$\begin{eqnarray}\frac{\text{d}Q}{\text{d}t}=-\int _{V}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{F})\,\text{d}V=-\int _{S}\boldsymbol{F}\boldsymbol{\cdot }\text{d}\boldsymbol{S}.\end{eqnarray}$$

If the boundary conditions on $S$ are such that this flux vanishes, then $Q$ is constant; otherwise, changes in $Q$ can be accounted for by the flux of $\boldsymbol{F}$ through $S$ . In this sense the quantity is said to be conserved. The prototype is mass, for which $q={\it\rho}$ and $\boldsymbol{F}={\it\rho}\boldsymbol{u}$ .

A material invariant is a scalar field $f(\boldsymbol{x},t)$ for which

(4.4) $$\begin{eqnarray}\frac{\text{D}f}{\text{D}t}=0,\end{eqnarray}$$

which implies that $f$ is constant for each fluid element, and is therefore conserved following the fluid motion. A simple example is the specific entropy in ideal fluid dynamics. When combined with mass conservation, this yields an equation in conservative form,

(4.5) $$\begin{eqnarray}\frac{\partial }{\partial t}({\it\rho}f)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}f\boldsymbol{u})=0.\end{eqnarray}$$

4.2 Synthesis of the total energy equation

Starting from the ideal MHD equations, we construct the total energy equation piece by piece.

Kinetic energy:

(4.6) $$\begin{eqnarray}{\it\rho}\frac{\text{D}}{\text{D}t}\left(\frac{1}{2}u^{2}\right)={\it\rho}\boldsymbol{u}\boldsymbol{\cdot }\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-{\it\rho}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Phi}-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p+\frac{1}{{\it\mu}_{0}}\boldsymbol{u}\boldsymbol{\cdot }\left[(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}\right].\end{eqnarray}$$

Gravitational energy (assuming initially that the system is non-self-gravitating and that ${\it\Phi}$ is independent of $t$ ):

(4.7) $$\begin{eqnarray}{\it\rho}\frac{\text{D}{\it\Phi}}{\text{D}t}={\it\rho}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Phi}.\end{eqnarray}$$

Internal (thermal) energy (using the fundamental thermodynamic identity $\text{d}e=T\,\text{d}s-p\,\text{d}v$ ):

(4.8) $$\begin{eqnarray}{\it\rho}\frac{\text{D}e}{\text{D}t}={\it\rho}T\frac{\text{D}s}{\text{D}t}+p\frac{\text{D}\ln {\it\rho}}{\text{D}t}=-p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}.\end{eqnarray}$$

Sum of these three:

(4.9) $$\begin{eqnarray}{\it\rho}\frac{\text{D}}{\text{D}t}\left(\frac{1}{2}u^{2}+{\it\Phi}+e\right)=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(p\boldsymbol{u})+\frac{1}{{\it\mu}_{0}}\boldsymbol{u}\boldsymbol{\cdot }\left[(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}\right].\end{eqnarray}$$

The last term can be rewritten as

(4.10) $$\begin{eqnarray}\frac{1}{{\it\mu}_{0}}\boldsymbol{u}\boldsymbol{\cdot }\left[(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}\right]=\frac{1}{{\it\mu}_{0}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\boldsymbol{\cdot }(-\boldsymbol{u}\times \boldsymbol{B})=\frac{1}{{\it\mu}_{0}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\boldsymbol{\cdot }\boldsymbol{E}.\end{eqnarray}$$

Using mass conservation:

(4.11) $$\begin{eqnarray}\frac{\partial }{\partial t}\left[{\it\rho}\left(\frac{1}{2}u^{2}+{\it\Phi}+e\right)\right]+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\rho}\boldsymbol{u}\left(\frac{1}{2}u^{2}+{\it\Phi}+e\right)+p\boldsymbol{u}\right]=\frac{1}{{\it\mu}_{0}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\boldsymbol{\cdot }\boldsymbol{E}.\end{eqnarray}$$

Magnetic energy:

(4.12) $$\begin{eqnarray}\frac{\partial }{\partial t}\left(\frac{B^{2}}{2{\it\mu}_{0}}\right)=\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }\frac{\partial \boldsymbol{B}}{\partial t}=-\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\times \boldsymbol{E}.\end{eqnarray}$$

Total energy:

(4.13) $$\begin{eqnarray}\frac{\partial }{\partial t}\left[{\it\rho}\left(\frac{1}{2}u^{2}+{\it\Phi}+e\right)+\frac{B^{2}}{2{\it\mu}_{0}}\right]+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\rho}\boldsymbol{u}\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)+\frac{\boldsymbol{E}\times \boldsymbol{B}}{{\it\mu}_{0}}\right]=0,\end{eqnarray}$$

where $h=e+p/{\it\rho}$ is the specific enthalpy and we have used the identity $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{E}\times \boldsymbol{B})=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\times \boldsymbol{E}-\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}$ . Note that $(\boldsymbol{E}\times \boldsymbol{B})/{\it\mu}_{0}$ is the Poynting vector, the electromagnetic energy flux density. The total energy is therefore conserved.

For a self-gravitating system satisfying Poisson’s equation, the gravitational energy density can instead be regarded as $-g^{2}/8{\rm\pi}G$ :

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(-\frac{g^{2}}{8{\rm\pi}G}\right)=-\frac{1}{4{\rm\pi}G}\boldsymbol{{\rm\nabla}}{\it\Phi}\boldsymbol{\cdot }\frac{\partial \boldsymbol{{\rm\nabla}}{\it\Phi}}{\partial t} & \displaystyle\end{eqnarray}$$

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(-\frac{g^{2}}{8{\rm\pi}G}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{{\it\Phi}}{4{\rm\pi}G}\frac{\partial \boldsymbol{{\rm\nabla}}{\it\Phi}}{\partial t}\right)=\frac{{\it\Phi}}{4{\rm\pi}G}\frac{\partial {\rm\nabla}^{2}{\it\Phi}}{\partial t}={\it\Phi}\frac{\partial {\it\rho}}{\partial t}=-{\it\Phi}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u}) & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(-\frac{g^{2}}{8{\rm\pi}G}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left({\it\rho}\boldsymbol{u}{\it\Phi}+\frac{{\it\Phi}}{4{\rm\pi}G}\frac{\partial \boldsymbol{{\rm\nabla}}{\it\Phi}}{\partial t}\right)={\it\rho}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Phi}. & \displaystyle\end{eqnarray}$$

The total energy equation is then

(4.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\partial }{\partial t}\left[{\it\rho}\left(\frac{1}{2}u^{2}+e\right)-\frac{g^{2}}{8{\rm\pi}G}+\frac{B^{2}}{2{\it\mu}_{0}}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\rho}\boldsymbol{u}\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)+\frac{{\it\Phi}}{4{\rm\pi}G}\frac{\partial \boldsymbol{{\rm\nabla}}{\it\Phi}}{\partial t}+\frac{\boldsymbol{E}\times \boldsymbol{B}}{{\it\mu}_{0}}\right]=0.\end{eqnarray}$$

It is important to note that some of the gravitational and magnetic energy of an astrophysical body is stored in the exterior region, even if the mass density vanishes there.

4.3 Other conservation laws in ideal MHD

In ideal fluid dynamics there are certain invariants with a geometrical or topological interpretation. In homentropic or barotropic flow, for example, vorticity (or, equivalently, circulation) and kinetic helicity are conserved, while, in non-barotropic flow, potential vorticity is conserved (see Example A.2). The Lorentz force breaks these conservation laws because the curl of the Lorentz force per unit mass does not vanish in general. However, some new topological invariants associated with the magnetic field appear.

The magnetic helicity in a volume $V$ with bounding surface $S$ is defined as

(4.18) $$\begin{eqnarray}H_{m}=\int _{V}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}V,\end{eqnarray}$$

where $\boldsymbol{A}$ is the magnetic vector potential, such that $\boldsymbol{B}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{A}$ . Now

(4.19) $$\begin{eqnarray}\frac{\partial \boldsymbol{A}}{\partial t}=-\boldsymbol{E}-\boldsymbol{{\rm\nabla}}{\it\Phi}_{e}=\boldsymbol{u}\times \boldsymbol{B}-\boldsymbol{{\rm\nabla}}{\it\Phi}_{e},\end{eqnarray}$$

where ${\it\Phi}_{e}$ is the electrostatic potential. This can be thought of as the ‘uncurl’ of the induction equation. Thus

(4.20) $$\begin{eqnarray}\frac{\partial }{\partial t}(\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B})=-\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Phi}_{e}+\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B}).\end{eqnarray}$$

In ideal MHD, therefore, magnetic helicity is conserved:

(4.21) $$\begin{eqnarray}\frac{\partial }{\partial t}(\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\Phi}_{e}\boldsymbol{B}+\boldsymbol{A}\times (\boldsymbol{u}\times \boldsymbol{B})\right]=0.\end{eqnarray}$$

However, care is needed because $\boldsymbol{A}$ is not uniquely defined. Under a gauge transformation $\boldsymbol{A}\mapsto \boldsymbol{A}+\boldsymbol{{\rm\nabla}}{\it\chi}$ , ${\it\Phi}_{e}\mapsto {\it\Phi}_{e}-\partial {\it\chi}/\partial t$ , where ${\it\chi}(\boldsymbol{x},t)$ is a scalar field, $\boldsymbol{E}$ and $\boldsymbol{B}$ are invariant, but $H_{m}$ changes by an amount

(4.22) $$\begin{eqnarray}\int _{V}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\chi}\,\text{d}V=\int _{V}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\chi}\boldsymbol{B})\,\text{d}V=\int _{S}{\it\chi}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S.\end{eqnarray}$$

Therefore $H_{m}$ is not uniquely defined unless $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ on the surface $S$ .

Magnetic helicity is a pseudoscalar quantity: it changes sign under a reflection of the spatial coordinates. Indeed, it is non-zero only when the magnetic field lacks reflectional symmetry. It can also be interpreted topologically in terms of the twistedness and knottedness of the magnetic field (see Example A.10). Since the field is ‘frozen in’ to the fluid and deformed continuously by it, the topological properties of the field are conserved. The equivalent conserved quantity in homentropic or barotropic ideal gas dynamics (without a magnetic field) is the kinetic helicity

(4.23) $$\begin{eqnarray}H_{k}=\int _{V}\boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u})\,\text{d}V.\end{eqnarray}$$

The cross-helicity in a volume $V$ is

(4.24) $$\begin{eqnarray}H_{c}=\int _{V}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}V.\end{eqnarray}$$

It is helpful here to write the equation of motion in ideal MHD in the form

(4.25) $$\begin{eqnarray}\frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u})\times \boldsymbol{u}+\boldsymbol{{\rm\nabla}}\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)=T\boldsymbol{{\rm\nabla}}s+\frac{1}{{\it\mu}_{0}{\it\rho}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B},\end{eqnarray}$$

using the relation $\text{d}h=T\,\text{d}s+v\,\text{d}p$ . Thus

(4.26) $$\begin{eqnarray}\frac{\partial }{\partial t}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[\boldsymbol{u}\times (\boldsymbol{u}\times \boldsymbol{B})+\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)\boldsymbol{B}\right]=T\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s,\end{eqnarray}$$

and so cross-helicity is conserved in ideal MHD in homentropic or barotropic flow.

Bernoulli’s theorem follows from the inner product of (4.25) with $\boldsymbol{u}$ . In steady flow

(4.27) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)=0,\end{eqnarray}$$

which implies that the Bernoulli function $(1/2)u^{2}+{\it\Phi}+h$ is constant along streamlines, but only if $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{F}_{m}=0$ (e.g. if $\boldsymbol{u}\,\Vert \,\boldsymbol{B}$ ), i.e. if the magnetic field does no work on the flow.

Related examples: A.10, A.11.

4.4 Symmetries of the equations

The equations of ideal gas dynamics and MHD have numerous symmetries. In the case of an isolated, self-gravitating system, these include:

1. (i) Translations of time and space, and rotations of space: related (via Noether’s theorem) to the conservation of energy, momentum and angular momentum.

2. (ii) Reversal of time: related to the absence of dissipation.

3. (iii) Reflections of space (but note that $\boldsymbol{B}$ is a pseudovector and behaves oppositely to  $\boldsymbol{u}$ under a reflection).

4. (iv) Galilean transformations.

5. (v) Reversal of the sign of $\boldsymbol{B}$ .

6. (vi) Similarity transformations (exercise): if space and time are rescaled by independent factors ${\it\lambda}$ and ${\it\mu}$ , i.e.

   (4.28a,b ) $$\begin{eqnarray}\boldsymbol{x}\mapsto {\it\lambda}\,\boldsymbol{x},\quad t\mapsto {\it\mu}\,t,\end{eqnarray}$$
   then
   (4.29a-e ) $$\begin{eqnarray}\boldsymbol{u}\mapsto {\it\lambda}{\it\mu}^{-1}\,\boldsymbol{u},\quad {\it\rho}\mapsto {\it\mu}^{-2}\,{\it\rho},\quad p\mapsto {\it\lambda}^{2}{\it\mu}^{-4}\,p,\quad {\it\Phi}\mapsto {\it\lambda}^{2}{\it\mu}^{-2}\,{\it\Phi},\quad \boldsymbol{B}\mapsto {\it\lambda}{\it\mu}^{-2}\,\boldsymbol{B}.\end{eqnarray}$$
   (This symmetry requires a perfect gas so that the thermodynamic relations are scale free.)

In the case of a non-isolated system with an external potential ${\it\Phi}_{\mathit{ext}}$ , these symmetries (other than $\boldsymbol{B}\mapsto -\boldsymbol{B}$ ) apply only if ${\it\Phi}_{\mathit{ext}}$ has them. However, in the approximation of a non-self-gravitating system, the mass can be rescaled by any factor ${\it\lambda}$ such that

(4.30a-c ) $$\begin{eqnarray}{\it\rho}\mapsto {\it\lambda}\,{\it\rho},\quad p\mapsto {\it\lambda}\,p,\quad \boldsymbol{B}\mapsto {\it\lambda}^{1/2}\,\boldsymbol{B}.\end{eqnarray}$$

(This symmetry also requires a perfect gas.)

4.5 Hyperbolic structure

Analysing the so-called hyperbolic structure of the equations of AFD is one way of understanding the wave modes of the system and the way in which information propagates in the fluid. It is fundamental to the construction of some types of numerical method for solving the equations. We temporarily neglect the gravitational force here, because in a Newtonian theory it involves instantaneous action at a distance and is not associated with a finite wave speed.

In ideal gas dynamics, the equation of mass conservation, the thermal energy equation and the equation of motion (omitting gravity) can be written as

(4.31) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial {\it\rho}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}+{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\\ \displaystyle \frac{\partial p}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p+{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\\ \displaystyle \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\frac{1}{{\it\rho}}\boldsymbol{{\rm\nabla}}p=\mathbf{0},\end{array}\right\}\end{eqnarray}$$

and then combined into the form

(4.32) $$\begin{eqnarray}\frac{\partial \boldsymbol{U}}{\partial t}+\unicode[STIX]{x1D63C}_{i}\frac{\partial \boldsymbol{U}}{\partial x_{i}}=\mathbf{0},\end{eqnarray}$$

where

(4.33) $$\begin{eqnarray}\boldsymbol{U}=\left[\begin{array}{@{}c@{}}{\it\rho}\\ p\\ u_{x}\\ u_{y}\\ u_{z}\end{array}\right]\end{eqnarray}$$

is a five-dimensional ‘state vector’ and $\unicode[STIX]{x1D63C}_{x}$ , $\unicode[STIX]{x1D63C}_{y}$ and $\unicode[STIX]{x1D63C}_{z}$ are the three $5\times 5$ matrices

(4.34a-c ) $$\begin{eqnarray}\left[\begin{array}{@{}ccccc@{}}u_{x} & 0 & {\it\rho} & 0 & 0\\ 0 & u_{x} & {\it\gamma}p & 0 & 0\\ 0 & \frac{1}{{\it\rho}} & u_{x} & 0 & 0\\ 0 & 0 & 0 & u_{x} & 0\\ 0 & 0 & 0 & 0 & u_{x}\end{array}\right],\quad \left[\begin{array}{@{}ccccc@{}}u_{y} & 0 & 0 & {\it\rho} & 0\\ 0 & u_{y} & 0 & {\it\gamma}p & 0\\ 0 & 0 & u_{y} & 0 & 0\\ 0 & \frac{1}{{\it\rho}} & 0 & u_{y} & 0\\ 0 & 0 & 0 & 0 & u_{y}\end{array}\right],\quad \left[\begin{array}{@{}ccccc@{}}u_{z} & 0 & 0 & 0 & {\it\rho}\\ 0 & u_{z} & 0 & 0 & {\it\gamma}p\\ 0 & 0 & u_{z} & 0 & 0\\ 0 & 0 & 0 & u_{z} & 0\\ 0 & \frac{1}{{\it\rho}} & 0 & 0 & u_{z}\end{array}\right].\end{eqnarray}$$

This works because every term in the equations involves a first derivative with respect to either time or space.

The system of equations is said to be hyperbolic if the eigenvalues of $\unicode[STIX]{x1D63C}_{i}n_{i}$ are real for any unit vector $\boldsymbol{n}$ and if the eigenvectors span the five-dimensional space. As will be seen in § 6.2, the eigenvalues can be identified as wave speeds, and the eigenvectors as wave modes, with $\boldsymbol{n}$ being the unit wavevector, locally normal to the wavefronts.

Taking $\boldsymbol{n}=\boldsymbol{e}_{x}$ without loss of generality, we find (exercise)

(4.35) $$\begin{eqnarray}\det (\unicode[STIX]{x1D63C}_{x}-v\unicode[STIX]{x1D644})=-(v-u_{x})^{3}\left[(v-u_{x})^{2}-v_{s}^{2}\right],\end{eqnarray}$$

where

(4.36) $$\begin{eqnarray}v_{s}=\left(\frac{{\it\gamma}p}{{\it\rho}}\right)^{1/2}\end{eqnarray}$$

is the adiabatic sound speed. The wave speeds $v$ are real and the system is indeed hyperbolic.

Two of the wave modes are sound waves (acoustic waves), which have wave speeds $v=u_{x}\pm v_{s}$ and therefore propagate at the sound speed relative to the moving fluid. Their eigenvectors are

(4.37) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}{\it\rho}\\ {\it\gamma}p\\ \pm v_{s}\\ 0\\ 0\end{array}\right]\end{eqnarray}$$

and involve perturbations of density, pressure and longitudinal velocity.

The remaining three wave modes have wave speed $v=u_{x}$ and do not propagate relative to the fluid. Their eigenvectors are

(4.38a-c ) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}1\\ 0\\ 0\\ 0\\ 0\end{array}\right],\quad \left[\begin{array}{@{}c@{}}0\\ 0\\ 0\\ 1\\ 0\end{array}\right],\quad \left[\begin{array}{@{}c@{}}0\\ 0\\ 0\\ 0\\ 1\end{array}\right].\end{eqnarray}$$

The first is the entropy wave, which involves only a density perturbation but no pressure perturbation. Since the entropy can be considered as a function of the density and pressure, this wave involves an entropy perturbation. It must therefore propagate at the fluid velocity because the entropy is a material invariant. The other two modes with $v=u_{x}$ are vortical waves, which involve perturbations of the transverse velocity components, and therefore of the vorticity. Conservation of vorticity implies that these waves propagate with the fluid velocity.

To extend the analysis to ideal MHD, we may consider the induction equation in the form

(4.39) $$\begin{eqnarray}\frac{\partial \boldsymbol{B}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}-\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{B}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u})=\mathbf{0},\end{eqnarray}$$

and include the Lorentz force in the equation of motion. Every new term involves a first derivative. So the equation of mass conservation, the thermal energy equation, the equation of motion and the induction equation can be written in the combined form

(4.40) $$\begin{eqnarray}\frac{\partial \boldsymbol{U}}{\partial t}+\unicode[STIX]{x1D63C}_{i}\frac{\partial \boldsymbol{U}}{\partial x_{i}}=\mathbf{0},\end{eqnarray}$$

where

(4.41) $$\begin{eqnarray}\boldsymbol{U}=\left[\begin{array}{@{}c@{}}{\it\rho}\\ p\\ u_{x}\\ u_{y}\\ u_{z}\\ B_{x}\\ B_{y}\\ B_{z}\end{array}\right]\end{eqnarray}$$

is now an eight-dimensional ‘state vector’ and the $\unicode[STIX]{x1D63C}_{i}$ are three $8\times 8$ matrices, e.g.

(4.42) $$\begin{eqnarray}\unicode[STIX]{x1D63C}_{x}=\left[\begin{array}{@{}cccccccc@{}}u_{x} & 0 & {\it\rho} & 0 & 0 & 0 & 0 & 0\\ 0 & u_{x} & {\it\gamma}p & 0 & 0 & 0 & 0 & 0\\ 0 & {\displaystyle \frac{1}{{\it\rho}}} & u_{x} & 0 & 0 & 0 & {\displaystyle \frac{B_{y}}{{\it\mu}_{0}{\it\rho}}} & {\displaystyle \frac{B_{z}}{{\it\mu}_{0}{\it\rho}}}\\ 0 & 0 & 0 & u_{x} & 0 & 0 & -{\displaystyle \frac{B_{x}}{{\it\mu}_{0}{\it\rho}}} & 0\\ 0 & 0 & 0 & 0 & u_{x} & 0 & 0 & -{\displaystyle \frac{B_{x}}{{\it\mu}_{0}{\it\rho}}}\\ 0 & 0 & 0 & 0 & 0 & u_{x} & 0 & 0\\ 0 & 0 & B_{y} & -B_{x} & 0 & 0 & u_{x} & 0\\ 0 & 0 & B_{z} & 0 & -B_{x} & 0 & 0 & u_{x}\end{array}\right].\end{eqnarray}$$

We now find, after some algebra,

(4.43) $$\begin{eqnarray}\det (\unicode[STIX]{x1D63C}_{x}-v\unicode[STIX]{x1D644})=(v-u_{x})^{2}\left[(v-u_{x})^{2}-v_{ax}^{2}\right]\left[(v-u_{x})^{4}-(v_{s}^{2}+v_{a}^{2})(v-u_{x})^{2}+v_{s}^{2}v_{ax}^{2}\right].\end{eqnarray}$$

The wave speeds $v$ are real and the system is indeed hyperbolic. The various MHD wave modes will be examined later (§ 5).

In this representation, there are two modes that have $v=u_{x}$ and do not propagate relative to the fluid. One is still the entropy wave, which is physical and involves only a density perturbation. The other is the ‘ $\text{div}\boldsymbol{B}$ ’ mode, which is unphysical and involves a perturbation of $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}$ (i.e. of $B_{x}$ , in the case $\boldsymbol{n}=\boldsymbol{e}_{x}$ ). This must be eliminated by imposing the constraint $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$ . (In fact the equations in the form we have written them imply that $(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B})/{\it\rho}$ is a material invariant and could be non-zero unless the initial condition $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$ is imposed.) The vortical waves are replaced by Alfvén waves with speeds $u_{x}\pm v_{ax}$ .

4.6 Stress tensor and virial theorem

In the absence of external forces, the equation of motion of a fluid can usually be written in the form

(4.44a,b ) $$\begin{eqnarray}{\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\quad \text{or}\quad {\it\rho}\frac{\text{D}u_{i}}{\text{D}t}=\frac{\partial T_{ji}}{\partial x_{j}},\end{eqnarray}$$

where $\unicode[STIX]{x1D64F}$ is the stress tensor, a symmetric second-rank tensor field. Using the equation of mass conservation, we can relate this to the conservative form of the momentum equation,

(4.45) $$\begin{eqnarray}\frac{\partial }{\partial t}({\it\rho}\boldsymbol{u})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u}\boldsymbol{u}-\unicode[STIX]{x1D64F})=\mathbf{0},\end{eqnarray}$$

which shows that $-\unicode[STIX]{x1D64F}$ is the momentum flux density excluding the advective flux of momentum.

For a self-gravitating system in ideal MHD, the stress tensor is

(4.46) $$\begin{eqnarray}\unicode[STIX]{x1D64F}=-p\,\unicode[STIX]{x1D644}-\frac{1}{4{\rm\pi}G}\left(\boldsymbol{g}\boldsymbol{g}-\frac{1}{2}g^{2}\,\unicode[STIX]{x1D644}\right)+\frac{1}{{\it\mu}_{0}}\left(\boldsymbol{B}\boldsymbol{B}-\frac{1}{2}B^{2}\,\unicode[STIX]{x1D644}\right),\end{eqnarray}$$

or, in Cartesian components,

(4.47) $$\begin{eqnarray}T_{ij}=-p\,{\it\delta}_{ij}-\frac{1}{4{\rm\pi}G}\left(g_{i}g_{j}-\frac{1}{2}g^{2}\,{\it\delta}_{ij}\right)+\frac{1}{{\it\mu}_{0}}\left(B_{i}B_{j}-\frac{1}{2}B^{2}\,{\it\delta}_{ij}\right).\end{eqnarray}$$

We have already identified the Maxwell stress tensor associated with the magnetic field. The idea of a gravitational stress tensor works for a self-gravitating system in which the gravitational field $\boldsymbol{g}=-\boldsymbol{{\rm\nabla}}{\it\Phi}$ and the density ${\it\rho}$ are related through Poisson’s equation $-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{g}={\rm\nabla}^{2}{\it\Phi}=4{\rm\pi}G{\it\rho}$ . In fact, for a general vector field $\boldsymbol{v}$ , it can be shown that (exercise)

(4.48) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{v}\boldsymbol{v}-{\textstyle \frac{1}{2}}v^{2}\,\unicode[STIX]{x1D644})=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})\boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v}-\boldsymbol{{\rm\nabla}}\left({\textstyle \frac{1}{2}}v^{2}\right)=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})\boldsymbol{v}+(\boldsymbol{{\rm\nabla}}\times \boldsymbol{v})\times \boldsymbol{v}.\end{eqnarray}$$

In the magnetic case ( $\boldsymbol{v}=\boldsymbol{B}$ ) the first term in the final expression vanishes, while in the gravitational case ( $\boldsymbol{v}=\boldsymbol{g}$ ) the second term vanishes, leaving $-4{\rm\pi}G{\it\rho}\boldsymbol{g}$ , which becomes the force per unit volume, ${\it\rho}\boldsymbol{g}$ , when divided by $-4{\rm\pi}G$ .

The virial equations are the spatial moments of the equation of motion, and provide integral measures of the balance of forces acting on the fluid. The first moments are generally the most useful. Consider

(4.49) $$\begin{eqnarray}{\it\rho}\frac{\text{D}^{2}}{\text{D}t^{2}}(x_{i}x_{j})={\it\rho}\frac{\text{D}}{\text{D}t}(u_{i}x_{j}+x_{i}u_{j})=2{\it\rho}u_{i}u_{j}+x_{j}\frac{\partial T_{ki}}{\partial x_{k}}+x_{i}\frac{\partial T_{kj}}{\partial x_{k}}.\end{eqnarray}$$

Integrate this equation over a material volume $V$ bounded by a surface $S$ (with material invariant mass element $\text{d}m={\it\rho}\,\text{d}V$ ):

(4.50) $$\begin{eqnarray}\displaystyle \frac{\text{d}^{2}}{\text{d}t^{2}}\int _{V}x_{i}x_{j}\,\text{d}m & = & \displaystyle \int _{V}\left(2{\it\rho}u_{i}u_{j}+x_{j}\frac{\partial T_{ki}}{\partial x_{k}}+x_{i}\frac{\partial T_{kj}}{\partial x_{k}}\right)\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{V}(2{\it\rho}u_{i}u_{j}-T_{ji}-T_{ij})\,\text{d}V+\int _{S}(x_{j}T_{ki}+x_{i}T_{kj})n_{k}\,\text{d}S,\end{eqnarray}$$

where we have integrated by parts using the divergence theorem. In the case of an isolated system with no external sources of gravity or magnetic field, $\boldsymbol{g}$ decays proportional to $|\boldsymbol{x}|^{-2}$ at large distance, and $\boldsymbol{B}$ decays faster. Therefore $T_{ij}$ decays proportional to $|\boldsymbol{x}|^{-4}$ and the surface integral can be eliminated if we let $V$ occupy the whole of space. We then obtain (after division by  $2$ ) the tensor virial theorem

(4.51) $$\begin{eqnarray}\frac{1}{2}\frac{\text{d}^{2}\unicode[STIX]{x1D610}_{ij}}{\text{d}t^{2}}=2\unicode[STIX]{x1D612}_{ij}-{\mathcal{T}}_{ij},\end{eqnarray}$$

where

(4.52) $$\begin{eqnarray}\unicode[STIX]{x1D610}_{ij}=\int x_{i}x_{j}\,\text{d}m\end{eqnarray}$$

is related to the inertia tensor of the system,

(4.53) $$\begin{eqnarray}\unicode[STIX]{x1D612}_{ij}=\int \frac{1}{2}u_{i}u_{j}\,\text{d}m\end{eqnarray}$$

is a kinetic energy tensor and

(4.54) $$\begin{eqnarray}{\mathcal{T}}_{ij}=\int T_{ij}\,\text{d}V\end{eqnarray}$$

is the integrated stress tensor. (If the conditions above are not satisfied, there will be an additional contribution from the surface integral.)

The scalar virial theorem is the trace of this tensor equation, which we write as

(4.55) $$\begin{eqnarray}\frac{1}{2}\frac{\text{d}^{2}I}{\text{d}t^{2}}=2K-{\mathcal{T}}.\end{eqnarray}$$

Note that $K$ is the total kinetic energy. Now

(4.56) $$\begin{eqnarray}-{\mathcal{T}}=\int \left(3p-\frac{g^{2}}{8{\rm\pi}G}+\frac{B^{2}}{2{\it\mu}_{0}}\right)\,\text{d}V=3({\it\gamma}-1)U+W+M,\end{eqnarray}$$

for a perfect gas with no external gravitational field, where $U$ , $W$ and $M$ are the total internal, gravitational and magnetic energies. Thus

(4.57) $$\begin{eqnarray}\frac{1}{2}\frac{\text{d}^{2}I}{\text{d}t^{2}}=2K+3({\it\gamma}-1)U+W+M.\end{eqnarray}$$

On the right-hand side, only $W$ is negative. For the system to be bound (i.e. not fly apart) the kinetic, internal and magnetic energies are limited by

(4.58) $$\begin{eqnarray}2K+3({\it\gamma}-1)U+M\leqslant |W|.\end{eqnarray}$$

In fact equality must hold, at least on average, unless the system is collapsing or contracting.

The tensor virial theorem provides more specific information relating to the energies associated with individual directions, and is particularly relevant in cases where anisotropy is introduced by rotation or a magnetic field. It has been used in estimating the conditions required for gravitational collapse in molecular clouds. A higher-order tensor virial method was used by Chandrasekhar and Lebovitz to study the equilibrium and stability of rotating ellipsoidal bodies (Chandrasekhar 1969).

5 Linear waves in homogeneous media

In ideal MHD the density, pressure and magnetic field evolve according to

(5.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial {\it\rho}}{\partial t}=-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u},\\ \displaystyle \frac{\partial p}{\partial t}=-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u},\\ \displaystyle \frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B}).\end{array}\right\}\end{eqnarray}$$

Consider a magnetostatic equilibrium in which the density, pressure and magnetic field are ${\it\rho}_{0}(\boldsymbol{x})$ , $p_{0}(\boldsymbol{x})$ and $\boldsymbol{B}_{0}(\boldsymbol{x})$ . The above equations are exactly satisfied in this basic state because $\boldsymbol{u}=\mathbf{0}$ and the time derivatives vanish. Now consider small perturbations from equilibrium, such that ${\it\rho}(\boldsymbol{x},t)={\it\rho}_{0}(\boldsymbol{x})+{\it\delta}{\it\rho}(\boldsymbol{x},t)$ with $|{\it\delta}{\it\rho}|\ll {\it\rho}_{0}$ , etc. The linearized equations are

(5.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial {\it\delta}{\it\rho}}{\partial t}=-{\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}_{0}-{\it\rho}_{0}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\delta}\boldsymbol{u},\\ \displaystyle \frac{\partial {\it\delta}p}{\partial t}=-{\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p_{0}-{\it\gamma}p_{0}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\delta}\boldsymbol{u},\\ \displaystyle \frac{\partial {\it\delta}\boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times ({\it\delta}\boldsymbol{u}\times \boldsymbol{B}_{0}).\end{array}\right\}\end{eqnarray}$$

By introducing the displacement ${\bf\xi}(\boldsymbol{x},t)$ such that ${\it\delta}\boldsymbol{u}=\partial {\bf\xi}/\partial t$ , we can integrate these equations to obtain

(5.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\delta}{\it\rho}=-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi},\\ \displaystyle {\it\delta}p=-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi},\\ \displaystyle {\it\delta}\boldsymbol{B}=\boldsymbol{{\rm\nabla}}\times ({\bf\xi}\times \boldsymbol{B})=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}-\boldsymbol{B}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi})-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}.\end{array}\right\}\end{eqnarray}$$

We have now dropped the subscript ‘0’ without danger of confusion.

(The above relations allow some freedom to add arbitrary functions of $\boldsymbol{x}$ . At least when studying wave-like solutions in which all variables have the same harmonic time dependence, such additional terms can be discarded.)

The linearized equation of motion is

(5.4) $$\begin{eqnarray}{\it\rho}\frac{\partial ^{2}{\bf\xi}}{\partial t^{2}}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}-{\it\delta}{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}-\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Pi}+\frac{1}{{\it\mu}_{0}}({\it\delta}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}+\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\delta}\boldsymbol{B}),\end{eqnarray}$$

where the perturbation of total pressure is

(5.5) $$\begin{eqnarray}{\it\delta}{\it\Pi}={\it\delta}p+\frac{\boldsymbol{B}\boldsymbol{\cdot }{\it\delta}\boldsymbol{B}}{{\it\mu}_{0}}=-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Pi}-\left({\it\gamma}p+\frac{B^{2}}{{\it\mu}_{0}}\right)\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}+\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }(\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}).\end{eqnarray}$$

The gravitational potential perturbation satisfies the linearized Poisson equation

(5.6) $$\begin{eqnarray}{\rm\nabla}^{2}{\it\delta}{\it\Phi}=4{\rm\pi}G\,{\it\delta}{\it\rho}.\end{eqnarray}$$

We consider a basic state of uniform density, pressure and magnetic field, in the absence of gravity. Such a system is homogeneous but anisotropic, because the uniform field distinguishes a particular direction. The problem simplifies to

(5.7) $$\begin{eqnarray}{\it\rho}\frac{\partial ^{2}{\bf\xi}}{\partial t^{2}}=-\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Pi}+\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\left[\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}-\boldsymbol{B}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi})\right],\end{eqnarray}$$

with

(5.8) $$\begin{eqnarray}{\it\delta}{\it\Pi}=-\left({\it\gamma}p+\frac{B^{2}}{{\it\mu}_{0}}\right)\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}+\frac{1}{{\it\mu}_{0}}\boldsymbol{B}\boldsymbol{\cdot }(\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}).\end{eqnarray}$$

Owing to the symmetries of the basic state, plane-wave solutions exist, of the form

(5.9) $$\begin{eqnarray}{\bf\xi}(\boldsymbol{x},t)=\text{Re}[\tilde{{\bf\xi}}\,\exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}-\text{i}{\it\omega}t)],\end{eqnarray}$$

where ${\it\omega}$ and $\boldsymbol{k}$ are the frequency and wavevector, and $\tilde{{\bf\xi}}$ is a constant vector representing the amplitude of the wave. For such solutions, (5.7) gives

(5.10) $$\begin{eqnarray}{\it\rho}{\it\omega}^{2}{\bf\xi}=\left[\left({\it\gamma}p+\frac{B^{2}}{{\it\mu}_{0}}\right)\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi}-\frac{1}{{\it\mu}_{0}}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B})\boldsymbol{B}\boldsymbol{\cdot }{\bf\xi}\right]\boldsymbol{k}+\frac{1}{{\it\mu}_{0}}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B})\left[(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B}){\bf\xi}-\boldsymbol{B}(\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi})\right],\end{eqnarray}$$

where we have changed the sign and omitted the tilde.

For transverse displacements that are orthogonal to both the wavevector and the magnetic field, i.e.  $\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi}=\boldsymbol{B}\boldsymbol{\cdot }{\bf\xi}=0$ , this equation simplifies to

(5.11) $$\begin{eqnarray}{\it\rho}{\it\omega}^{2}{\bf\xi}=\frac{1}{{\it\mu}_{0}}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B})^{2}{\bf\xi}.\end{eqnarray}$$

Such solutions are called Alfvén waves. Their dispersion relation is

(5.12) $$\begin{eqnarray}{\it\omega}^{2}=(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{a})^{2}.\end{eqnarray}$$

Given the dispersion relation ${\it\omega}(\boldsymbol{k})$ of any wave mode, the phase and group velocities of the wave can be identified as

(5.13a,b ) $$\begin{eqnarray}\boldsymbol{v}_{p}=\frac{{\it\omega}}{k}\,\hat{\boldsymbol{k}},\quad \boldsymbol{v}_{g}=\frac{\partial {\it\omega}}{\partial \boldsymbol{k}}=\boldsymbol{{\rm\nabla}}_{\!\boldsymbol{k}}{\it\omega},\end{eqnarray}$$

where $\hat{\boldsymbol{k}}=\boldsymbol{k}/k$ . The phase velocity is that with which the phase of the wave travels, while the group velocity is that which the energy of the wave (or the centre of a wavepacket) is transported.

For Alfvén waves, therefore,

(5.14a,b ) $$\begin{eqnarray}\boldsymbol{v}_{p}=\pm v_{a}\cos {\it\theta}\,\hat{\boldsymbol{k}},\quad \boldsymbol{v}_{g}=\pm \boldsymbol{v}_{a},\end{eqnarray}$$

where ${\it\theta}$ is the angle between $\boldsymbol{k}$ and $\boldsymbol{B}$ .

To find the other solutions, we take the inner product of (5.10) with $\boldsymbol{k}$ and then with $\boldsymbol{B}$ to obtain first

(5.15) $$\begin{eqnarray}{\it\rho}{\it\omega}^{2}\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi}=\left[\left({\it\gamma}p+\frac{B^{2}}{{\it\mu}_{0}}\right)\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi}-\frac{1}{{\it\mu}_{0}}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B})\boldsymbol{B}\boldsymbol{\cdot }{\bf\xi}\right]k^{2}\end{eqnarray}$$

and then

(5.16) $$\begin{eqnarray}{\it\rho}{\it\omega}^{2}\boldsymbol{B}\boldsymbol{\cdot }{\bf\xi}={\it\gamma}p(\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi})\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B}.\end{eqnarray}$$

These equations can be written in the form

(5.17) $$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}{\it\rho}{\it\omega}^{2}-\left({\it\gamma}p+{\displaystyle \frac{B^{2}}{{\it\mu}_{0}}}\right)k^{2} & {\displaystyle \frac{1}{{\it\mu}_{0}}}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B})k^{2}\\ -{\it\gamma}p(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B}) & {\it\rho}{\it\omega}^{2}\end{array}\right]\left[\begin{array}{@{}c@{}}\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi}\\ \boldsymbol{B}\boldsymbol{\cdot }{\bf\xi}\end{array}\right]=\left[\begin{array}{@{}c@{}}0\\ 0\end{array}\right].\end{eqnarray}$$

The ‘trivial solution’ $\boldsymbol{k}\boldsymbol{\cdot }{\bf\xi}=\boldsymbol{B}\boldsymbol{\cdot }{\bf\xi}=0$ corresponds to the Alfvén wave that we have already identified. The other solutions satisfy

(5.18) $$\begin{eqnarray}{\it\rho}{\it\omega}^{2}\left[{\it\rho}{\it\omega}^{2}-\left({\it\gamma}p+\frac{B^{2}}{{\it\mu}_{0}}\right)k^{2}\right]+{\it\gamma}pk^{2}\frac{1}{{\it\mu}_{0}}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{B})^{2}=0,\end{eqnarray}$$

which simplifies to

(5.19) $$\begin{eqnarray}v_{p}^{4}-(v_{s}^{2}+v_{a}^{2})v_{p}^{2}+v_{s}^{2}v_{a}^{2}\cos ^{2}{\it\theta}=0.\end{eqnarray}$$

The two solutions

(5.20) $$\begin{eqnarray}v_{p}^{2}={\textstyle \frac{1}{2}}(v_{s}^{2}+v_{a}^{2})\pm \left[{\textstyle \frac{1}{4}}(v_{s}^{2}+v_{a}^{2})^{2}-v_{s}^{2}v_{a}^{2}\cos ^{2}{\it\theta}\right]^{1/2}\end{eqnarray}$$

are called fast and slow magnetoacoustic (or magnetosonic) waves, respectively.

In the special case ${\it\theta}=0$ ( $\boldsymbol{k}\Vert \boldsymbol{B}$ ), we have

(5.21) $$\begin{eqnarray}v_{p}^{2}=v_{s}^{2}\quad \text{or}\quad v_{a}^{2},\end{eqnarray}$$

together with $v_{p}^{2}=v_{a}^{2}$ for the Alfvén wave. Note that the fast wave could be either $v_{p}^{2}=v_{s}^{2}$ or $v_{p}^{2}=v_{a}^{2}$ , whichever is greater.

Figure 4. Polar plots of the phase velocity (a,c) and group velocity (b,d) of MHD waves for the cases $v_{a}=0.7\,v_{s}$ (a,b) and $v_{s}=0.7\,v_{a}$ (c,d) with a magnetic field in the horizontal direction. (The group velocity plot for the Alfvén wave consists of the two points $(\pm v_{a},0)$ .)

In the special case ${\it\theta}={\rm\pi}/2$ ( $\boldsymbol{k}\bot \boldsymbol{B}$ ), we have

(5.22) $$\begin{eqnarray}v_{p}^{2}=v_{s}^{2}+v_{a}^{2}\quad \text{or}\quad 0,\end{eqnarray}$$

together with $v_{p}^{2}=0$ for the Alfvén wave.

The effects of the magnetic field on wave propagation can be understood as resulting from the two aspects of the Lorentz force. The magnetic tension gives rise to Alfvén waves, which are similar to waves on an elastic string, and are trivial in the absence of the magnetic field. In addition, the magnetic pressure affects the response of the fluid to compression, and therefore modifies the propagation of acoustic waves.

The phase and group velocity for the full range of ${\it\theta}$ are usually exhibited in Friedrichs diagrams ⁶ (figure 4).

We can interpret:

1. (i) the fast wave as a quasi-isotropic acoustic-type wave in which both gas and magnetic pressure contribute;

2. (ii) the slow wave as an acoustic-type wave that is strongly guided by the magnetic field;

3. (iii) the Alfvén wave as analogous to a wave on an elastic string, propagating by means of magnetic tension and perfectly guided by the magnetic field.

Related example: A.12.

6 Nonlinear waves, shocks and other discontinuities

6.1 One-dimensional gas dynamics

6.1.1 Riemann’s analysis

The equations of mass conservation and motion in one dimension are

(6.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial {\it\rho}}{\partial t}+u\frac{\partial {\it\rho}}{\partial x}=-{\it\rho}\frac{\partial u}{\partial x},\\ \displaystyle \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial x}.\end{array}\right\}\end{eqnarray}$$

We assume the gas is homentropic ( $s=\text{const.}$ ) and perfect. (This eliminates the entropy wave and leaves only the two sound waves.) Then $p\propto {\it\rho}^{{\it\gamma}}$ and $v_{s}^{2}={\it\gamma}p/{\it\rho}\propto {\it\rho}^{{\it\gamma}-1}$ . It is convenient to use $v_{s}$ as a variable in place of ${\it\rho}$ or $p$ :

(6.2a,b ) $$\begin{eqnarray}\text{d}p=v_{s}^{2}\,\text{d}{\it\rho},\quad \text{d}{\it\rho}=\frac{{\it\rho}}{v_{s}}\left(\frac{2\,\text{d}v_{s}}{{\it\gamma}-1}\right).\end{eqnarray}$$

Then

(6.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v_{s}\frac{\partial }{\partial x}\left(\frac{2v_{s}}{{\it\gamma}-1}\right)=0,\\ \displaystyle \frac{\partial }{\partial t}\left(\frac{2v_{s}}{{\it\gamma}-1}\right)+u\frac{\partial }{\partial x}\left(\frac{2v_{s}}{{\it\gamma}-1}\right)+v_{s}\frac{\partial u}{\partial x}=0.\end{array}\right\}\end{eqnarray}$$

We add and subtract to obtain

(6.4) $$\begin{eqnarray}\left[\frac{\partial }{\partial t}+(u+v_{s})\frac{\partial }{\partial x}\right]\left(u+\frac{2v_{s}}{{\it\gamma}-1}\right)=0,\end{eqnarray}$$

(6.5) $$\begin{eqnarray}\left[\frac{\partial }{\partial t}+(u-v_{s})\frac{\partial }{\partial x}\right]\left(u-\frac{2v_{s}}{{\it\gamma}-1}\right)=0.\end{eqnarray}$$

Define the two Riemann invariants

(6.6) $$\begin{eqnarray}R_{\pm }=u\pm \frac{2v_{s}}{{\it\gamma}-1}.\end{eqnarray}$$

Then we deduce that $R_{\pm }=\text{const.}$ along a characteristic (curve) of gradient $\text{d}x/\text{d}t=u\pm v_{s}$ in the $(x,t)$ plane. The $+$ and $-$ characteristics form an interlocking web covering the space–time diagram (figure 5).

Figure 5. Characteristic curves in the space–time diagram.

Note that both Riemann invariants are needed to reconstruct the solution ( $u$ and $v_{s}$ ). Half of the information is propagated along one set of characteristics and half along the other.

In general the characteristics are not known in advance but must be determined along with the solution. The $+$ and $-$ characteristics propagate at the speed of sound to the right and left, respectively, with respect to the motion of the fluid.

This concept generalizes to nonlinear waves the solution of the classical wave equation for acoustic waves on a uniform and static background, which is of the form $f(x-v_{s}t)+g(x+v_{s}t)$ .

6.1.2 Method of characteristics

A numerical method of solution can be based on the following idea:

1. (i) Start with the initial data ( $u$ and $v_{s}$ ) for all relevant $x$ at $t=0$ .

2. (ii) Determine the characteristic slopes at $t=0$ .

3. (iii) Propagate the $R_{\pm }$ information for a small increment of time, neglecting the variation of the characteristic slopes.

4. (iv) Combine the $R_{\pm }$ information to find $u$ and $v_{s}$ at each $x$ at the new value of $t$ .

5. (v) Re-evaluate the slopes and repeat.

The domain of dependence of a point $P$ in the space–time diagram is that region of the diagram bounded by the $\pm$ characteristics through $P$ and located in the past of $P$ . The solution at $P$ cannot depend on anything that occurs outside the domain of dependence. Similarly, the domain of influence of $P$ is the region in the future of $P$ bounded by the characteristics through $P$ (figure 6).

Figure 6. Domains of dependence and of influence.

6.1.3 A simple wave

Suppose that $R_{-}$ is uniform, having the same constant value on every characteristic emanating from an undisturbed region to the right. Its value everywhere is that of the undisturbed region:

(6.7) $$\begin{eqnarray}u-\frac{2v_{s}}{{\it\gamma}-1}=u_{0}-\frac{2v_{s0}}{{\it\gamma}-1}.\end{eqnarray}$$

Then, along the $+$ characteristics, both $R_{+}$ and $R_{-}$ , and therefore $u$ and $v_{s}$ , must be constant. The $+$ characteristics therefore have constant slope $v=u+v_{s}$ , so they are straight lines.

The statement that the wave speed $v$ is constant on the family of straight lines $\text{d}x/\text{d}t=v$ is expressed by the equation

(6.8) $$\begin{eqnarray}\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}=0.\end{eqnarray}$$

This is known as the inviscid Burgers equation ⁷ or the nonlinear advection equation.

The inviscid Burgers equation has only one set of characteristics, with slope $\text{d}x/\text{d}t=v$ . It is easily solved by the method of characteristics. The initial data define $v_{0}(x)=v(x,0)$ and the characteristics are straight lines. In regions where $\text{d}v_{0}/\text{d}x>0$ the characteristics diverge in the future. In regions where $\text{d}v_{0}/\text{d}x<0$ the characteristics converge and will form a shock at some point. Contradictory information arrives at the same point in the space–time diagram, leading to a breakdown of the solution (figure 7).

Figure 7. Formation of a shock from intersecting characteristics.

Another viewpoint is that of wave steepening. The graph of $v$ versus $x$ evolves in time by moving each point at its wave speed $v$ . The crest of the wave moves fastest and eventually overtakes the trough to the right of it. The profile would become multiple valued, but this is physically meaningless and the wave breaks, forming a discontinuity (figure 8).

Figure 8. Wave steepening and shock formation. The dotted profile is multiple valued and is replaced in practice with a discontinuous profile including a shock.

Indeed, the formal solution of the inviscid Burgers equation is

(6.9) $$\begin{eqnarray}v(x,t)=v_{0}(x_{0})\quad \text{with }x=x_{0}+v_{0}(x_{0})t.\end{eqnarray}$$

By the chain rule, $\partial v/\partial x=v_{0}^{\prime }/(1+v_{0}^{\prime }t)$ , which diverges first at the breaking time $t=1/\max (-v_{0}^{\prime })$ .

The crest of a sound wave moves faster than the trough for two reasons. It is partly because the crest is denser and hotter, so the sound speed is higher (unless the gas is isothermal), but it is also because of the self-advection of the wave (recall that the wave speed is $u+v_{s}$ ). The breaking time depends on the amplitude and wavelength of the wave.

6.2 General analysis of simple nonlinear waves

Recall the hyperbolic structure of the equations of AFD (§ 4.5):

(6.10) $$\begin{eqnarray}\frac{\partial \boldsymbol{U}}{\partial t}+\unicode[STIX]{x1D63C}_{i}\frac{\partial \boldsymbol{U}}{\partial x_{i}}=\mathbf{0},\quad \boldsymbol{U}=[{\it\rho},p,\boldsymbol{u},\boldsymbol{B}]^{\text{T}}.\end{eqnarray}$$

The system is hyperbolic because the eigenvalues of $\unicode[STIX]{x1D63C}_{i}n_{i}$ are real for any unit vector $n_{i}$ . The eigenvalues are identified as wave speeds, and the corresponding eigenvectors as wave modes.

In a simple wave propagating in the $x$ -direction, all physical quantities are functions of a single variable, the phase ${\it\varphi}(x,t)$ . Then $\boldsymbol{U}=\boldsymbol{U}({\it\varphi})$ and so

(6.11) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{U}}{\text{d}{\it\varphi}}\frac{\partial {\it\varphi}}{\partial t}+\unicode[STIX]{x1D63C}_{x}\frac{\text{d}\boldsymbol{U}}{\text{d}{\it\varphi}}\frac{\partial {\it\varphi}}{\partial x}=0.\end{eqnarray}$$

This equation is satisfied if $\text{d}\boldsymbol{U}/\text{d}{\it\varphi}$ is an eigenvector of the hyperbolic system and if

(6.12) $$\begin{eqnarray}\frac{\partial {\it\varphi}}{\partial t}+v\frac{\partial {\it\varphi}}{\partial x}=0,\end{eqnarray}$$

where $v$ is the corresponding wave speed. But since $v=v({\it\varphi})$ we again find

(6.13) $$\begin{eqnarray}\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}=0,\end{eqnarray}$$

the inviscid Burgers equation.

Wave steepening is therefore generic for simple waves. However, waves do not always steepen in practice. For example, linear dispersion arising from Coriolis or buoyancy forces (see § 11) can counteract nonlinear wave steepening. Waves propagating on a non-uniform background are not simple waves. In addition, waves may be damped by diffusive processes (viscosity, thermal conduction or resistivity) before they can steepen.

Furthermore, even some simple waves do not undergo steepening, in spite of the above argument. This happens if the wave speed $v$ does not depend on the variables that actually vary in the wave mode. One example is the entropy wave in hydrodynamics, in which the density varies but not the pressure or the velocity. The wave speed is the fluid velocity, which does not vary in this wave; therefore the relevant solution of the inviscid Burgers equation is just $v=\text{const}$ . Another example is the Alfvén wave, which involves variations in transverse velocity and magnetic field components, but whose speed depends on the longitudinal components and the density. The slow and fast magnetoacoustic waves, though, are ‘genuinely nonlinear’ and undergo steepening.

6.3 Shocks and other discontinuities

6.3.1 Jump conditions

Discontinuities are resolved in reality by diffusive processes (viscosity, thermal conduction or resistivity) that become more important on smaller length scales. Properly, we should solve an enhanced set of equations to resolve the internal structure of a shock. This internal solution would then be matched on to the external solution in which diffusion can be neglected. However, the matching conditions can in fact be determined from general principles without resolving the internal structure.

Without loss of generality, we consider a shock front at rest at $x=0$ (making a Galilean transformation if necessary). We look for a stationary, one-dimensional solution in which gas flows from left to right. On the left is upstream, pre-shock material ( ${\it\rho}_{1}$ , $p_{1}$ , etc.). On the right is downstream, post-shock material ( ${\it\rho}_{2}$ , $p_{2}$ , etc.) (figure 9).

Figure 9. A shock front in its rest frame.

Consider any equation in conservative form

(6.14) $$\begin{eqnarray}\frac{\partial q}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{F}=0.\end{eqnarray}$$

For a stationary solution in one dimension,

(6.15) $$\begin{eqnarray}\frac{\text{d}F_{x}}{\text{d}x}=0,\end{eqnarray}$$

which implies that the flux density $F_{x}$ has the same value on each side of the shock. We write the matching condition as

(6.16) $$\begin{eqnarray}[F_{x}]_{1}^{2}=F_{x2}-F_{x1}=0.\end{eqnarray}$$

Including additional physics means that additional diffusive fluxes (not of mass but of momentum, energy, magnetic flux, etc.) are present. These fluxes are negligible outside the shock, so they do not affect the jump conditions. This approach is permissible provided that the new physics does not introduce any source terms in the equations. So the total energy is a properly conserved quantity but not the entropy (see later).

From mass conservation:

(6.17) $$\begin{eqnarray}[{\it\rho}u_{x}]_{1}^{2}=0.\end{eqnarray}$$

From momentum conservation:

(6.18) $$\begin{eqnarray}\left[{\it\rho}u_{x}^{2}+{\it\Pi}-\frac{B_{x}^{2}}{{\it\mu}_{0}}\right]_{1}^{2}=0,\end{eqnarray}$$

(6.19) $$\begin{eqnarray}\left[{\it\rho}u_{x}u_{y}-\frac{B_{x}B_{y}}{{\it\mu}_{0}}\right]_{1}^{2}=0,\end{eqnarray}$$

(6.20) $$\begin{eqnarray}\left[{\it\rho}u_{x}u_{z}-\frac{B_{x}B_{z}}{{\it\mu}_{0}}\right]_{1}^{2}=0.\end{eqnarray}$$

From $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$ :

(6.21) $$\begin{eqnarray}[B_{x}]_{1}^{2}=0.\end{eqnarray}$$

From $\partial \boldsymbol{B}/\partial t+\boldsymbol{{\rm\nabla}}\times \boldsymbol{E}=\mathbf{0}$ :

(6.22) $$\begin{eqnarray}[u_{x}B_{y}-u_{y}B_{x}]_{1}^{2}=-[E_{z}]_{1}^{2}=0,\end{eqnarray}$$

(6.23) $$\begin{eqnarray}[u_{x}B_{z}-u_{z}B_{x}]_{1}^{2}=[E_{y}]_{1}^{2}=0.\end{eqnarray}$$

(These are the standard electromagnetic conditions at an interface: the normal component of $\boldsymbol{B}$ and the tangential components of $\boldsymbol{E}$ are continuous.) From total energy conservation:

(6.24) $$\begin{eqnarray}\left[{\it\rho}u_{x}\left(\frac{1}{2}u^{2}+h\right)+\frac{1}{{\it\mu}_{0}}(E_{y}B_{z}-E_{z}B_{y})\right]_{1}^{2}=0.\end{eqnarray}$$

Note that the conservative form of the momentum equation used above is (cf. 4.45)

(6.25) $$\begin{eqnarray}\frac{\partial }{\partial t}({\it\rho}u_{i})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left({\it\rho}u_{i}\boldsymbol{u}+{\it\Pi}\,\boldsymbol{e}_{i}-\frac{B_{i}\boldsymbol{B}}{{\it\mu}_{0}}\right)=0.\end{eqnarray}$$

Including gravity makes no difference to the shock relations because ${\it\Phi}$ is always continuous (it satisfies ${\rm\nabla}^{2}{\it\Phi}=4{\rm\pi}G{\it\rho}$ ).

Although the entropy in ideal MHD satisfies an equation of conservative form,

(6.26) $$\begin{eqnarray}\frac{\partial }{\partial t}({\it\rho}s)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}s\boldsymbol{u})=0,\end{eqnarray}$$

the dissipation of energy within the shock provides a source term for entropy. Therefore the entropy flux is not continuous across the shock.

6.3.2 Non-magnetic shocks

First consider a normal shock ( $u_{y}=u_{z}=0$ ) with no magnetic field. We obtain the Rankine–Hugionot relations ⁸

(6.27a-c ) $$\begin{eqnarray}[{\it\rho}u_{x}]_{1}^{2}=0,\quad [{\it\rho}u_{x}^{2}+p]_{1}^{2}=0,\quad \left[{\it\rho}u_{x}\left({\textstyle \frac{1}{2}}u_{x}^{2}+h\right)\right]_{1}^{2}=0.\end{eqnarray}$$

The specific enthalpy of a perfect gas is

(6.28) $$\begin{eqnarray}h=\left(\frac{{\it\gamma}}{{\it\gamma}-1}\right)\frac{p}{{\it\rho}}\end{eqnarray}$$

and these equations can be solved algebraically (see Example A.13). Introduce the upstream Mach number (the shock Mach number)

(6.29) $$\begin{eqnarray}{\mathcal{M}}_{1}=\frac{u_{x1}}{v_{\text{s}1}}>0.\end{eqnarray}$$

Then we find

(6.30a,b ) $$\begin{eqnarray}\frac{{\it\rho}_{2}}{{\it\rho}_{1}}=\frac{u_{x1}}{u_{x2}}=\frac{({\it\gamma}+1){\mathcal{M}}_{1}^{2}}{({\it\gamma}-1){\mathcal{M}}_{1}^{2}+2},\quad \frac{p_{2}}{p_{1}}=\frac{2{\it\gamma}{\mathcal{M}}_{1}^{2}-({\it\gamma}-1)}{({\it\gamma}+1)},\end{eqnarray}$$

and

(6.31) $$\begin{eqnarray}{\mathcal{M}}_{2}^{2}=\frac{2+({\it\gamma}-1){\mathcal{M}}_{1}^{2}}{2{\it\gamma}{\mathcal{M}}_{1}^{2}-({\it\gamma}-1)}.\end{eqnarray}$$

Note that ${\it\rho}_{2}/{\it\rho}_{1}$ and $p_{2}/p_{1}$ are increasing functions of ${\mathcal{M}}_{1}$ . The case ${\mathcal{M}}_{1}=1$ is trivial as it corresponds to ${\it\rho}_{2}/{\it\rho}_{1}=p_{2}/p_{1}=1$ . The other two cases are the compression shock ( ${\mathcal{M}}_{1}>1$ , ${\mathcal{M}}_{2}<1$ , ${\it\rho}_{2}>{\it\rho}_{1}$ , $p_{2}>p_{1}$ ) and the rarefaction shock ( ${\mathcal{M}}_{1}<1$ , ${\mathcal{M}}_{2}>1$ , ${\it\rho}_{2}<{\it\rho}_{1}$ , $p_{2}<p_{1}$ ).

It is shown in Example A.13 that the entropy change in passing through the shock is positive for compression shocks and negative for rarefaction shocks. Therefore only compression shocks are physically realizable. Rarefaction shocks are excluded by the second law of thermodynamics. All shocks involve dissipation and irreversibility.

The fact that ${\mathcal{M}}_{1}>1$ while ${\mathcal{M}}_{2}<1$ means that the shock travels supersonically relative to the upstream gas and subsonically relative to the downstream gas.

In the weak shock limit ${\mathcal{M}}_{1}-1\ll 1$ the relative velocity of the fluid and the shock is close to the sound speed on both sides.

In the strong shock limit ${\mathcal{M}}_{1}\gg 1$ , common in astrophysical applications, we have

(6.32a-c ) $$\begin{eqnarray}\frac{{\it\rho}_{2}}{{\it\rho}_{1}}=\frac{u_{x1}}{u_{x2}}\rightarrow \frac{{\it\gamma}+1}{{\it\gamma}-1},\quad \frac{p_{2}}{p_{1}}\gg 1,\quad {\mathcal{M}}_{2}^{2}\rightarrow \frac{{\it\gamma}-1}{2{\it\gamma}}.\end{eqnarray}$$

Note that the compression ratio ${\it\rho}_{2}/{\it\rho}_{1}$ is finite (and equal to $4$ when ${\it\gamma}=5/3$ ). In the rest frame of the undisturbed (upstream) gas the shock speed is $u_{sh}=-u_{x1}$ . The downstream density, velocity (in that frame) and pressure in the limit of a strong shock are (as will be used in § 7)

(6.33a-c ) $$\begin{eqnarray}{\it\rho}_{2}=\left(\frac{{\it\gamma}+1}{{\it\gamma}-1}\right){\it\rho}_{1},\quad u_{x2}-u_{x1}=\frac{2u_{sh}}{{\it\gamma}+1},\quad p_{2}=\frac{2{\it\rho}_{1}u_{sh}^{2}}{{\it\gamma}+1}.\end{eqnarray}$$

A significant amount of thermal energy is generated out of kinetic energy by the passage of a strong shock:

(6.34) $$\begin{eqnarray}e_{2}=\frac{2u_{sh}^{2}}{({\it\gamma}+1)^{2}}.\end{eqnarray}$$

6.3.3 Oblique shocks

When $u_{y}$ or $u_{z}$ is non-zero, we have the additional relations

(6.35) $$\begin{eqnarray}[{\it\rho}u_{x}u_{y}]_{1}^{2}=[{\it\rho}u_{x}u_{z}]_{1}^{2}=0.\end{eqnarray}$$

Since ${\it\rho}u_{x}$ is continuous across the shock (and non-zero), we deduce that $[u_{y}]_{1}^{2}=[u_{z}]_{1}^{2}=0$ . Momentum and energy conservation apply as before, and we recover the Rankine–Hugoniot relations. (See Example A.14.)

6.3.4 Other discontinuities

The discontinuity is not called a shock if there is no normal flow ( $u_{x}=0$ ). In this case we can deduce only that $[p]_{1}^{2}=0$ . Arbitrary discontinuities are allowed in ${\it\rho}$ , $u_{y}$ and $u_{z}$ . These are related to the entropy and vortical waves. If there is a jump in ${\it\rho}$ we have a contact discontinuity. If there is a jump in $u_{y}$ or $u_{z}$ we have a tangential discontinuity or vortex sheet (the vorticity being proportional to ${\it\delta}(x)$ ). Note that these discontinuities are not produced naturally by wave steepening, because the entropy and vortical waves do not steepen. However they do appear in the Riemann problem (§ 6.3.6) and other situations with discontinuous initial conditions.

6.3.5 MHD shocks and discontinuities

When a magnetic field is included, the jump conditions allow a wider variety of solutions. There are different types of discontinuity associated with the three MHD waves (Alfvén, slow and fast), which we will not discuss here. Since the parallel components of $\boldsymbol{B}$ need not be continuous, it is possible for them to ‘switch on’ or ‘switch off’ on passage through a shock.

A current sheet is a tangential discontinuity in the magnetic field. A classic case would be where $B_{y}$ , say, changes sign across the interface, with $B_{x}=0$ . The current density is then proportional to ${\it\delta}(x)$ .

6.3.6 The Riemann problem

The Riemann problem is a fundamental initial value problem for a hyperbolic system and plays a central role in some numerical methods for solving the equations of AFD.

The initial condition at $t=0$ consists of two uniform states separated by a discontinuity at $x=0$ . In the case of one-dimensional gas dynamics, we have

(6.36a-c ) $$\begin{eqnarray}{\it\rho}=\left\{\begin{array}{@{}ll@{}}{\it\rho}_{L},\quad & x<0\\ {\it\rho}_{R},\quad & x>0\end{array}\right.,\quad p=\left\{\begin{array}{@{}ll@{}}p_{L},\quad & x<0\\ p_{R},\quad & x>0\end{array}\right.,\quad u_{x}=\left\{\begin{array}{@{}ll@{}}u_{L},\quad & x<0\\ u_{R},\quad & x>0\end{array}\right.,\end{eqnarray}$$

where ‘L’ and ‘R’ denote the left and right states. A simple example is a ‘shock-tube’ problem in which gas at different pressures is at rest on either side of a partition, which is released at $t=0$ .

It can be shown that the initial discontinuity resolves generically into three simple waves. The inner one is a contact discontinuity while the outer ones are shocks or rarefaction waves (see below).

The initial data define no natural length scale for the Riemann problem, but they do allow a characteristic velocity scale $c$ to be defined (although not uniquely). The result is a similarity solution in which variables depend on $x$ and $t$ only through the dimensionless combination ${\it\xi}=x/ct$ .

Unlike the unphysical rarefaction shock, the rarefaction wave (or expansion wave) is a non-dissipative, homentropic, continuous simple wave in which $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}>0$ . If we seek a similarity solution $v=v({\it\xi})$ of the inviscid Burgers equation $v_{t}+vv_{x}=0$ we find $v=x/t$ (or the trivial solution $v=\text{const.}$ ). The characteristics form an expansion fan (figure 10).

Figure 10. Expansion fan of characteristics in a rarefaction wave.

The ‘ $+$ ’ rarefaction wave has $u+v_{s}=x/t$ and $R_{-}=u-2v_{s}/({\it\gamma}-1)=\text{const.}$ , determined by the undisturbed right-hand state. The ‘ $-$ ’ rarefaction wave has $u-v_{s}=x/t$ and $R_{+}=u+2v_{s}/({\it\gamma}-1)=\text{const.}$ , determined by the undisturbed left-hand state. In each case $u$ and $v_{s}$ are linear functions of $x/t$ and $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=2t^{-1}/({\it\gamma}+1)>0$ .

A typical outcome of a shock-tube problem consists of (from left to right): undisturbed region, rarefaction wave, uniform region, contact discontinuity, uniform region, shock, undisturbed region (figure 11).

Figure 11. Typical outcome of a shock-tube problem. The two uniform regions are separated by a contact discontinuity. The other discontinuity is a shock.

In Godunov’s method and related computational algorithms, the equations of AFD are advanced in time by solving (either exactly or approximately) a Riemann problem at each cell boundary.

Related examples: A.13–A.16.

7 Spherical blast waves: supernovae

Note: in this section $(r,{\it\theta},{\it\phi})$ are spherical polar coordinates.

7.1 Introduction

In a supernova, an energy of order $10^{51}~\text{erg}$ ( $10^{44}~\text{J}$ ) is released into the interstellar medium. An expanding spherical blast wave is formed as the explosion sweeps up the surrounding gas. Several good examples of these supernova remnants are observed in the Galaxy, e.g. Tycho’s supernova of 1572 and Kepler’s supernova of 1604 ⁹ .

The effect is similar to a bomb. When photographs ¹⁰ (complete with length and time scales) were released of the first atomic bomb test in New Mexico in 1945, both Sedov ¹¹ in the Soviet Union and Taylor ¹² in the UK were able to work out the energy of the bomb (equivalent to approximately 20 kilotons of TNT), which was supposed to be a secret.

We suppose that an energy $E$ is released at $t=0$ , $r=0$ and that the explosion is spherically symmetric. The external medium has density ${\it\rho}_{0}$ and pressure $p_{0}$ . In the Sedov–Taylor phase of the explosion, the pressure $p\gg p_{0}$ . Then a strong shock is formed and the external pressure $p_{0}$ can be neglected (formally set to zero). Gravity is also negligible in the dynamics.

7.2 Governing equations

The equations governing the spherically symmetric flow of a perfect gas, with radial velocity $u_{r}=u(r,t)$ , may be written as

(7.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \left(\frac{\partial }{\partial t}+u\frac{\partial }{\partial r}\right){\it\rho}=-\frac{{\it\rho}}{r^{2}}\frac{\partial }{\partial r}(r^{2}u),\\ \displaystyle \left(\frac{\partial }{\partial t}+u\frac{\partial }{\partial r}\right)u=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial r},\\ \displaystyle \left(\frac{\partial }{\partial t}+u\frac{\partial }{\partial r}\right)\ln (p{\it\rho}^{-{\it\gamma}})=0.\end{array}\right\}\end{eqnarray}$$

These imply the total energy equation

(7.2) $$\begin{eqnarray}\frac{\partial }{\partial t}\left(\frac{1}{2}{\it\rho}u^{2}+\frac{p}{{\it\gamma}-1}\right)+\frac{1}{r^{2}}\frac{\partial }{\partial r}\left[r^{2}\left(\frac{1}{2}{\it\rho}u^{2}+\frac{{\it\gamma}p}{{\it\gamma}-1}\right)u\right]=0.\end{eqnarray}$$

The shock is at $r=R(t)$ , and the shock speed is ${\dot{R}}$ . The equations are to be solved in $0<r<R$ with the strong shock conditions (6.33) at $r=R$ :

(7.3a-c ) $$\begin{eqnarray}{\it\rho}=\left(\frac{{\it\gamma}+1}{{\it\gamma}-1}\right){\it\rho}_{0},\quad u=\frac{2{\dot{R}}}{{\it\gamma}+1},\quad p=\frac{2{\it\rho}_{0}{\dot{R}}^{2}}{{\it\gamma}+1}.\end{eqnarray}$$

The total energy of the explosion is

(7.4) $$\begin{eqnarray}E=\int _{0}^{R}\left(\frac{1}{2}{\it\rho}u^{2}+\frac{p}{{\it\gamma}-1}\right)\,4{\rm\pi}r^{2}\,\text{d}r=\text{const.},\end{eqnarray}$$

the thermal energy of the external medium being negligible.

7.3 Dimensional analysis

The dimensional parameters of the problem on which the solution might depend are $E$ and ${\it\rho}_{0}$ . Their dimensions are

(7.5a,b ) $$\begin{eqnarray}[E]=ML^{2}T^{-2},\quad [{\it\rho}_{0}]=ML^{-3}.\end{eqnarray}$$

Together, they do not define a characteristic length scale, so the explosion is ‘scale free’ or ‘self-similar’. If the dimensional analysis includes the time $t$ since the explosion, however, we can find a time-dependent characteristic length scale. The radius of the shock must be

(7.6) $$\begin{eqnarray}R={\it\alpha}\left(\frac{Et^{2}}{{\it\rho}_{0}}\right)^{1/5},\end{eqnarray}$$

where ${\it\alpha}$ is a dimensionless constant to be determined.

7.4 Similarity solution

The self-similarity of the explosion is expressed using the dimensionless similarity variable ${\it\xi}=r/R(t)$ . The solution has the form

(7.7a-c ) $$\begin{eqnarray}{\it\rho}={\it\rho}_{0}\,\tilde{{\it\rho}}({\it\xi}),\quad u={\dot{R}}\,\tilde{u} ({\it\xi}),\quad p={\it\rho}_{0}{\dot{R}}^{2}\,\tilde{p}({\it\xi}),\end{eqnarray}$$

where $\tilde{{\it\rho}}({\it\xi})$ , $\tilde{u} ({\it\xi})$ and $\tilde{p}({\it\xi})$ are dimensionless functions to be determined. The meaning of this type of solution is that the graph of $u$ versus $r$ , for example, has a constant shape but both axes of the graph are rescaled as time proceeds and the shock expands.

7.5 Dimensionless equations

We substitute these forms into (7.1) and cancel the dimensional factors to obtain

(7.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (\tilde{u} -{\it\xi})\tilde{{\it\rho}}^{\prime }=-\frac{\tilde{{\it\rho}}}{{\it\xi}^{2}}\frac{\text{d}}{\text{d}{\it\xi}}({\it\xi}^{2}\tilde{u} ),\\ \displaystyle (\tilde{u} -{\it\xi})\tilde{u} ^{\prime }-\frac{3}{2}\tilde{u} =-\frac{\tilde{p}^{\prime }}{\tilde{{\it\rho}}},\\ \displaystyle (\tilde{u} -{\it\xi})\left(\frac{\tilde{p}^{\prime }}{\tilde{p}}-\frac{{\it\gamma}\tilde{{\it\rho}}^{\prime }}{\tilde{{\it\rho}}}\right)-3=0.\end{array}\right\}\end{eqnarray}$$

Similarly, the strong shock conditions are that

(7.9a-c ) $$\begin{eqnarray}\tilde{{\it\rho}}=\frac{{\it\gamma}+1}{{\it\gamma}-1},\quad \tilde{u} =\frac{2}{{\it\gamma}+1},\quad \tilde{p}=\frac{2}{{\it\gamma}+1}\end{eqnarray}$$

at ${\it\xi}=1$ , while the total energy integral provides a normalization condition,

(7.10) $$\begin{eqnarray}1=\frac{16{\rm\pi}}{25}{\it\alpha}^{5}\int _{0}^{1}\left(\frac{1}{2}\tilde{{\it\rho}}\tilde{u} ^{2}+\frac{\tilde{p}}{{\it\gamma}-1}\right)\,{\it\xi}^{2}\,\text{d}{\it\xi},\end{eqnarray}$$

that will ultimately determine the value of ${\it\alpha}$ .

7.6 First integral

The one-dimensional conservative form of the total energy equation (7.2) is

(7.11) $$\begin{eqnarray}\frac{\partial q}{\partial t}+\frac{\partial F}{\partial r}=0,\end{eqnarray}$$

where

(7.12a,b ) $$\begin{eqnarray}q=r^{2}\left(\frac{1}{2}{\it\rho}u^{2}+\frac{p}{{\it\gamma}-1}\right),\quad F=r^{2}\left(\frac{1}{2}{\it\rho}u^{2}+\frac{{\it\gamma}p}{{\it\gamma}-1}\right)u.\end{eqnarray}$$

In dimensionless form,

(7.13a,b ) $$\begin{eqnarray}q={\it\rho}_{0}R^{2}{\dot{R}}^{2}\,\tilde{q}({\it\xi}),\quad F={\it\rho}_{0}R^{2}{\dot{R}}^{3}\,\tilde{F}({\it\xi}),\end{eqnarray}$$

with

(7.14a,b ) $$\begin{eqnarray}\tilde{q}={\it\xi}^{2}\left(\frac{1}{2}\tilde{{\it\rho}}\tilde{u} ^{2}+\frac{\tilde{p}}{{\it\gamma}-1}\right),\quad \tilde{F}={\it\xi}^{2}\left(\frac{1}{2}\tilde{{\it\rho}}\tilde{u} ^{2}+\frac{{\it\gamma}\tilde{p}}{{\it\gamma}-1}\right)\tilde{u} .\end{eqnarray}$$

We substitute the forms (7.13) into the energy equation (7.11) to find

(7.15) $$\begin{eqnarray}-{\it\xi}\tilde{q}^{\prime }-\tilde{q}+\tilde{F}^{\prime }=0,\end{eqnarray}$$

which implies

(7.16) $$\begin{eqnarray}\frac{\text{d}}{\text{d}{\it\xi}}(\tilde{F}-{\it\xi}\tilde{q})=0.\end{eqnarray}$$

Thus

(7.17) $$\begin{eqnarray}\tilde{F}-{\it\xi}\tilde{q}=\text{const.}=0\end{eqnarray}$$

for a solution that is finite at ${\it\xi}=0$ . This equation can be solved for $\tilde{p}$ :

(7.18) $$\begin{eqnarray}\tilde{p}=\frac{({\it\gamma}-1)\tilde{{\it\rho}}\tilde{u} ^{2}({\it\xi}-\tilde{u} )}{2({\it\gamma}\tilde{u} -{\it\xi})}.\end{eqnarray}$$

Note that this solution is compatible with the shock boundary conditions. Having found a first integral, we can now dispense with (e.g.) the thermal energy equation.

Let $\tilde{u} =v{\it\xi}$ . We now have

(7.19) $$\begin{eqnarray}(v-1)\frac{\text{d}\ln \tilde{{\it\rho}}}{\text{d}\ln {\it\xi}}=-\frac{\text{d}v}{\text{d}\ln {\it\xi}}-3v,\end{eqnarray}$$

(7.20) $$\begin{eqnarray}(v-1)\frac{\text{d}v}{\text{d}\ln {\it\xi}}+\frac{1}{\tilde{{\it\rho}}{\it\xi}^{2}}\frac{\text{d}}{\text{d}\ln {\it\xi}}\left[\frac{({\it\gamma}-1)\tilde{{\it\rho}}{\it\xi}^{2}v^{2}(1-v)}{2({\it\gamma}v-1)}\right]=\frac{3}{2}v.\end{eqnarray}$$

Eliminate $\tilde{{\it\rho}}$ :

(7.21) $$\begin{eqnarray}\frac{\text{d}v}{\text{d}\ln {\it\xi}}=\frac{v({\it\gamma}v-1)[5-(3{\it\gamma}-1)v]}{{\it\gamma}({\it\gamma}+1)v^{2}-2({\it\gamma}+1)v+2}.\end{eqnarray}$$

Invert and split into partial fractions:

(7.22) $$\begin{eqnarray}\frac{\text{d}\ln {\it\xi}}{\text{d}v}=-\frac{2}{5v}+\frac{{\it\gamma}({\it\gamma}-1)}{(2{\it\gamma}+1)({\it\gamma}v-1)}+\frac{13{\it\gamma}^{2}-7{\it\gamma}+12}{5(2{\it\gamma}+1)[5-(3{\it\gamma}-1)v]}.\end{eqnarray}$$

The solution is

(7.23) $$\begin{eqnarray}{\it\xi}\propto v^{-2/5}({\it\gamma}v-1)^{({\it\gamma}-1)/(2{\it\gamma}+1)}[5-(3{\it\gamma}-1)v]^{-(13{\it\gamma}^{2}-7{\it\gamma}+12)/5(2{\it\gamma}+1)(3{\it\gamma}-1)}.\end{eqnarray}$$

Now

(7.24) $$\begin{eqnarray}\displaystyle \frac{\text{d}\ln \tilde{{\it\rho}}}{\text{d}v} & = & \displaystyle -\frac{1}{v-1}-\frac{3v}{v-1}\frac{\text{d}\ln {\it\xi}}{\text{d}v}\nonumber\\ \displaystyle & = & \displaystyle \frac{2}{(2-{\it\gamma})(1-v)}+\frac{3{\it\gamma}}{(2{\it\gamma}+1)({\it\gamma}v-1)}-\frac{13{\it\gamma}^{2}-7{\it\gamma}+12}{(2-{\it\gamma})(2{\it\gamma}+1)[5-(3{\it\gamma}-1)v]}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The solution is

(7.25) $$\begin{eqnarray}\tilde{{\it\rho}}\propto (1-v)^{-2/(2-{\it\gamma})}({\it\gamma}v-1)^{3/(2{\it\gamma}+1)}[5-(3{\it\gamma}-1)v]^{(13{\it\gamma}^{2}-7{\it\gamma}+12)/(2-{\it\gamma})(2{\it\gamma}+1)(3{\it\gamma}-1)}.\end{eqnarray}$$

For example, in the case ${\it\gamma}=5/3$ :

(7.26) $$\begin{eqnarray}{\it\xi}\propto v^{-2/5}\left(\frac{5v}{3}-1\right)^{2/13}(5-4v)^{-82/195},\end{eqnarray}$$

(7.27) $$\begin{eqnarray}\tilde{{\it\rho}}\propto (1-v)^{-6}\left(\frac{5v}{3}-1\right)^{9/13}(5-4v)^{82/13}.\end{eqnarray}$$

To satisfy $v=2/({\it\gamma}+1)=3/4$ and $\tilde{{\it\rho}}=({\it\gamma}+1)/({\it\gamma}-1)=4$ at ${\it\xi}=1$ :

(7.28) $$\begin{eqnarray}{\it\xi}=\left(\frac{4v}{3}\right)^{-2/5}\left(\frac{20v}{3}-4\right)^{2/13}\left(\frac{5}{2}-2v\right)^{-82/195},\end{eqnarray}$$

(7.29) $$\begin{eqnarray}\tilde{{\it\rho}}=4\left(4-4v\right)^{-6}\left(\frac{20v}{3}-4\right)^{9/13}\left(\frac{5}{2}-2v\right)^{82/13}.\end{eqnarray}$$

Then, from the first integral,

(7.30) $$\begin{eqnarray}\tilde{p}=\frac{3}{4}\left(\frac{4v}{3}\right)^{6/5}\left(4-4v\right)^{-5}\left(\frac{5}{2}-2v\right)^{82/15}.\end{eqnarray}$$

In this solution (figure 12), ${\it\xi}$ ranges from $0$ to $1$ , and $v$ from $3/5$ to $3/4$ . The normalization integral (numerically) yields ${\it\alpha}\approx 1.152$ .

Figure 12. Sedov’s solution for a spherical blast wave in the case ${\it\gamma}=5/3$ .

7.7 Applications

Some rough estimates are as follows:

1. (i) For a supernova: $E\sim 10^{51}~\text{erg}$ , ${\it\rho}_{0}\sim 10^{-24}~\text{g}~\text{cm}^{-3}$ . Then $R\approx 6~\text{pc}$ and ${\dot{R}}\approx 2000~\text{km}~\text{s}^{-1}$ at $t=1000~\text{year}$ .

2. (ii) For the 1945 New Mexico explosion: $E\approx 8\times 10^{20}~\text{erg}$ , ${\it\rho}_{0}\approx 1.2\times 10^{-3}~\text{g}~\text{cm}^{-3}$ . Then $R\approx 100~\text{m}$ and ${\dot{R}}\approx 4~\text{km}~\text{s}^{-1}$ at $t=0.01~\text{s}$ .

The similarity method is useful in a very wide range of nonlinear problems. In this case it reduced partial differential equations to integrable ordinary differential equations.

Related example: A.17.

8 Spherically symmetric steady flows: stellar winds and accretion

Note: in this section $(r,{\it\theta},{\it\phi})$ are spherical polar coordinates.

8.1 Introduction

Many stars, including the Sun, lose mass through a stellar wind. The gas must be sufficiently hot to escape from the star’s gravitational field. Gravitating bodies can also accrete gas from the interstellar medium. The simplest models of these processes neglect the effects of rotation or magnetic fields and involve a steady, spherically symmetric flow.

8.2 Basic equations

We consider a purely radial flow, either away from or towards a body of mass $M$ . The gas is perfect and non-self-gravitating, so ${\it\Phi}=-GM/r$ . The fluid variables are functions of $r$ only, and the only velocity component is $u_{r}=u(r)$ .

Mass conservation for such a flow implies that the mass flux

(8.1) $$\begin{eqnarray}4{\rm\pi}r^{2}{\it\rho}u=-{\dot{M}}=\text{const.}\end{eqnarray}$$

If $u>0$ (a stellar wind), $-{\dot{M}}$ is the mass loss rate. If $u<0$ (an accretion flow), ${\dot{M}}$ is the mass accretion rate. We ignore the secular change in the mass $M$ , which would otherwise violate the steady nature of the flow.

The thermal energy equation (assuming $u\neq 0$ ) implies homentropic flow:

(8.2) $$\begin{eqnarray}p=K{\it\rho}^{{\it\gamma}},\quad K=\text{const.}\end{eqnarray}$$

The equation of motion has only one component:

(8.3) $$\begin{eqnarray}{\it\rho}u\frac{\text{d}u}{\text{d}r}=-{\it\rho}\frac{\text{d}{\it\Phi}}{\text{d}r}-\frac{\text{d}p}{\text{d}r}.\end{eqnarray}$$

Alternatively, we can use the integral form (Bernoulli’s equation):

(8.4) $$\begin{eqnarray}\frac{1}{2}u^{2}+{\it\Phi}+h=B=\text{const.},\quad h=\left(\frac{{\it\gamma}}{{\it\gamma}-1}\right)\frac{p}{{\it\rho}}=\frac{v_{s}^{2}}{{\it\gamma}-1}.\end{eqnarray}$$

In highly subsonic flow the $u^{2}/2$ term on the left-hand side of Bernoulli’s equation is negligible and the gas is quasi-hydrostatic. In highly supersonic flow the $h$ term is negligible and the flow is quasi-ballistic (freely falling). As discussed below, we are usually interested in transonic solutions that pass smoothly from subsonic to supersonic flow.

Our aim is to solve for $u(r)$ , and to determine ${\dot{M}}$ if possible. At what rate does a star lose mass through a wind, or a black hole accrete mass from the surrounding medium?

8.3 First treatment

We first use the differential form of the equation of motion. Rewrite the pressure gradient using the other two equations:

(8.5) $$\begin{eqnarray}-\frac{\text{d}p}{\text{d}r}=-p\frac{\text{d}\ln p}{\text{d}r}=-{\it\gamma}p\frac{\text{d}\ln {\it\rho}}{\text{d}r}={\it\rho}v_{s}^{2}\left(\frac{2}{r}+\frac{1}{u}\frac{\text{d}u}{\text{d}r}\right).\end{eqnarray}$$

Equation (8.3), multiplied by $u/{\it\rho}$ , becomes

(8.6) $$\begin{eqnarray}(u^{2}-v_{s}^{2})\frac{\text{d}u}{\text{d}r}=u\left(\frac{2v_{s}^{2}}{r}-\frac{\text{d}{\it\Phi}}{\text{d}r}\right).\end{eqnarray}$$

A critical point (sonic point) occurs at any radius $r=r_{s}$ where $|u|=v_{s}$ . For the flow to pass smoothly from subsonic to supersonic, the right-hand side must vanish at the sonic point:

(8.7) $$\begin{eqnarray}\frac{2v_{ss}^{2}}{r_{s}}-\frac{GM}{r_{s}^{2}}=0.\end{eqnarray}$$

Evaluate Bernoulli’s equation (8.4) at the sonic point:

(8.8) $$\begin{eqnarray}\left(\frac{1}{2}+\frac{1}{{\it\gamma}-1}\right)v_{ss}^{2}-\frac{GM}{r_{s}}=B.\end{eqnarray}$$

We deduce that

(8.9a,b ) $$\begin{eqnarray}v_{ss}^{2}=\frac{2({\it\gamma}-1)}{(5-3{\it\gamma})}B,\quad r_{s}=\frac{(5-3{\it\gamma})}{4({\it\gamma}-1)}\frac{GM}{B}.\end{eqnarray}$$

There is a unique transonic solution, which exists only for $1\leqslant {\it\gamma}<5/3$ . (The case ${\it\gamma}=1$ can be treated separately or by taking a limit.)

Now evaluate ${\dot{M}}$ at the sonic point:

(8.10) $$\begin{eqnarray}|{\dot{M}}|=4{\rm\pi}r_{s}^{2}{\it\rho}_{s}v_{ss}.\end{eqnarray}$$

8.4 Second treatment

We now use Bernoulli’s equation instead of the equation of motion.

Introduce the local Mach number ${\mathcal{M}}=|u|/v_{s}$ . Then

(8.11a,b ) $$\begin{eqnarray}4{\rm\pi}r^{2}{\it\rho}v_{s}{\mathcal{M}}=|{\dot{M}}|,\quad v_{s}^{2}={\it\gamma}K{\it\rho}^{{\it\gamma}-1}.\end{eqnarray}$$

Eliminate ${\it\rho}$ to obtain

(8.12) $$\begin{eqnarray}v_{s}=({\it\gamma}K)^{1/({\it\gamma}+1)}\left(\frac{|{\dot{M}}|}{4{\rm\pi}r^{2}{\mathcal{M}}}\right)^{({\it\gamma}-1)/({\it\gamma}+1)}.\end{eqnarray}$$

Bernoulli’s equation (8.4) is

(8.13) $$\begin{eqnarray}\frac{1}{2}v_{s}^{2}{\mathcal{M}}^{2}-\frac{GM}{r}+\frac{v_{s}^{2}}{{\it\gamma}-1}=B.\end{eqnarray}$$

Substitute for $v_{s}$ and separate the variables:

(8.14) $$\begin{eqnarray}\displaystyle & & \displaystyle ({\it\gamma}K)^{2/({\it\gamma}+1)}\left(\frac{|{\dot{M}}|}{4{\rm\pi}}\right)^{2({\it\gamma}-1)/({\it\gamma}+1)}\left[\frac{{\mathcal{M}}^{4/({\it\gamma}+1)}}{2}+\frac{{\mathcal{M}}^{-2({\it\gamma}-1)/({\it\gamma}+1)}}{{\it\gamma}-1}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad =Br^{4({\it\gamma}-1)/({\it\gamma}+1)}+GMr^{-(5-3{\it\gamma})/({\it\gamma}+1)}.\end{eqnarray}$$

This equation is of the form $f({\mathcal{M}})=g(r)$ . Assume that $1<{\it\gamma}<5/3$ and $B>0$ . (If $B<0$ then the flow cannot reach infinity.) Then each of $f$ and $g$ is the sum of a positive power and a negative power, with positive coefficients. $f({\mathcal{M}})$ has a minimum at ${\mathcal{M}}=1$ , while $g(r)$ has a minimum at

(8.15) $$\begin{eqnarray}r=\frac{(5-3{\it\gamma})}{4({\it\gamma}-1)}\frac{GM}{B},\end{eqnarray}$$

which is the sonic radius $r_{s}$ identified previously. A smooth passage through the sonic point is possible only if $|{\dot{M}}|$ has a special value, so that the minima of $f$ and $g$ are equal. If $|{\dot{M}}|$ is too large then the solution does not work for all $r$ . If it is too small then the solution remains subsonic (or supersonic) for all $r$ , which may not agree with the boundary conditions (figure 13).

Figure 13. Shapes of the functions $f({\mathcal{M}})$ and $g(r)$ for the case ${\it\gamma}=4/3$ . Only if ${\dot{M}}$ is equal to the critical value at which the minima of $f$ and $g$ coincide (solid line, left panel) does a smooth transonic solution exist.

The $(r,{\mathcal{M}})$ plane shows an $\mathit{X}$ -type critical point at $(r_{s},1)$ (figure 14).

Figure 14. Solution curves for a stellar wind or accretion flow in the case ${\it\gamma}=4/3$ , showing an $\mathit{X}$ -type critical point at the sonic radius and at Mach number ${\mathcal{M}}=1$ .

For $r\ll r_{s}$ the subsonic solution is close to a hydrostatic atmosphere. The supersonic solution is close to free fall.

For $r\gg r_{s}$ the subsonic solution approaches a uniform state ( $p=\text{const.}$ , ${\it\rho}=\text{const.}$ ). The supersonic solution is close to $u=\text{const.}$ (so ${\it\rho}\propto r^{-2}$ ).

8.5 Stellar wind

For a stellar wind the appropriate solution is subsonic (quasi-hydrostatic) at small $r$ and supersonic (coasting) at large $r$ . Parker (1958) ¹³ first presented this simplified model for the solar wind. The mass loss rate can be determined from the properties of the quasi-hydrostatic part of the solution, e.g. the density and temperature at the base of the solar corona. A completely hydrostatic solution is unacceptable unless the external medium can provide a significant non-zero pressure. Subsonic solutions with $|{\dot{M}}|$ less than the critical value are usually unacceptable for similar reasons. (In fact the interstellar medium does arrest the supersonic solar wind in a termination shock well beyond Pluto’s orbit.)

8.6 Accretion

In spherical or Bondi (1952) ¹⁴ accretion we consider a gas that is uniform and at rest at infinity (with pressure $p_{0}$ and density ${\it\rho}_{0}$ ). Then $B=v_{s0}^{2}/({\it\gamma}-1)$ and $v_{ss}^{2}=2v_{s0}^{2}/(5-3{\it\gamma})$ . The appropriate solution is subsonic (uniform) at large $r$ and supersonic (freely falling) at small $r$ . If the accreting object has a dense surface (a star rather than a black hole) then the accretion flow will be arrested by a shock above the surface.

The accretion rate of the critical solution is

(8.16) $$\begin{eqnarray}{\dot{M}}=4{\rm\pi}r_{s}^{2}{\it\rho}_{s}v_{ss}=4{\rm\pi}r_{s}^{2}{\it\rho}_{0}v_{s0}\left(\frac{v_{ss}}{v_{s0}}\right)^{({\it\gamma}+1)/({\it\gamma}-1)}=f({\it\gamma}){\dot{M}}_{B},\end{eqnarray}$$

where

(8.17) $$\begin{eqnarray}{\dot{M}}_{B}=\frac{{\rm\pi}G^{2}M^{2}{\it\rho}_{0}}{v_{s0}^{3}}=4{\rm\pi}r_{a}^{2}{\it\rho}_{0}v_{s0}\end{eqnarray}$$

is the characteristic Bondi accretion rate and

(8.18) $$\begin{eqnarray}f({\it\gamma})=\left(\frac{2}{5-3{\it\gamma}}\right)^{(5-3{\it\gamma})/2({\it\gamma}-1)}\end{eqnarray}$$

is a dimensionless factor. Here

(8.19) $$\begin{eqnarray}r_{a}=\frac{GM}{2v_{s0}^{2}}\end{eqnarray}$$

is the nominal accretion radius, roughly the radius within which the mass $M$ captures the surrounding medium into a supersonic inflow.

Exercise: show that

(8.20a,b ) $$\begin{eqnarray}\lim _{{\it\gamma}\rightarrow 1}f({\it\gamma})=\text{e}^{3/2},\quad \lim _{{\it\gamma}\rightarrow 5/3}f({\it\gamma})=1.\end{eqnarray}$$

(However, although the case ${\it\gamma}=1$ admits a sonic point, the important case ${\it\gamma}=5/3$ does not.)

At different times in its life a star may gain mass from, or lose mass to, its environment. Currently the Sun is losing mass at an average rate of approximately $2\times 10^{-14}~M_{\odot }~\text{yr}^{-1}$ . If it were not doing so, it could theoretically accrete at the Bondi rate of approximately $3\times 10^{-15}~M_{\odot }~\text{yr}^{-1}$ from the interstellar medium.

Related examples: A.17–A.19.

9 Axisymmetric rotating magnetized flows: astrophysical jets

Note: in this section $(r,{\it\phi},z)$ are cylindrical polar coordinates.

9.1 Introduction

Stellar winds and jets from accretion discs are examples of outflows in which rotation and magnetic fields have important or essential roles. Using cylindrical polar coordinates $(r,{\it\phi},z)$ , we examine steady ( $\partial /\partial t=0$ ), axisymmetric ( $\partial /\partial {\it\phi}=0$ ) models based on the equations of ideal MHD.

9.2 Representation of an axisymmetric magnetic field

The solenoidal condition for an axisymmetric magnetic field is

(9.1) $$\begin{eqnarray}\frac{1}{r}\frac{\partial }{\partial r}(rB_{r})+\frac{\partial B_{z}}{\partial z}=0.\end{eqnarray}$$

We may write

(9.2a,b ) $$\begin{eqnarray}B_{r}=-\frac{1}{r}\frac{\partial {\it\psi}}{\partial z},\quad B_{z}=\frac{1}{r}\frac{\partial {\it\psi}}{\partial r},\end{eqnarray}$$

where ${\it\psi}(r,z)$ is the magnetic flux function (figure 15). This is related to the magnetic vector potential by ${\it\psi}=rA_{{\it\phi}}$ . The magnetic flux contained inside the circle ( $r=\text{const.}$ , $z=\text{const.}$ ) is

(9.3) $$\begin{eqnarray}\int _{0}^{r}B_{z}(r^{\prime },z)\,2{\rm\pi}r^{\prime }\,\text{d}r^{\prime }=2{\rm\pi}{\it\psi}(r,z),\end{eqnarray}$$

plus an arbitrary constant that can be set to zero.

Figure 15. Magnetic flux function and poloidal magnetic field.

Since $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}=0$ , ${\it\psi}$ labels the magnetic field lines or their surfaces of revolution, known as magnetic surfaces. The magnetic field may be written in the form

(9.4) $$\begin{eqnarray}\boldsymbol{B}=\boldsymbol{{\rm\nabla}}{\it\psi}\times \boldsymbol{{\rm\nabla}}{\it\phi}+B_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}}=\left[-\frac{1}{r}\boldsymbol{e}_{{\it\phi}}\times \boldsymbol{{\rm\nabla}}{\it\psi}\right]+\left[B_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}}\right].\end{eqnarray}$$

The two square brackets represent the poloidal (meridional) and toroidal (azimuthal) parts of the magnetic field:

(9.5) $$\begin{eqnarray}\boldsymbol{B}=\boldsymbol{B}_{p}+B_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}}.\end{eqnarray}$$

Note that

(9.6) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}_{p}=0.\end{eqnarray}$$

Similarly, one can write the velocity in the form

(9.7) $$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{u}_{p}+u_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}},\end{eqnarray}$$

although $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{p}\neq 0$ in general.

9.3 Mass loading and angular velocity

The steady induction equation in ideal MHD,

(9.8) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B})=\mathbf{0},\end{eqnarray}$$

implies

(9.9) $$\begin{eqnarray}\boldsymbol{u}\times \boldsymbol{B}=-\boldsymbol{E}=\boldsymbol{{\rm\nabla}}{\it\Phi}_{e},\end{eqnarray}$$

where ${\it\Phi}_{e}$ is the electrostatic potential. Now

(9.10) $$\begin{eqnarray}\displaystyle \boldsymbol{u}\times \boldsymbol{B} & = & \displaystyle (\boldsymbol{u}_{p}+u_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}})\times (\boldsymbol{B}_{p}+B_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}})\nonumber\\ \displaystyle & = & \displaystyle \left[\boldsymbol{e}_{{\it\phi}}\times (u_{{\it\phi}}\boldsymbol{B}_{p}-B_{{\it\phi}}\boldsymbol{u}_{p})\right]+\left[\boldsymbol{u}_{p}\times \boldsymbol{B}_{p}\right].\end{eqnarray}$$

For an axisymmetric solution with $\partial {\it\Phi}_{\text{e}}/\partial {\it\phi}=0$ , we have

(9.11) $$\begin{eqnarray}\boldsymbol{u}_{p}\times \boldsymbol{B}_{p}=\mathbf{0},\end{eqnarray}$$

i.e. the poloidal velocity is parallel to the poloidal magnetic field ¹⁵ . Let

(9.12) $$\begin{eqnarray}{\it\rho}\boldsymbol{u}_{p}=k\boldsymbol{B}_{p},\end{eqnarray}$$

where $k$ is the mass loading, i.e. the ratio of mass flux to magnetic flux.

The steady equation of mass conservation is

(9.13) $$\begin{eqnarray}0=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u})=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u}_{p})=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(k\boldsymbol{B}_{p})=\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}k.\end{eqnarray}$$

Therefore

(9.14) $$\begin{eqnarray}k=k({\it\psi}),\end{eqnarray}$$

i.e.  $k$ is a surface function, constant on each magnetic surface.

We now have

(9.15) $$\begin{eqnarray}\boldsymbol{u}\times \boldsymbol{B}=\boldsymbol{e}_{{\it\phi}}\times (u_{{\it\phi}}\boldsymbol{B}_{p}-B_{{\it\phi}}\boldsymbol{u}_{p})=\left(\frac{u_{{\it\phi}}}{r}-\frac{kB_{{\it\phi}}}{r{\it\rho}}\right)\boldsymbol{{\rm\nabla}}{\it\psi}.\end{eqnarray}$$

Taking the curl of this equation, we find

(9.16) $$\begin{eqnarray}\mathbf{0}=\boldsymbol{{\rm\nabla}}\left(\frac{u_{{\it\phi}}}{r}-\frac{kB_{{\it\phi}}}{r{\it\rho}}\right)\times \boldsymbol{{\rm\nabla}}{\it\psi}.\end{eqnarray}$$

Therefore

(9.17) $$\begin{eqnarray}\frac{u_{{\it\phi}}}{r}-\frac{kB_{{\it\phi}}}{r{\it\rho}}={\it\omega},\end{eqnarray}$$

where ${\it\omega}({\it\psi})$ is another surface function, known as the angular velocity of the magnetic surface.

The complete velocity field may be written in the form

(9.18) $$\begin{eqnarray}\boldsymbol{u}=\frac{k\boldsymbol{B}}{{\it\rho}}+r{\it\omega}\boldsymbol{e}_{{\it\phi}},\end{eqnarray}$$

i.e. the total velocity is parallel to the total magnetic field in a frame of reference rotating with angular velocity ${\it\omega}$ . It is useful to think of the fluid being constrained to move along the field line like a bead on a rotating wire.

9.4 Entropy

The steady thermal energy equation,

(9.19) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s=0,\end{eqnarray}$$

implies that $\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s=0$ and so

(9.20) $$\begin{eqnarray}s=s({\it\psi})\end{eqnarray}$$

is another surface function.

9.5 Angular momentum

The azimuthal component of the equation of motion is

(9.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\rho}\left(\boldsymbol{u}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}u_{{\it\phi}}+\frac{u_{r}u_{{\it\phi}}}{r}\right)=\frac{1}{{\it\mu}_{0}}\left(\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}B_{{\it\phi}}+\frac{B_{r}B_{{\it\phi}}}{r}\right)\\ \displaystyle \frac{1}{r}{\it\rho}\boldsymbol{u}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(ru_{{\it\phi}})-\frac{1}{{\it\mu}_{0}r}\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(rB_{{\it\phi}})=0\\ \displaystyle \frac{1}{r}\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\left(kru_{{\it\phi}}-\frac{rB_{{\it\phi}}}{{\it\mu}_{0}}\right)=0,\end{array}\right\}\end{eqnarray}$$

and so

(9.22) $$\begin{eqnarray}ru_{{\it\phi}}=\frac{rB_{{\it\phi}}}{{\it\mu}_{0}k}+\ell ,\end{eqnarray}$$

where

(9.23) $$\begin{eqnarray}\ell =\ell ({\it\psi})\end{eqnarray}$$

is another surface function, the angular momentum invariant. This is the angular momentum removed in the outflow per unit mass, although part of the torque is carried by the magnetic field.

9.6 The Alfvén surface

Define the poloidal Alfvén number (cf. the Mach number)

(9.24) $$\begin{eqnarray}A=\frac{u_{p}}{v_{ap}}.\end{eqnarray}$$

Then

(9.25) $$\begin{eqnarray}A^{2}=\frac{{\it\mu}_{0}{\it\rho}u_{p}^{2}}{B_{p}^{2}}=\frac{{\it\mu}_{0}k^{2}}{{\it\rho}},\end{eqnarray}$$

and so $A\propto {\it\rho}^{-1/2}$ on each magnetic surface.

Consider the two equations

(9.26a,b ) $$\begin{eqnarray}\frac{u_{{\it\phi}}}{r}=\frac{kB_{{\it\phi}}}{r{\it\rho}}+{\it\omega},\quad ru_{{\it\phi}}=\frac{rB_{{\it\phi}}}{{\it\mu}_{0}k}+\ell .\end{eqnarray}$$

Eliminate $B_{{\it\phi}}$ to obtain

(9.27) $$\begin{eqnarray}u_{{\it\phi}}=\frac{r^{2}{\it\omega}-A^{2}\ell }{r(1-A^{2})}=\left(\frac{1}{1-A^{2}}\right)r{\it\omega}+\left(\frac{A^{2}}{A^{2}-1}\right)\frac{\ell }{r}.\end{eqnarray}$$

For $A\ll 1$ we have

(9.28) $$\begin{eqnarray}u_{{\it\phi}}\approx r{\it\omega},\end{eqnarray}$$

i.e. the fluid is in uniform rotation, corotating with the magnetic surface. For $A\gg 1$ we have

(9.29) $$\begin{eqnarray}u_{{\it\phi}}\approx \frac{\ell }{r},\end{eqnarray}$$

i.e. the fluid conserves its specific angular momentum. The point $r=r_{a}({\it\psi})$ where $A=1$ is the Alfvén point. The locus of Alfvén points for different magnetic surfaces forms the Alfvén surface. To avoid a singularity there we require

(9.30) $$\begin{eqnarray}\ell =r_{a}^{2}{\it\omega}.\end{eqnarray}$$

Typically the outflow will start at low velocity in high-density material, where $A\ll 1$ . We can therefore identify ${\it\omega}$ as the angular velocity $u_{{\it\phi}}/r={\it\Omega}_{0}$ of the footpoint $r=r_{0}$ of the magnetic field line at the source of the outflow. It will then accelerate smoothly through an Alfvén surface and become super-Alfvénic ( $A>1$ ). If mass is lost at a rate ${\dot{M}}$ in the outflow, angular momentum is lost at a rate ${\dot{M}}\ell ={\dot{M}}r_{a}^{2}{\it\Omega}_{0}$ . In contrast, in a hydrodynamic outflow, angular momentum is conserved by fluid elements and is therefore lost at a rate ${\dot{M}}r_{0}^{2}{\it\Omega}_{0}$ . A highly efficient removal of angular momentum occurs if the Alfvén radius is large compared to the footpoint radius. This effect is the magnetic lever arm. The loss of angular momentum through a stellar wind is called magnetic braking (figure 16). In the case of the Sun, the Alfvén radius is approximately between $20$ and $30\,R_{\odot }$ .

Figure 16. Acceleration through an Alfvén point along a poloidal magnetic field line, leading to angular momentum loss and magnetic braking.

9.7 The Bernoulli function

The total energy equation for a steady flow is

(9.31) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\rho}\boldsymbol{u}\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)+\frac{\boldsymbol{E}\times \boldsymbol{B}}{{\it\mu}_{0}}\right]=0.\end{eqnarray}$$

Now since

(9.32) $$\begin{eqnarray}\boldsymbol{u}=\frac{k\boldsymbol{B}}{{\it\rho}}+r{\it\omega}\boldsymbol{e}_{{\it\phi}},\end{eqnarray}$$

we have

(9.33) $$\begin{eqnarray}\boldsymbol{E}=-\boldsymbol{u}\times \boldsymbol{B}=-r{\it\omega}\boldsymbol{e}_{{\it\phi}}\times \boldsymbol{B}=-r{\it\omega}\boldsymbol{e}_{{\it\phi}}\times \boldsymbol{B}_{p},\end{eqnarray}$$

which is purely poloidal. Thus

(9.34) $$\begin{eqnarray}(\boldsymbol{E}\times \boldsymbol{B})_{p}=\boldsymbol{E}\times (B_{{\it\phi}}\,\boldsymbol{e}_{{\it\phi}})=-r{\it\omega}B_{{\it\phi}}\boldsymbol{B}_{p}.\end{eqnarray}$$

The total energy equation is therefore

(9.35) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[k\boldsymbol{B}_{p}\left(\frac{1}{2}u^{2}+{\it\Phi}+h\right)-\frac{r{\it\omega}B_{{\it\phi}}}{{\it\mu}_{0}}\boldsymbol{B}_{p}\right]=0\\ \displaystyle \boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\left[k\left(\frac{1}{2}u^{2}+{\it\Phi}+h-\frac{r{\it\omega}B_{{\it\phi}}}{{\it\mu}_{0}k}\right)\right]=0\\ \displaystyle \frac{1}{2}u^{2}+{\it\Phi}+h-\frac{r{\it\omega}B_{{\it\phi}}}{{\it\mu}_{0}k}={\it\varepsilon},\end{array}\right\}\end{eqnarray}$$

where

(9.36) $$\begin{eqnarray}{\it\varepsilon}={\it\varepsilon}({\it\psi})\end{eqnarray}$$

is another surface function, the energy invariant.

An alternative invariant is

(9.37) $$\begin{eqnarray}\displaystyle \tilde{{\it\varepsilon}} & = & \displaystyle {\it\varepsilon}-\ell {\it\omega}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}u^{2}+{\it\Phi}+h-\frac{r{\it\omega}B_{{\it\phi}}}{{\it\mu}_{0}k}-\left(ru_{{\it\phi}}-\frac{rB_{{\it\phi}}}{{\it\mu}_{0}k}\right){\it\omega}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}u^{2}+{\it\Phi}+h-ru_{{\it\phi}}{\it\omega}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}u_{p}^{2}+\frac{1}{2}(u_{{\it\phi}}-r{\it\omega})^{2}+{\it\Phi}_{cg}+h,\end{eqnarray}$$

where

(9.38) $$\begin{eqnarray}{\it\Phi}_{cg}={\it\Phi}-{\textstyle \frac{1}{2}}{\it\omega}^{2}r^{2}\end{eqnarray}$$

is the centrifugal–gravitational potential associated with the magnetic surface. One can then see that $\tilde{{\it\varepsilon}}$ is the Bernoulli function of the flow in the frame rotating with angular velocity ${\it\omega}$ . In this frame the flow is strictly parallel to the magnetic field and the field therefore does no work because $\boldsymbol{J}\times \boldsymbol{B}\bot \boldsymbol{B}$ and so $\boldsymbol{J}\times \boldsymbol{B}\bot (\boldsymbol{u}-r{\it\omega}\boldsymbol{e}_{{\it\phi}})$ .

9.8 Summary

We have been able to integrate almost all of the MHD equations, reducing them to a set of algebraic relations on each magnetic surface. If the poloidal magnetic field $\boldsymbol{B}_{p}$ (or, equivalently, the flux function ${\it\psi}$ ) is specified in advance, these algebraic equations are sufficient to determine the complete solution on each magnetic surface separately, although we must also (i) specify the initial conditions at the source of the outflow and (ii) ensure that the solution passes smoothly through critical points where the flow speed matches the speeds of slow and fast magnetoacoustic waves (see Example A.21).

The component of the equation of motion perpendicular to the magnetic surfaces is the only piece of information not yet used. This ‘transfield’ or ‘Grad–Shafranov’ equation ultimately determines the equilibrium shape of the magnetic surfaces. It is a very complicated nonlinear partial differential equation for ${\it\psi}(r,z)$ and cannot be reduced to simple terms. We do not consider it here.

9.9 Acceleration from the surface of an accretion disc

We now consider the launching of an outflow from a thin accretion disc. The angular velocity ${\it\Omega}(r)$ of the disc corresponds approximately to circular Keplerian orbital motion around a central mass  $M$ :

(9.39) $$\begin{eqnarray}{\it\Omega}\approx \left(\frac{GM}{r^{3}}\right)^{1/2}.\end{eqnarray}$$

If the flow starts essentially from rest in high-density material ( $A\ll 1$ ), we have

(9.40) $$\begin{eqnarray}{\it\omega}\approx {\it\Omega},\end{eqnarray}$$

i.e. the angular velocity of the magnetic surface is the angular velocity of the disc at the footpoint of the field line. In the sub-Alfvénic region we have

(9.41) $$\begin{eqnarray}\tilde{{\it\varepsilon}}\approx {\textstyle \frac{1}{2}}u_{p}^{2}+{\it\Phi}_{cg}+h.\end{eqnarray}$$

As in the case of stellar winds, if the gas is hot (comparable to the escape temperature) an outflow can be driven by thermal pressure. Of more interest here is the possibility of a dynamically driven outflow. For a ‘cold’ wind the enthalpy makes a negligible contribution in this equation. Whether the flow accelerates or not above the disc then depends on the variation of the centrifugal–gravitational potential along the field line.

Consider a Keplerian disc in a point-mass potential. Let the footpoint of the field line be at $r=r_{0}$ , and let the angular velocity of the field line be

(9.42) $$\begin{eqnarray}{\it\omega}={\it\Omega}_{0}=\left(\frac{GM}{r_{0}^{3}}\right)^{1/2},\end{eqnarray}$$

as argued above. Then

(9.43) $$\begin{eqnarray}{\it\Phi}_{cg}=-GM(r^{2}+z^{2})^{-1/2}-\frac{1}{2}\frac{GM}{r_{0}^{3}}r^{2}.\end{eqnarray}$$

Figure 17. Contours of ${\it\Phi}_{cg}$ , in units such that $r_{0}=1$ . The downhill directions are indicated by dotted contours. If the inclination of the poloidal magnetic field to the vertical direction at the surface of the disc exceeds $30^{\circ }$ , gas is accelerated along the field lines away from the disc.

In units such that $r_{0}=1$ , the equation of the equipotential passing through the footpoint $(r_{0},z)$ is

(9.44) $$\begin{eqnarray}(r^{2}+z^{2})^{-1/2}+\frac{r^{2}}{2}=\frac{3}{2}.\end{eqnarray}$$

This can be rearranged into the form

(9.45) $$\begin{eqnarray}z^{2}=\frac{(2-r)(r-1)^{2}(r+1)^{2}(r+2)}{(3-r^{2})^{2}}.\end{eqnarray}$$

Close to the footpoint $(1,0)$ we have

(9.46) $$\begin{eqnarray}z^{2}\approx 3(r-1)^{2}\quad \Rightarrow \quad z\approx \pm \sqrt{3}(r-1).\end{eqnarray}$$

The footpoint lies at a saddle point of ${\it\Phi}_{cg}$ (figure 17). If the inclination of the field line to the vertical, $i$ , at the surface of the disc exceeds $30^{\circ }$ , the flow is accelerated without thermal assistance ¹⁶ . This is magnetocentrifugal acceleration.

The critical equipotential has an asymptote at $r=r_{0}\sqrt{3}$ . The field line must continue to expand sufficiently in the radial direction in order to sustain the magnetocentrifugal acceleration.

9.10 Magnetically driven accretion

To allow a quantity of mass ${\rm\Delta}M_{acc}$ to be accreted from radius $r_{0}$ , its orbital angular momentum $r_{0}^{2}{\it\Omega}_{0}\,{\rm\Delta}M_{acc}$ must be removed. The angular momentum removed by a quantity of mass ${\rm\Delta}M_{jet}$ flowing out in a magnetized jet from radius $r_{0}$ is $\ell \,{\rm\Delta}M_{jet}=r_{a}^{2}{\it\Omega}_{0}\,{\rm\Delta}M_{jet}$ . Therefore accretion can in principle be driven by an outflow, with

(9.47) $$\begin{eqnarray}\frac{{\dot{M}}_{acc}}{{\dot{M}}_{jet}}\approx \frac{r_{a}^{2}}{r_{0}^{2}}.\end{eqnarray}$$

The magnetic lever arm allows an efficient removal of angular momentum if the Alfvén radius is large compared to the footpoint radius.

Related examples: A.20, A.21.

10 Lagrangian formulation of ideal MHD

10.1 The Lagrangian viewpoint

In § 11 we will discuss waves and instabilities in differentially rotating astrophysical bodies. Here we develop a general theory of disturbances to fluid flows that makes use of the conserved quantities in ideal fluids and takes a Lagrangian approach.

The flow of a fluid can be considered as a time-dependent map,

(10.1) $$\begin{eqnarray}\boldsymbol{a}\mapsto \boldsymbol{x}(\boldsymbol{a},t),\end{eqnarray}$$

where $\boldsymbol{a}$ is the position vector of a fluid element at some initial time $t_{0}$ , and $\boldsymbol{x}$ is its position vector at time $t$ . The Cartesian components of $\boldsymbol{a}$ are examples of Lagrangian variables, labelling the fluid element. The components of $\boldsymbol{x}$ are Eulerian variables, labelling a fixed point in space. Any fluid property $X$ (scalar, vector or tensor) can be regarded as a function of either Lagrangian or Eulerian variables:

(10.2) $$\begin{eqnarray}X=X^{L}(\boldsymbol{a},t)=X^{E}(\boldsymbol{x},t).\end{eqnarray}$$

The Lagrangian time-derivative is simply

(10.3) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}=\left(\frac{\partial }{\partial t}\right)_{\boldsymbol{a}},\end{eqnarray}$$

and the velocity of the fluid is

(10.4) $$\begin{eqnarray}\boldsymbol{u}=\frac{\text{D}\boldsymbol{x}}{\text{D}t}=\left(\frac{\partial \boldsymbol{x}}{\partial t}\right)_{\boldsymbol{a}}.\end{eqnarray}$$

The aim of a Lagrangian formulation of ideal MHD is to derive a nonlinear evolutionary equation for the function $\boldsymbol{x}(\boldsymbol{a},t)$ . The dynamics is Hamiltonian in character and can be derived from a Lagrangian function or action principle. There are many similarities with classical field theories.

10.2 The deformation tensor

Introduce the deformation tensor of the flow,

(10.5) $$\begin{eqnarray}\unicode[STIX]{x1D60D}_{ij}=\frac{\partial x_{i}}{\partial a_{j}},\end{eqnarray}$$

and its determinant

(10.6) $$\begin{eqnarray}F=\det (\unicode[STIX]{x1D60D}_{ij})\end{eqnarray}$$

and inverse

(10.7) $$\begin{eqnarray}\unicode[STIX]{x1D60E}_{ij}=\frac{\partial a_{i}}{\partial x_{j}}.\end{eqnarray}$$

We note the following properties. First, the derivative

(10.8) $$\begin{eqnarray}\frac{\partial F}{\partial \unicode[STIX]{x1D60D}_{ij}}=\unicode[STIX]{x1D60A}_{ij}=F\unicode[STIX]{x1D60E}_{ji}=\frac{1}{2}{\it\epsilon}_{ik\ell }{\it\epsilon}_{jmn}\unicode[STIX]{x1D60D}_{km}F_{\ell n}\end{eqnarray}$$

is equal to the cofactor $\unicode[STIX]{x1D60A}_{ij}$ of the matrix element $\unicode[STIX]{x1D60D}_{ij}$ . (This follows from the fact that the determinant can be expanded as the sum of the products of any row’s elements with their cofactors, which do not depend on that row’s elements.) Second, the matrix of cofactors has zero divergence on its second index:

(10.9) $$\begin{eqnarray}\frac{\partial \unicode[STIX]{x1D60A}_{ij}}{\partial a_{j}}=\frac{\partial }{\partial a_{j}}\left(\frac{1}{2}{\it\epsilon}_{ik\ell }{\it\epsilon}_{jmn}\frac{\partial x_{k}}{\partial a_{m}}\frac{\partial x_{\ell }}{\partial a_{n}}\right)=0.\end{eqnarray}$$

(This follows because the resulting derivative involves the contraction of the antisymmetric tensor ${\it\epsilon}_{jmn}$ with expressions that are symmetric in either $jm$ or $jn$ .)

Now

(10.10) $$\begin{eqnarray}\frac{\text{D}\unicode[STIX]{x1D60D}_{ij}}{\text{D}t}=\frac{\partial u_{i}}{\partial a_{j}}\end{eqnarray}$$

and, according to (10.8),

(10.11) $$\begin{eqnarray}\frac{\text{D}\ln F}{\text{D}t}=\unicode[STIX]{x1D60E}_{ji}\frac{\text{D}\unicode[STIX]{x1D60D}_{ij}}{\text{D}t}=\frac{\partial a_{j}}{\partial x_{i}}\frac{\partial u_{i}}{\partial a_{j}}=\frac{\partial u_{i}}{\partial x_{i}}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}.\end{eqnarray}$$

10.3 Geometrical conservation laws

The equations of ideal MHD comprise the equation of motion and three ‘geometrical’ conservation laws. These are the conservation of specific entropy (thermal energy equation),

(10.12) $$\begin{eqnarray}\frac{\text{D}s}{\text{D}t}=0,\end{eqnarray}$$

the conservation of mass,

(10.13) $$\begin{eqnarray}\frac{\text{D}{\it\rho}}{\text{D}t}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u},\end{eqnarray}$$

and the conservation of magnetic flux (induction equation),

(10.14) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}\left(\frac{\boldsymbol{B}}{{\it\rho}}\right)=\left(\frac{\boldsymbol{B}}{{\it\rho}}\right)\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}.\end{eqnarray}$$

These equations describe the pure advection of fluid properties in a manner equivalent to the advection of various geometrical objects. The specific entropy is advected as a simple scalar, so that its numerical value is conserved by material points. The specific volume $v=1/{\it\rho}$ is advected in the same way as a material volume element $\text{d}V$ . The quantity $\boldsymbol{B}/{\it\rho}$ is advected in the same way as a material line element ${\it\delta}\boldsymbol{x}$ . Equivalently, the mass element ${\it\delta}m={\it\rho}{\it\delta}V$ and the magnetic flux element ${\it\delta}{\it\Phi}=\boldsymbol{B}\boldsymbol{\cdot }{\it\delta}\boldsymbol{S}$ satisfy $\text{D}{\it\delta}m/\text{D}t=0$ and $\text{D}{\it\delta}{\it\Phi}/\text{D}t=0$ , where ${\it\delta}\boldsymbol{S}$ is a material surface element. All three conservation laws can be integrated exactly in Lagrangian variables.

The exact solutions of (10.12)–(10.14) are then

(10.15a-c ) $$\begin{eqnarray}s^{L}(\boldsymbol{a},t)=s_{0}(\boldsymbol{a}),\quad {\it\rho}^{L}(\boldsymbol{a},t)=F^{-1}{\it\rho}_{0}(\boldsymbol{a}),\quad B_{i}^{L}(\boldsymbol{a},t)=F^{-1}\unicode[STIX]{x1D60D}_{ij}B_{j0}(\boldsymbol{a}),\end{eqnarray}$$

where $s_{0}$ , ${\it\rho}_{0}$ , and $\boldsymbol{B}_{0}$ are the initial values at time $t_{0}$ . The verification of (10.14) is

(10.16) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}\left(\frac{B_{i}}{{\it\rho}}\right)=\frac{\text{D}}{\text{D}t}\left(\frac{\unicode[STIX]{x1D60D}_{ij}B_{j0}}{{\it\rho}_{0}}\right)=\frac{\partial u_{i}}{\partial a_{j}}\frac{B_{j0}}{{\it\rho}_{0}}=\frac{\partial u_{i}}{\partial x_{k}}\unicode[STIX]{x1D60D}_{kj}\frac{B_{j0}}{{\it\rho}_{0}}=\left(\frac{B_{k}}{{\it\rho}}\right)\frac{\partial u_{i}}{\partial x_{k}}.\end{eqnarray}$$

Note that the advected quantities at time $t$ depend only on the initial values and on the instantaneous mapping $\boldsymbol{a}\mapsto \boldsymbol{x}$ , not on the intermediate history of the flow. The ‘memory’ of an ideal fluid is perfect.

10.4 The Lagrangian of ideal MHD

Newtonian dynamics can be formulated using Hamilton’s principle of stationary action,

(10.17) $$\begin{eqnarray}{\it\delta}\int L\,\text{d}t=0,\end{eqnarray}$$

where the Lagrangian $L$ is the difference between the kinetic energy and the potential energy of the system. By analogy, we may expect the Lagrangian of ideal MHD to take the form

(10.18) $$\begin{eqnarray}L=\int {\mathcal{L}}\,\text{d}V,\end{eqnarray}$$

where (for a non-self-gravitating fluid)

(10.19) $$\begin{eqnarray}{\mathcal{L}}={\it\rho}\left(\frac{1}{2}u^{2}-{\it\Phi}-e-\frac{B^{2}}{2{\it\mu}_{0}{\it\rho}}\right)\end{eqnarray}$$

is the Lagrangian density.

To verify this, we assume that the equation of state can be written in the form $e=e(v,s)$ , where $v={\it\rho}^{-1}$ is the specific volume. Since $\text{d}e=T\,\text{d}s-p\,\text{d}v$ , we have

(10.20a,b ) $$\begin{eqnarray}\left(\frac{\partial e}{\partial v}\right)_{s}=-p,\quad \left(\frac{\partial ^{2}e}{\partial v^{2}}\right)_{s}=\frac{{\it\gamma}p}{v}\end{eqnarray}$$

(strictly, ${\it\gamma}$ should be ${\it\Gamma}_{1}$ here).

We then write the action using Lagrangian variables,

(10.21) $$\begin{eqnarray}S[\boldsymbol{x}]=\iint \tilde{{\mathcal{L}}}(\boldsymbol{x},\boldsymbol{u},\unicode[STIX]{x1D641})\,\text{d}^{3}\boldsymbol{a}\,\text{d}t,\end{eqnarray}$$

with

(10.22) $$\begin{eqnarray}\tilde{{\mathcal{L}}}={\it\rho}_{0}\left[\frac{1}{2}u^{2}-{\it\Phi}(\boldsymbol{x})-e(F{\it\rho}_{0}^{-1},s_{0})-\frac{F^{-1}\unicode[STIX]{x1D60D}_{ij}B_{j0}\unicode[STIX]{x1D60D}_{ik}B_{k0}}{2{\it\mu}_{0}{\it\rho}_{0}}\right].\end{eqnarray}$$

This uses the fact that $F$ is the Jacobian determinant of the transformation $\boldsymbol{a}\mapsto \boldsymbol{x}$ , or, equivalently, that ${\it\rho}\,\text{d}^{3}\boldsymbol{x}={\it\rho}_{0}\,\text{d}^{3}\boldsymbol{a}=\text{d}m$ is an invariant mass measure. $\tilde{{\mathcal{L}}}$ is now expressed in terms of the function $\boldsymbol{x}(\boldsymbol{a},t)$ and its derivatives with respect to time ( $\boldsymbol{u}$ ) and space ( $\unicode[STIX]{x1D641}$ ). The Euler–Lagrange equation for the variational principle ${\it\delta}S=0$ is

(10.23) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}\frac{\partial \tilde{{\mathcal{L}}}}{\partial u_{i}}+\frac{\partial }{\partial a_{j}}\frac{\partial \tilde{{\mathcal{L}}}}{\partial \unicode[STIX]{x1D60D}_{ij}}-\frac{\partial \tilde{{\mathcal{L}}}}{\partial x_{i}}=0.\end{eqnarray}$$

The straightforward terms are

(10.24a,b ) $$\begin{eqnarray}\frac{\partial \tilde{{\mathcal{L}}}}{\partial u_{i}}={\it\rho}_{0}u_{i},\quad \frac{\partial \tilde{{\mathcal{L}}}}{\partial x_{i}}=-{\it\rho}_{0}\frac{\partial {\it\Phi}}{\partial x_{i}}.\end{eqnarray}$$

Now

(10.25) $$\begin{eqnarray}\displaystyle \frac{\partial \tilde{{\mathcal{L}}}}{\partial \unicode[STIX]{x1D60D}_{ij}} & = & \displaystyle \left(p+\frac{B^{2}}{2{\it\mu}_{0}}\right)\frac{\partial F}{\partial \unicode[STIX]{x1D60D}_{ij}}-\frac{F^{-1}B_{j0}\unicode[STIX]{x1D60D}_{ik}B_{k0}}{{\it\mu}_{0}}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D60A}_{ij}\left(p+\frac{B^{2}}{2{\it\mu}_{0}}\right)-\frac{1}{{\it\mu}_{0}}\unicode[STIX]{x1D60A}_{kj}B_{i}B_{k}\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D60A}_{kj}\unicode[STIX]{x1D61D}_{ik},\end{eqnarray}$$

where

(10.26) $$\begin{eqnarray}\unicode[STIX]{x1D61D}_{ik}=-\left(p+\frac{B^{2}}{2{\it\mu}_{0}}\right){\it\delta}_{ik}+\frac{B_{i}B_{k}}{{\it\mu}_{0}}\end{eqnarray}$$

is the stress tensor due to pressure and the magnetic field.

The Euler–Lagrange equation is therefore

(10.27) $$\begin{eqnarray}{\it\rho}_{0}\frac{\text{D}u_{i}}{\text{D}t}=-{\it\rho}_{0}\frac{\partial {\it\Phi}}{\partial x_{i}}+\frac{\partial }{\partial a_{j}}(\unicode[STIX]{x1D60A}_{kj}\unicode[STIX]{x1D61D}_{ik}).\end{eqnarray}$$

Using (10.9) we note that

(10.28) $$\begin{eqnarray}\frac{\partial }{\partial a_{j}}(\unicode[STIX]{x1D60A}_{kj}\unicode[STIX]{x1D61D}_{ik})=\unicode[STIX]{x1D60A}_{kj}\frac{\partial \unicode[STIX]{x1D61D}_{ik}}{\partial a_{j}}=F\unicode[STIX]{x1D60E}_{jk}\frac{\partial \unicode[STIX]{x1D61D}_{ik}}{\partial a_{j}}=F\frac{\partial \unicode[STIX]{x1D61D}_{ik}}{\partial x_{k}}.\end{eqnarray}$$

On dividing through by $F$ , the Euler–Lagrange equation becomes the desired equation of motion,

(10.29) $$\begin{eqnarray}{\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D651}.\end{eqnarray}$$

In this construction, the fluid flow is viewed as a field $\boldsymbol{x}(\boldsymbol{a},t)$ on the initial state space. Ideal MHD is seen as a nonlinear field theory derived from an action principle. When considering stability problems, it is useful to generalize this concept and to view a perturbed flow as a field on an unperturbed flow.

10.5 The Lagrangian displacement

Now consider two different flows, $\boldsymbol{x}(\boldsymbol{a},t)$ and $\hat{\boldsymbol{x}}(\boldsymbol{a},t)$ , for which the initial values of the advected quantities, $s_{0}$ , ${\it\rho}_{0}$ and $\boldsymbol{B}_{0}$ , are the same. The two deformation tensors are related by the chain rule,

(10.30) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D60D}}_{ij}=\unicode[STIX]{x1D611}_{ik}\unicode[STIX]{x1D60D}_{kj},\end{eqnarray}$$

where

(10.31) $$\begin{eqnarray}\unicode[STIX]{x1D611}_{ik}=\frac{\partial \hat{x}_{i}}{\partial x_{k}}\end{eqnarray}$$

is the Jacobian matrix of the map $\boldsymbol{x}\mapsto \hat{\boldsymbol{x}}$ . Similarly,

(10.32) $$\begin{eqnarray}\hat{F}=JF,\end{eqnarray}$$

where

(10.33) $$\begin{eqnarray}J=\det (\unicode[STIX]{x1D611}_{ij})\end{eqnarray}$$

is the Jacobian determinant. The advected quantities in the two flows are therefore related by the composition of maps,

(10.34) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\hat{s}}^{L}(\boldsymbol{a},t)=s^{L}(\boldsymbol{a},t),\\ \displaystyle \hat{{\it\rho}}^{L}(\boldsymbol{a},t)=J^{-1}{\it\rho}^{L}(\boldsymbol{a},t),\\ \displaystyle \hat{B}_{i}^{L}(\boldsymbol{a},t)=J^{-1}\unicode[STIX]{x1D611}_{ij}B_{j}^{L}(\boldsymbol{a},t).\end{array}\right\}\end{eqnarray}$$

The Lagrangian displacement (figure 18) is the relative displacement of the fluid element in the two flows,

(10.35) $$\begin{eqnarray}{\bf\xi}=\hat{\boldsymbol{x}}-\boldsymbol{x}.\end{eqnarray}$$

Figure 18. The Lagrangian displacement of a fluid element.

Thus (with ${\it\xi}_{i,j}=\partial {\it\xi}_{i}/\partial x_{j}$ )

(10.36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D611}_{ij} & = & \displaystyle {\it\delta}_{ij}+{\it\xi}_{i,j},\nonumber\\ \displaystyle J & = & \displaystyle {\textstyle \frac{1}{6}}{\it\epsilon}_{ijk}{\it\epsilon}_{lmn}\unicode[STIX]{x1D611}_{il}\unicode[STIX]{x1D611}_{jm}\unicode[STIX]{x1D611}_{kn}\nonumber\\ \displaystyle & = & \displaystyle {\textstyle \frac{1}{6}}{\it\epsilon}_{ijk}{\it\epsilon}_{lmn}({\it\delta}_{il}+{\it\xi}_{i,l})({\it\delta}_{jm}+{\it\xi}_{j,m})({\it\delta}_{kn}+{\it\xi}_{k,n})\nonumber\\ \displaystyle & = & \displaystyle {\textstyle \frac{1}{6}}{\it\epsilon}_{ijk}{\it\epsilon}_{ijk}+{\textstyle \frac{1}{2}}{\it\epsilon}_{ijk}{\it\epsilon}_{ijn}{\it\xi}_{k,n}+{\textstyle \frac{1}{2}}{\it\epsilon}_{ijk}{\it\epsilon}_{imn}{\it\xi}_{j,m}{\it\xi}_{k,n}+O({\it\xi}^{3})\nonumber\\ \displaystyle & = & \displaystyle 1+{\it\xi}_{k,k}+{\textstyle \frac{1}{2}}({\it\xi}_{j,j}{\it\xi}_{k,k}-{\it\xi}_{j,k}{\it\xi}_{k,j})+O({\it\xi}^{3}).\end{eqnarray}$$

From the binomial theorem,

(10.37) $$\begin{eqnarray}J^{-1}=1-{\it\xi}_{k,k}+{\textstyle \frac{1}{2}}({\it\xi}_{j,j}{\it\xi}_{k,k}+{\it\xi}_{j,k}{\it\xi}_{k,j})+O({\it\xi}^{3}).\end{eqnarray}$$

10.6 The Lagrangian for a perturbed flow

We now use the action principle to construct a theory for the displacement as a field on the unperturbed flow: ${\bf\xi}={\bf\xi}(\boldsymbol{x},t)$ . The action for the perturbed flow is

(10.38) $$\begin{eqnarray}{\hat{S}}[{\bf\xi}]=\int \!\!\int \hat{{\mathcal{L}}}\left({\bf\xi},\frac{\partial {\bf\xi}}{\partial t},\boldsymbol{{\rm\nabla}}{\bf\xi}\right)\,\text{d}^{3}\boldsymbol{x}\,\text{d}t,\end{eqnarray}$$

where

(10.39) $$\begin{eqnarray}\displaystyle \displaystyle \hat{{\mathcal{L}}} & = & \displaystyle {\it\rho}\left(\frac{1}{2}\hat{u} ^{2}-\hat{{\it\Phi}}-\hat{e}-\frac{\hat{B}^{2}}{2{\it\mu}_{0}\hat{{\it\rho}}}\right)\nonumber\\ \displaystyle & = & \displaystyle {\it\rho}\left[\frac{1}{2}\hat{u} ^{2}-{\it\Phi}(\boldsymbol{x}+{\bf\xi})-e(J{\it\rho}^{-1},s)-\frac{J^{-1}\unicode[STIX]{x1D611}_{ij}B_{j}\unicode[STIX]{x1D611}_{ik}B_{k}}{2{\it\mu}_{0}{\it\rho}}\right],\end{eqnarray}$$

with

(10.40) $$\begin{eqnarray}\hat{\boldsymbol{u}}=\frac{\text{D}\hat{\boldsymbol{x}}}{\text{D}t}=\boldsymbol{u}+\frac{\partial {\bf\xi}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}.\end{eqnarray}$$

The Euler–Lagrange equation for the variational principle ${\it\delta}{\hat{S}}=0$ is

(10.41) $$\begin{eqnarray}\frac{\partial }{\partial t}\frac{\partial \hat{{\mathcal{L}}}}{\partial \dot{{\it\xi}}_{i}}+\frac{\partial }{\partial x_{j}}\frac{\partial \hat{{\mathcal{L}}}}{\partial {\it\xi}_{i,j}}-\frac{\partial \hat{{\mathcal{L}}}}{\partial {\it\xi}_{i}}=0,\end{eqnarray}$$

where $\dot{{\it\xi}}_{i}=\partial {\it\xi}_{i}/\partial t$ . We expand the various terms of $\hat{{\mathcal{L}}}$ in powers of ${\bf\xi}$ :

(10.42) $$\begin{eqnarray}\displaystyle \frac{1}{2}{\it\rho}\hat{u} ^{2}\! & = & \displaystyle \!\frac{1}{2}{\it\rho}u^{2}+{\it\rho}u_{i}\!\left(\frac{\partial {\it\xi}_{i}}{\partial t}+u_{j}\frac{\partial {\it\xi}_{i}}{\partial x_{j}}\right)+\frac{1}{2}{\it\rho}\!\left(\frac{\partial {\it\xi}_{i}}{\partial t}+u_{j}\frac{\partial {\it\xi}_{i}}{\partial x_{j}}\right)\!\!\left(\frac{\partial {\it\xi}_{i}}{\partial t}+u_{k}\frac{\partial {\it\xi}_{i}}{\partial x_{k}}\right)\!,\nonumber\\ \displaystyle -\,{\it\rho}{\it\Phi}(\boldsymbol{x}+{\bf\xi})\! & = & \displaystyle \!-{\it\rho}\left[{\it\Phi}(\boldsymbol{x})+{\it\xi}_{i}\frac{\partial {\it\Phi}}{\partial x_{i}}+\frac{1}{2}{\it\xi}_{i}{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+O({\it\xi}^{3})\right]\!,\nonumber\\ \displaystyle -\,{\it\rho}\left(\hat{e}+\frac{\hat{B}^{2}}{2{\it\mu}_{0}\hat{{\it\rho}}}\right)\! & = & \displaystyle \!-\left({\it\rho}e+\frac{B^{2}}{2{\it\mu}_{0}}\right)-\unicode[STIX]{x1D61D}_{ij}\frac{\partial {\it\xi}_{i}}{\partial x_{j}}-\frac{1}{2}\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{i}}{\partial x_{j}}\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}+O({\it\xi}^{3}).\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

This last expression uses the fact that $\hat{e}$ , $\hat{\boldsymbol{B}}$ and $\hat{{\it\rho}}$ depend only on $\boldsymbol{{\rm\nabla}}{\bf\xi}$ (through $J$ and $\unicode[STIX]{x1D611}_{ij}$ ) and can therefore be expanded in a Taylor series in this quantity. A short calculation of this expansion gives

(10.43) $$\begin{eqnarray}\unicode[STIX]{x1D61D}_{ij}=-\left(p+\frac{B^{2}}{2{\it\mu}_{0}}\right){\it\delta}_{ij}+\frac{B_{i}B_{j}}{{\it\mu}_{0}},\end{eqnarray}$$

which is the stress tensor used above, and

(10.44) $$\begin{eqnarray}\unicode[STIX]{x1D61D}_{ijk\ell }=\left[({\it\gamma}-1)p+\frac{B^{2}}{2{\it\mu}_{0}}\right]{\it\delta}_{ij}{\it\delta}_{k\ell }+\left(p+\frac{B^{2}}{2{\it\mu}_{0}}\right){\it\delta}_{i\ell }{\it\delta}_{jk}+\frac{1}{{\it\mu}_{0}}B_{j}B_{\ell }{\it\delta}_{ik}-\frac{1}{{\it\mu}_{0}}(B_{i}B_{j}{\it\delta}_{k\ell }+B_{k}B_{\ell }{\it\delta}_{ij}),\end{eqnarray}$$

which has the symmetry

(10.45) $$\begin{eqnarray}\unicode[STIX]{x1D61D}_{ijk\ell }=V_{k\ell ij}\end{eqnarray}$$

necessitated by its function in the Taylor series. We now have

(10.46) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial \hat{{\mathcal{L}}}}{\partial \dot{{\it\xi}}_{i}}={\it\rho}u_{i}+{\it\rho}\frac{\text{D}{\it\xi}_{i}}{\text{D}t},\\ \displaystyle \frac{\partial \hat{{\mathcal{L}}}}{\partial {\it\xi}_{i,j}}={\it\rho}u_{i}u_{j}+{\it\rho}u_{j}\frac{\text{D}{\it\xi}_{i}}{\text{D}t}-\unicode[STIX]{x1D61D}_{ij}-\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}+O({\it\xi}^{2}),\\ \displaystyle \frac{\partial \hat{{\mathcal{L}}}}{\partial {\it\xi}_{i}}=-{\it\rho}\frac{\partial {\it\Phi}}{\partial x_{i}}-{\it\rho}{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+O({\it\xi}^{2}).\end{array}\right\}\end{eqnarray}$$

Since ${\bf\xi}=\mathbf{0}$ must be a solution of the Euler–Lagrange equation, it is no surprise that the terms independent of ${\bf\xi}$ cancel by virtue of the equation of motion of the unperturbed flow,

(10.47) $$\begin{eqnarray}\frac{\partial }{\partial t}({\it\rho}u_{i})+\frac{\partial }{\partial x_{j}}({\it\rho}u_{i}u_{j}-\unicode[STIX]{x1D61D}_{ij})+{\it\rho}\frac{\partial {\it\Phi}}{\partial x_{i}}=0.\end{eqnarray}$$

The remaining terms are

(10.48) $$\begin{eqnarray}\frac{\partial }{\partial t}\left({\it\rho}\frac{\text{D}{\it\xi}_{i}}{\text{D}t}\right)+\frac{\partial }{\partial x_{j}}\left({\it\rho}u_{j}\frac{\text{D}{\it\xi}_{i}}{\text{D}t}-\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}\right)+{\it\rho}{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+O({\it\xi}^{2})=0,\end{eqnarray}$$

or (making use of the equation of mass conservation)

(10.49) $$\begin{eqnarray}{\it\rho}\frac{\text{D}^{2}{\it\xi}_{i}}{\text{D}t^{2}}=\frac{\partial }{\partial x_{j}}\left(\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}\right)-{\it\rho}{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+O({\it\xi}^{2}).\end{eqnarray}$$

This equation, which can be extended to any order in ${\bf\xi}$ , provides the basis for a nonlinear perturbation theory for any flow in ideal MHD. In a linear theory we would neglect terms of $O({\it\xi}^{2})$ .

10.7 Notes on linear perturbations

The Lagrangian perturbation ${\rm\Delta}X$ of a quantity $X$ is the difference in the values of the quantity in the two flows for the same fluid element,

(10.50) $$\begin{eqnarray}{\rm\Delta}X=\hat{X}^{L}(\boldsymbol{a},t)-X^{L}(\boldsymbol{a},t).\end{eqnarray}$$

It follows that

(10.51a-c ) $$\begin{eqnarray}{\rm\Delta}s=0,\quad {\rm\Delta}{\it\rho}=-{\it\rho}\frac{\partial {\it\xi}_{i}}{\partial x_{i}}+O({\it\xi}^{2}),\quad {\rm\Delta}B_{i}=B_{j}\frac{\partial {\it\xi}_{i}}{\partial x_{j}}-B_{i}\frac{\partial {\it\xi}_{j}}{\partial x_{j}}+O({\it\xi}^{2}),\end{eqnarray}$$

and

(10.52) $$\begin{eqnarray}{\rm\Delta}u_{i}=\frac{\text{D}{\it\xi}_{i}}{\text{D}t}.\end{eqnarray}$$

In linear theory, $\boldsymbol{{\rm\nabla}}{\bf\xi}$ is small and terms higher than the first order are neglected. Thus

(10.53a-c ) $$\begin{eqnarray}{\rm\Delta}s=0,\quad {\rm\Delta}{\it\rho}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi},\quad {\rm\Delta}\boldsymbol{B}=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}-(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi})\boldsymbol{B}.\end{eqnarray}$$

In linear theory, ${\rm\Delta}s=0$ implies

(10.54) $$\begin{eqnarray}{\rm\Delta}p=\frac{{\it\gamma}p}{{\it\rho}}{\rm\Delta}{\it\rho}=-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}.\end{eqnarray}$$

The Eulerian perturbation ${\it\delta}X$ of a quantity $X$ is the difference in the values of the quantity in the two flows at the same point in space,

(10.55) $$\begin{eqnarray}{\it\delta}X=\hat{X}^{E}(\boldsymbol{x},t)-X^{E}(\boldsymbol{x},t).\end{eqnarray}$$

By Taylor’s theorem,

(10.56) $$\begin{eqnarray}{\rm\Delta}X={\it\delta}X+{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}X+O({\it\xi}^{2}),\end{eqnarray}$$

and so, in linear theory,

(10.57) $$\begin{eqnarray}{\it\delta}X={\rm\Delta}X-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}X.\end{eqnarray}$$

Thus

(10.58) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\delta}{\it\rho}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho},\\ {\it\delta}p=-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p,\\ {\it\delta}\boldsymbol{B}=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\xi}-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}-(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi})\boldsymbol{B},\end{array}\right\}\end{eqnarray}$$

exactly as was obtained in § 5 for perturbations of magnetostatic equilibria.

The relation

(10.59) $$\begin{eqnarray}{\it\delta}\boldsymbol{u}=\frac{\text{D}{\bf\xi}}{\text{D}t}-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}\end{eqnarray}$$

can be used to introduce the Lagrangian displacement into a linear theory derived using Eulerian perturbations. Only in the case of a static basic state, $\boldsymbol{u}=\mathbf{0}$ , does this reduce to the simple relation ${\it\delta}\boldsymbol{u}=\partial {\bf\xi}/\partial t$ .

11 Waves and instabilities in stratified rotating astrophysical bodies

11.1 The energy principle

For linear perturbations to a static equilibrium ( $\boldsymbol{u}=\mathbf{0}$ ), the displacement satisfies

(11.1) $$\begin{eqnarray}{\it\rho}\frac{\partial ^{2}{\it\xi}_{i}}{\partial t^{2}}=-{\it\rho}\frac{\partial {\it\delta}{\it\Phi}}{\partial x_{i}}-{\it\rho}{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+\frac{\partial }{\partial x_{j}}\left(\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}\right),\end{eqnarray}$$

where we now allow for self-gravitation through

(11.2) $$\begin{eqnarray}{\rm\nabla}^{2}\,{\it\delta}{\it\Phi}=4{\rm\pi}G\,{\it\delta}{\it\rho}=-4{\rm\pi}G\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}{\bf\xi}).\end{eqnarray}$$

We may write (11.1) in the form

(11.3) $$\begin{eqnarray}\frac{\partial ^{2}{\bf\xi}}{\partial t^{2}}={\mathcal{F}}{\bf\xi},\end{eqnarray}$$

where ${\mathcal{F}}$ is a linear differential operator (or integro-differential if self-gravitation is taken into account). The force operator ${\mathcal{F}}$ can be shown to be self-adjoint with respect to the inner product,

(11.4) $$\begin{eqnarray}\langle \boldsymbol{{\it\eta}},{\bf\xi}\rangle =\int {\it\rho}\boldsymbol{{\it\eta}}^{\ast }\boldsymbol{\cdot }{\bf\xi}\,\text{d}V\end{eqnarray}$$

if appropriate boundary conditions apply to the vector fields ${\bf\xi}$ and $\boldsymbol{{\it\eta}}$ . Let ${\it\delta}{\it\Psi}$ be the gravitational potential perturbation associated with the displacement $\boldsymbol{{\it\eta}}$ , so ${\rm\nabla}^{2}{\it\delta}{\it\Psi}=-4{\rm\pi}G\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{{\it\eta}})$ . Then

(11.5) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{{\it\eta}},{\mathcal{F}}{\bf\xi}\rangle & = & \displaystyle \int \left[-{\it\rho}{\it\eta}_{i}^{\ast }\frac{\partial {\it\delta}{\it\Phi}}{\partial x_{i}}-{\it\rho}{\it\eta}_{i}^{\ast }{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+{\it\eta}_{i}^{\ast }\frac{\partial }{\partial x_{j}}\left(\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}\right)\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int \left[-{\it\delta}{\it\Phi}\frac{{\rm\nabla}^{2}{\it\delta}{\it\Psi}^{\ast }}{4{\rm\pi}G}-{\it\rho}{\it\xi}_{i}{\it\eta}_{j}^{\ast }\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}-\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}\frac{\partial {\it\eta}_{i}^{\ast }}{\partial x_{j}}\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int \left[\frac{\boldsymbol{{\rm\nabla}}({\it\delta}{\it\Phi})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}({\it\delta}{\it\Psi}^{\ast })}{4{\rm\pi}G}-{\it\rho}{\it\xi}_{i}{\it\eta}_{j}^{\ast }\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+{\it\xi}_{k}\frac{\partial }{\partial x_{\ell }}\left(\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\eta}_{i}^{\ast }}{\partial x_{j}}\right)\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int \left[-{\it\delta}{\it\Psi}^{\ast }\frac{{\rm\nabla}^{2}{\it\delta}{\it\Phi}}{4{\rm\pi}G}-{\it\rho}{\it\xi}_{i}{\it\eta}_{j}^{\ast }\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+{\it\xi}_{i}\frac{\partial }{\partial x_{j}}\left(V_{k\ell ij}\frac{\partial {\it\eta}_{k}^{\ast }}{\partial x_{\ell }}\right)\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int \left[-{\it\rho}{\it\xi}_{i}\frac{\partial {\it\delta}{\it\Psi}^{\ast }}{\partial x_{i}}-{\it\rho}{\it\xi}_{i}{\it\eta}_{j}^{\ast }\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+{\it\xi}_{i}\frac{\partial }{\partial x_{j}}\left(\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\eta}_{k}^{\ast }}{\partial x_{\ell }}\right)\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \langle {\mathcal{F}}\boldsymbol{{\it\eta}},{\bf\xi}\rangle .\end{eqnarray}$$

Here the integrals are over all space. We assume that the exterior of the body is a medium of zero density in which the force-free limit of MHD holds and $\boldsymbol{B}$ decays sufficiently fast as $|\boldsymbol{x}|\rightarrow \infty$ that we may integrate freely by parts (using the divergence theorem) and ignore surface terms. We also assume that the body is isolated and self-gravitating, so that ${\it\delta}{\it\Phi}=O(r^{-1})$ , or in fact $O(r^{-2})$ if ${\it\delta}M=0$ . We have used the symmetry properties of $\partial ^{2}{\it\Phi}/\partial x_{i}\partial x_{j}$ and $\unicode[STIX]{x1D61D}_{ijk\ell }$ .

The functional

(11.6) $$\begin{eqnarray}W[{\bf\xi}]=-\frac{1}{2}\langle {\bf\xi},{\mathcal{F}}{\bf\xi}\rangle =\frac{1}{2}\int \left(-\frac{|\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}|^{2}}{4{\rm\pi}G}+{\it\rho}{\it\xi}_{i}^{\ast }{\it\xi}_{j}\frac{\partial ^{2}{\it\Phi}}{\partial x_{i}\partial x_{j}}+\unicode[STIX]{x1D61D}_{ijk\ell }\frac{\partial {\it\xi}_{i}^{\ast }}{\partial x_{j}}\frac{\partial {\it\xi}_{k}}{\partial x_{\ell }}\right)\,\text{d}V\end{eqnarray}$$

is therefore real and represents the change in potential energy associated with the displacement ${\bf\xi}$ .

If the basic state is static, we may consider normal mode solutions of the form

(11.7) $$\begin{eqnarray}{\bf\xi}=\text{Re}\left[\tilde{{\bf\xi}}(\boldsymbol{x})\exp (-\text{i}{\it\omega}t)\right],\end{eqnarray}$$

for which we obtain

(11.8) $$\begin{eqnarray}-{\it\omega}^{2}\tilde{{\bf\xi}}={\mathcal{F}}\tilde{{\bf\xi}}\end{eqnarray}$$

and

(11.9) $$\begin{eqnarray}{\it\omega}^{2}=-\frac{\langle \tilde{{\bf\xi}},{\mathcal{F}}\tilde{{\bf\xi}}\rangle }{\langle \tilde{{\bf\xi}},\tilde{{\bf\xi}}\rangle }=\frac{2W[\tilde{{\bf\xi}}]}{\langle \tilde{{\bf\xi}},\tilde{{\bf\xi}}\rangle }.\end{eqnarray}$$

Therefore ${\it\omega}^{2}$ is real and we have either oscillations ( ${\it\omega}^{2}>0$ ) or instability ( ${\it\omega}^{2}<0$ ).

The above expression for ${\it\omega}^{2}$ satisfies the usual Rayleigh–Ritz variational principle for self-adjoint eigenvalue problems. The eigenvalues ${\it\omega}^{2}$ are the stationary values of $2W[{\bf\xi}]/\langle {\bf\xi},{\bf\xi}\rangle$ among trial displacements ${\bf\xi}$ satisfying the boundary conditions. In particular, the lowest eigenvalue is the global minimum value of $2W[{\bf\xi}]/\langle {\bf\xi},{\bf\xi}\rangle$ . Therefore the equilibrium is unstable if and only if $W[{\bf\xi}]$ can be made negative by a trial displacement  ${\bf\xi}$ satisfying the boundary conditions. This is called the energy principle.

This discussion is incomplete because it assumes that the eigenfunctions form a complete set. In general a continuous spectrum, not associated with square-integrable modes, is also present. However, it can be shown that a necessary and sufficient condition for instability is that $W[{\bf\xi}]$ can be made negative as described above. Consider the equation for twice the energy of the perturbation,

(11.10) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\left(\langle \dot{{\bf\xi}},\dot{{\bf\xi}}\rangle +2W[{\bf\xi}]\right) & = & \displaystyle \langle \ddot{{\bf\xi}},\dot{{\bf\xi}}\rangle +\langle \dot{{\bf\xi}},\ddot{{\bf\xi}}\rangle -\langle \dot{{\bf\xi}},{\mathcal{F}}{\bf\xi}\rangle -\langle {\bf\xi},{\mathcal{F}}\dot{{\bf\xi}}\rangle \nonumber\\ \displaystyle & = & \displaystyle \langle {\mathcal{F}}{\bf\xi},\dot{{\bf\xi}}\rangle +\langle \dot{{\bf\xi}},{\mathcal{F}}{\bf\xi}\rangle -\langle \dot{{\bf\xi}},{\mathcal{F}}{\bf\xi}\rangle -\langle {\mathcal{F}}{\bf\xi},\dot{{\bf\xi}}\rangle \nonumber\\ \displaystyle & = & \displaystyle 0.\end{eqnarray}$$

Therefore

(11.11) $$\begin{eqnarray}\langle \dot{{\bf\xi}},\dot{{\bf\xi}}\rangle +2W[{\bf\xi}]=2E=\text{const.},\end{eqnarray}$$

where $E$ is determined by the initial data ${\bf\xi}_{0}$ and $\dot{{\bf\xi}}_{0}$ . If $W$ is positive definite then the equilibrium is stable because ${\bf\xi}$ is limited by the constraint $W[{\bf\xi}]\leqslant E$ .

Suppose that a (real) trial displacement $\boldsymbol{{\it\eta}}$ can be found for which

(11.12) $$\begin{eqnarray}\frac{2W[\boldsymbol{{\it\eta}}]}{\langle \boldsymbol{{\it\eta}},\boldsymbol{{\it\eta}}\rangle }=-{\it\gamma}^{2},\end{eqnarray}$$

where ${\it\gamma}>0$ . Then let the initial conditions be ${\bf\xi}_{0}=\boldsymbol{{\it\eta}}$ and $\dot{{\bf\xi}}_{0}={\it\gamma}\boldsymbol{{\it\eta}}$ so that

(11.13) $$\begin{eqnarray}\langle \dot{{\bf\xi}},\dot{{\bf\xi}}\rangle +2W[{\bf\xi}]=2E=0.\end{eqnarray}$$

Now let

(11.14) $$\begin{eqnarray}a(t)=\ln \left(\frac{\langle {\bf\xi},{\bf\xi}\rangle }{\langle \boldsymbol{{\it\eta}},\boldsymbol{{\it\eta}}\rangle }\right)\end{eqnarray}$$

so that

(11.15) $$\begin{eqnarray}\frac{\text{d}a}{\text{d}t}=\frac{2\langle {\bf\xi},\dot{{\bf\xi}}\rangle }{\langle {\bf\xi},{\bf\xi}\rangle }\end{eqnarray}$$

and

(11.16) $$\begin{eqnarray}\displaystyle \frac{\text{d}^{2}a}{\text{d}t^{2}} & = & \displaystyle \frac{2(\langle {\bf\xi},{\mathcal{F}}{\bf\xi}\rangle +\langle \dot{{\bf\xi}},\dot{{\bf\xi}}\rangle )\langle {\bf\xi},{\bf\xi}\rangle -4\langle {\bf\xi},\dot{{\bf\xi}}\rangle ^{2}}{\langle {\bf\xi},{\bf\xi}\rangle ^{2}}\nonumber\\ \displaystyle & = & \displaystyle \frac{2(-2W[{\bf\xi}]+\langle \dot{{\bf\xi}},\dot{{\bf\xi}}\rangle )\langle {\bf\xi},{\bf\xi}\rangle -4\langle {\bf\xi},\dot{{\bf\xi}}\rangle ^{2}}{\langle {\bf\xi},{\bf\xi}\rangle ^{2}}\nonumber\\ \displaystyle & = & \displaystyle \frac{4(\langle \dot{{\bf\xi}},\dot{{\bf\xi}}\rangle \langle {\bf\xi},{\bf\xi}\rangle -\langle {\bf\xi},\dot{{\bf\xi}}\rangle ^{2})}{\langle {\bf\xi},{\bf\xi}\rangle ^{2}}\nonumber\\ \displaystyle & {\geqslant} & \displaystyle 0\end{eqnarray}$$

by the Cauchy–Schwarz inequality. Thus

(11.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{d}a}{\text{d}t}\geqslant {\dot{a}}_{0}=2{\it\gamma}\\ \displaystyle a\geqslant 2{\it\gamma}t+a_{0}=2{\it\gamma}t.\end{array}\right\}\end{eqnarray}$$

Therefore the disturbance with these initial conditions grows at least as fast as $\exp ({\it\gamma}t)$ and the equilibrium is unstable.

11.2 Spherically symmetric star

The simplest model of a star neglects rotation and magnetic fields and assumes a spherically symmetric hydrostatic equilibrium in which ${\it\rho}(r)$ and $p(r)$ satisfy

(11.18) $$\begin{eqnarray}\frac{\text{d}p}{\text{d}r}=-{\it\rho}g,\end{eqnarray}$$

with inward radial gravitational acceleration

(11.19) $$\begin{eqnarray}g(r)=\frac{\text{d}{\it\Phi}}{\text{d}r}=\frac{G}{r^{2}}\int _{0}^{r}{\it\rho}(r^{\prime })\,4{\rm\pi}r^{\prime 2}\,\text{d}r^{\prime }.\end{eqnarray}$$

The stratification induced by gravity provides a non-uniform background for wave propagation.

In this case the linearized equation of motion is (cf. § 5)

(11.20) $$\begin{eqnarray}{\it\rho}\frac{\partial ^{2}{\bf\xi}}{\partial t^{2}}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}-{\it\delta}{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}-\boldsymbol{{\rm\nabla}}{\it\delta}p,\end{eqnarray}$$

with ${\it\delta}{\it\rho}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}{\bf\xi})$ , ${\rm\nabla}^{2}{\it\delta}{\it\Phi}=4{\rm\pi}G\,{\it\delta}{\it\rho}$ and ${\it\delta}p=-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p$ . For normal modes $\propto \exp (-\text{i}{\it\omega}t)$ ,

(11.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\rho}{\it\omega}^{2}{\bf\xi}={\it\rho}\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}+{\it\delta}{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}+\boldsymbol{{\rm\nabla}}{\it\delta}p\\ \displaystyle {\it\omega}^{2}\int _{V}{\it\rho}|{\bf\xi}|^{2}\,\text{d}V=\int _{V}{\bf\xi}^{\ast }\boldsymbol{\cdot }({\it\rho}\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}+{\it\delta}{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}+\boldsymbol{{\rm\nabla}}{\it\delta}p)\,\text{d}V,\end{array}\right\}\end{eqnarray}$$

where $V$ is the volume of the star. At the surface $S$ of the star, we assume that ${\it\rho}$ and $p$ vanish. Then ${\it\delta}p$ also vanishes on $S$ (assuming that ${\bf\xi}$ and its derivatives are bounded).

The ${\it\delta}p$ term can be integrated by parts as follows:

(11.22) $$\begin{eqnarray}\displaystyle \int _{V}{\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\delta}p\,\text{d}V & = & \displaystyle -\int _{V}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi})^{\ast }{\it\delta}p\,\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{V}\frac{1}{{\it\gamma}p}({\it\delta}p+{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p)^{\ast }{\it\delta}p\,\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{V}\left[\frac{|{\it\delta}p|^{2}}{{\it\gamma}p}+\frac{1}{{\it\gamma}p}({\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p)(-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi})\right]\,\text{d}V.\end{eqnarray}$$

The ${\it\delta}{\it\rho}$ term partially cancels with the above:

(11.23) $$\begin{eqnarray}\displaystyle \int _{V}{\bf\xi}^{\ast }\boldsymbol{\cdot }({\it\delta}{\it\rho}\boldsymbol{{\rm\nabla}}{\it\Phi}) & = & \displaystyle \int _{V}(-{\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p)\frac{{\it\delta}{\it\rho}}{{\it\rho}}\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{V}({\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p)(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}+{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln {\it\rho})\,\text{d}V.\end{eqnarray}$$

Finally, the ${\it\delta}{\it\Phi}$ term can be transformed as in § 11.1 to give

(11.24) $$\begin{eqnarray}\int _{V}{\it\rho}{\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}\,\text{d}V=-\int _{\infty }\frac{|\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}|^{2}}{4{\rm\pi}G}\,\text{d}V,\end{eqnarray}$$

where the integral on the right-hand side is over all space. Thus

(11.25) $$\begin{eqnarray}\displaystyle {\it\omega}^{2}\int _{V}{\it\rho}|{\bf\xi}|^{2}\,\text{d}V & = & \displaystyle -\int _{\infty }\frac{|\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}|^{2}}{4{\rm\pi}G}\,\text{d}V\nonumber\\ \displaystyle & & \displaystyle +\,\int _{V}\left[\frac{|{\it\delta}p|^{2}}{{\it\gamma}p}-({\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p)\boldsymbol{\cdot }\left(\frac{1}{{\it\gamma}}{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln p-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln {\it\rho}\right)\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle -\int _{\infty }\frac{|\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}|^{2}}{4{\rm\pi}G}\,\text{d}V+\int _{V}\left(\frac{|{\it\delta}p|^{2}}{{\it\gamma}p}+{\it\rho}N^{2}|{\it\xi}_{r}|^{2}\right)\,\text{d}V,\end{eqnarray}$$

where $N(r)$ is the Brunt–Väisälä frequency ¹⁷ (or buoyancy frequency) given by

(11.26) $$\begin{eqnarray}N^{2}=g\left(\frac{1}{{\it\gamma}}\frac{\text{d}\ln p}{\text{d}r}-\frac{\text{d}\ln {\it\rho}}{\text{d}r}\right)\propto g\frac{\text{d}s}{\text{d}r}.\end{eqnarray}$$

$N$ is the frequency of oscillation of a fluid element that is displaced vertically in a stably stratified atmosphere if it maintains pressure equilibrium with its surroundings. The stratification is stable if the specific entropy increases outwards.

The integral expression for ${\it\omega}^{2}$ satisfies the energy principle. There are three contributions to ${\it\omega}^{2}$ : the self-gravitational term (destabilizing), the acoustic term (stabilizing) and the buoyancy term (stabilizing if $N^{2}>0$ ).

If $N^{2}<0$ for any interval of $r$ , a trial displacement can always be found such that ${\it\omega}^{2}<0$ . This is done by localizing ${\it\xi}_{r}$ in that interval and arranging the other components of ${\bf\xi}$ such that ${\it\delta}p=0$ . Therefore the star is unstable if $\partial s/\partial r<0$ anywhere. This is Schwarzschild’s criterion ¹⁸ for convective instability.

11.3 Modes of an incompressible sphere

Note: in this subsection $(r,{\it\theta},{\it\phi})$ are spherical polar coordinates.

Analytical solutions can be obtained in the case of a homogeneous incompressible ‘star’ of mass $M$ and radius $R$ which has

(11.27) $$\begin{eqnarray}{\it\rho}=\left(\frac{3M}{4{\rm\pi}R^{3}}\right)H(R-r),\end{eqnarray}$$

where $H$ is the Heaviside step function. For $r\leqslant R$ we have

(11.28a,b ) $$\begin{eqnarray}g=\frac{GMr}{R^{3}},\quad p=\frac{3GM^{2}(R^{2}-r^{2})}{8{\rm\pi}R^{6}}.\end{eqnarray}$$

For an incompressible fluid,

(11.29) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}=0,\end{eqnarray}$$

(11.30) $$\begin{eqnarray}{\it\delta}{\it\rho}=-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}={\it\xi}_{r}\left(\frac{3M}{4{\rm\pi}R^{3}}\right){\it\delta}(r-R)\end{eqnarray}$$

and

(11.31) $$\begin{eqnarray}{\rm\nabla}^{2}{\it\delta}{\it\Phi}=4{\rm\pi}G\,{\it\delta}{\it\rho}={\it\xi}_{r}\left(\frac{3GM}{R^{3}}\right){\it\delta}(r-R),\end{eqnarray}$$

while ${\it\delta}p$ is indeterminate and is a variable independent of ${\bf\xi}$ . The linearized equation of motion is

(11.32) $$\begin{eqnarray}-{\it\rho}{\it\omega}^{2}{\bf\xi}=-{\it\rho}\boldsymbol{{\rm\nabla}}{\it\delta}{\it\Phi}-\boldsymbol{{\rm\nabla}}{\it\delta}p.\end{eqnarray}$$

Therefore we have potential flow: ${\bf\xi}=\boldsymbol{{\rm\nabla}}U$ , with ${\rm\nabla}^{2}U=0$ and $-{\it\rho}{\it\omega}^{2}U=-{\it\rho}{\it\delta}{\it\Phi}-{\it\delta}p$ in $r\leqslant R$ . Appropriate solutions of Laplace’s equation regular at $r=0$ are the solid spherical harmonics (with arbitrary normalization)

(11.33) $$\begin{eqnarray}U=r^{\ell }Y_{\ell }^{m}({\it\theta},{\it\phi}),\end{eqnarray}$$

where $\ell$ and $m$ are integers with $\ell \geqslant |m|$ . Equation (11.31) also implies

(11.34) $$\begin{eqnarray}{\it\delta}{\it\Phi}=\left\{\begin{array}{@{}ll@{}}Ar^{\ell }Y_{\ell }^{m},\quad & r<R,\\ Br^{-\ell -1}Y_{\ell }^{m},\quad & r>R,\end{array}\right.\end{eqnarray}$$

where $A$ and $B$ are constants to be determined. The matching conditions from (11.31) at $r=R$ are

(11.35a,b ) $$\begin{eqnarray}[{\it\delta}{\it\Phi}]=0,\quad \left[\frac{\partial {\it\delta}{\it\Phi}}{\partial r}\right]={\it\xi}_{r}\left(\frac{3GM}{R^{3}}\right).\end{eqnarray}$$

Thus

(11.36a,b ) $$\begin{eqnarray}BR^{-\ell -1}-AR^{\ell }=0,\quad -(\ell +1)BR^{-\ell -2}-\ell AR^{\ell -1}=\ell R^{\ell -1}\left(\frac{3GM}{R^{3}}\right),\end{eqnarray}$$

with solution

(11.37a,b ) $$\begin{eqnarray}A=-\frac{\ell }{2\ell +1}\left(\frac{3GM}{R^{3}}\right),\quad B=AR^{2\ell +1}.\end{eqnarray}$$

At $r=R$ the Lagrangian pressure perturbation should vanish:

(11.38) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\rm\Delta}p={\it\delta}p+{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p=0\\ \displaystyle \left(\frac{3M}{4{\rm\pi}R^{3}}\right)\left[{\it\omega}^{2}R^{\ell }+\left(\frac{\ell }{2\ell +1}\right)\left(\frac{3GM}{R^{3}}\right)R^{\ell }\right]-\frac{3GM^{2}}{4{\rm\pi}R^{5}}\ell R^{\ell -1}=0\\ \displaystyle {\it\omega}^{2}=\left(\ell -\frac{3\ell }{2\ell +1}\right)\frac{GM}{R^{2}}=\frac{2\ell (\ell -1)}{2\ell +1}\frac{GM}{R^{3}}.\end{array}\right\}\end{eqnarray}$$

This result was obtained by Lord Kelvin. Since ${\it\omega}^{2}\geqslant 0$ the star is stable. Note that $\ell =0$ corresponds to ${\bf\xi}=\mathbf{0}$ and $\ell =1$ corresponds to ${\bf\xi}=\text{const.}$ , which is a translational mode of zero frequency. The remaining modes are non-trivial and are called $f$  modes (fundamental modes). These can be thought of as surface gravity waves, related to ocean waves for which ${\it\omega}^{2}=gk$ . In the first expression for ${\it\omega}^{2}$ above, the first term in brackets derives from surface gravity, while the second derives from self-gravity.

11.4 The plane-parallel atmosphere

The local dynamics of a stellar atmosphere can be studied in a Cartesian (‘plane-parallel’) approximation. The gravitational acceleration is taken to be constant (appropriate to an atmosphere) and in the $-z$ direction. For hydrostatic equilibrium,

(11.39) $$\begin{eqnarray}\frac{\text{d}p}{\text{d}z}=-{\it\rho}g.\end{eqnarray}$$

A simple example is an isothermal atmosphere in which $p=c_{s}^{2}{\it\rho}$ with $c_{s}=\text{const.}$ :

(11.40a,b ) $$\begin{eqnarray}{\it\rho}={\it\rho}_{0}\,\text{e}^{-z/H},\quad p=p_{0}\,\text{e}^{-z/H}.\end{eqnarray}$$

$H=c_{s}^{2}/g$ is the isothermal scale height. The Brunt–Väisälä frequency in an isothermal atmosphere is given by

(11.41) $$\begin{eqnarray}N^{2}=g\left(\frac{1}{{\it\gamma}}\frac{\text{d}\ln p}{\text{d}z}-\frac{\text{d}\ln {\it\rho}}{\text{d}z}\right)=\left(1-\frac{1}{{\it\gamma}}\right)\frac{g}{H},\end{eqnarray}$$

which is constant and is positive for ${\it\gamma}>1$ . An isothermal atmosphere is stably (subadiabatically) stratified if ${\it\gamma}>1$ and neutrally (adiabatically) stratified if ${\it\gamma}=1$ .

A further example is a polytropic atmosphere in which $p\propto {\it\rho}^{1+1/m}$ in the undisturbed state, where $m$ is a positive constant. In general $1+1/m$ differs from the adiabatic exponent ${\it\gamma}$ of the gas. For hydrostatic equilibrium,

(11.42) $$\begin{eqnarray}{\it\rho}^{1/m}\frac{\text{d}{\it\rho}}{\text{d}z}\propto -{\it\rho}g\quad \Rightarrow \quad {\it\rho}^{1/m}\propto -z,\end{eqnarray}$$

if the top of the atmosphere is located at $z=0$ , with vacuum above. Let

(11.43) $$\begin{eqnarray}{\it\rho}={\it\rho}_{0}\left(-\frac{z}{H}\right)^{m},\end{eqnarray}$$

for $z<0$ , where ${\it\rho}_{0}$ and $H$ are constants. Then

(11.44) $$\begin{eqnarray}p=p_{0}\left(-\frac{z}{H}\right)^{m+1},\end{eqnarray}$$

where

(11.45) $$\begin{eqnarray}p_{0}=\frac{{\it\rho}_{0}gH}{m+1}\end{eqnarray}$$

to satisfy $\text{d}p/\text{d}z=-{\it\rho}g$ . In this case

(11.46) $$\begin{eqnarray}N^{2}=\left(m-\frac{m+1}{{\it\gamma}}\right)\frac{g}{-z}.\end{eqnarray}$$

We return to the linearized equations, looking for solutions of the form

(11.47) $$\begin{eqnarray}{\bf\xi}=\text{Re}\left[\tilde{{\bf\xi}}(z)\exp (-\text{i}{\it\omega}t+\text{i}\boldsymbol{k}_{h}\boldsymbol{\cdot }\boldsymbol{x})\right],\quad \text{etc.},\end{eqnarray}$$

where ‘ $h$ ’ denotes a horizontal vector (having only $x$ and $y$ components). Then

(11.48) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle -{\it\rho}{\it\omega}^{2}{\bf\xi}_{h}=-\text{i}\boldsymbol{k}_{h}\,{\it\delta}p,\\ \displaystyle -{\it\rho}{\it\omega}^{2}{\it\xi}_{z}=-g\,{\it\delta}{\it\rho}-\frac{\text{d}\,{\it\delta}p}{\text{d}z},\\ \displaystyle {\it\delta}{\it\rho}=-{\it\xi}_{z}\frac{\text{d}{\it\rho}}{\text{d}z}-{\it\rho}{\it\Delta},\\ \displaystyle {\it\delta}p=-{\it\xi}_{z}\frac{\text{d}p}{\text{d}z}-{\it\gamma}p{\it\Delta},\end{array}\right\}\end{eqnarray}$$

where

(11.49) $$\begin{eqnarray}{\it\Delta}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}=\text{i}\boldsymbol{k}_{h}\boldsymbol{\cdot }{\bf\xi}_{h}+\frac{\text{d}{\it\xi}_{z}}{\text{d}z}.\end{eqnarray}$$

The self-gravitation of the perturbation is neglected in the atmosphere: ${\it\delta}{\it\Phi}=0$ (the Cowling approximation). Note that only two $z$ -derivatives of perturbation quantities occur: $\text{d}\,{\it\delta}p/\text{d}z$ and $\text{d}{\it\xi}_{z}/\text{d}z$ . This is a second-order system of ordinary differential equations (ODEs), combined with algebraic equations.

We can easily eliminate ${\bf\xi}_{h}$ to obtain

(11.50) $$\begin{eqnarray}{\it\Delta}=-\frac{k_{h}^{2}}{{\it\rho}{\it\omega}^{2}}\,{\it\delta}p+\frac{\text{d}{\it\xi}_{z}}{\text{d}z},\end{eqnarray}$$

where $k_{h}=|\boldsymbol{k}_{h}|$ , and eliminate ${\it\delta}{\it\rho}$ to obtain

(11.51) $$\begin{eqnarray}-{\it\rho}{\it\omega}^{2}{\it\xi}_{z}=g{\it\xi}_{z}\frac{\text{d}{\it\rho}}{\text{d}z}+{\it\rho}g{\it\Delta}-\frac{\text{d}\,{\it\delta}p}{\text{d}z}.\end{eqnarray}$$

We consider these two differential equations in combination with the remaining algebraic equation

(11.52) $$\begin{eqnarray}{\it\delta}p={\it\rho}g{\it\xi}_{z}-{\it\gamma}p{\it\Delta}.\end{eqnarray}$$

A first approach is to solve the algebraic equation for ${\it\Delta}$ and substitute to obtain the two coupled ODEs

(11.53) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{d}{\it\xi}_{z}}{\text{d}z}=\frac{g}{v_{s}^{2}}{\it\xi}_{z}+\frac{1}{{\it\rho}v_{s}^{2}}\left(\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}-1\right){\it\delta}p,\\ \displaystyle \frac{\text{d}\,{\it\delta}p}{\text{d}z}={\it\rho}({\it\omega}^{2}-N^{2}){\it\xi}_{z}-\frac{g}{v_{s}^{2}}{\it\delta}p.\end{array}\right\}\end{eqnarray}$$

Note that $v_{s}^{2}k_{h}^{2}$ is the square of the ‘Lamb frequency’, i.e. the ( $z$ -dependent) frequency of a horizontal sound wave of wavenumber $k_{h}$ . In a short-wavelength (WKB) approximation, where ${\it\xi}_{z}\propto \exp \left[\text{i}\int k_{z}(z)\,\text{d}z\right]$ with $k_{z}\gg g/v_{s}^{2}$ , the local dispersion relation derived from these ODEs is

(11.54) $$\begin{eqnarray}v_{s}^{2}k_{z}^{2}=({\it\omega}^{2}-N^{2})\left(1-\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}\right).\end{eqnarray}$$

Propagating waves ( $k_{z}^{2}>0$ ) are possible when

(11.55a,b ) $$\begin{eqnarray}\text{either }{\it\omega}^{2}>\max (v_{s}^{2}k_{h}^{2},N^{2})\quad \text{or}\quad 0<{\it\omega}^{2}<\min (v_{s}^{2}k_{h}^{2},N^{2}).\end{eqnarray}$$

The high-frequency branch describes $p$ modes (acoustic waves: ‘ $p$ ’ for pressure) while the low-frequency branch describes $g$ modes (internal gravity waves: ‘ $g$ ’ for gravity).

There is a special incompressible solution in which ${\it\Delta}=0$ , i.e. ${\it\delta}p={\it\rho}g{\it\xi}_{z}$ . This satisfies

(11.56a,b ) $$\begin{eqnarray}\frac{\text{d}{\it\xi}_{z}}{\text{d}z}=\frac{gk_{h}^{2}}{{\it\omega}^{2}}{\it\xi}_{z},\quad \frac{\text{d}{\it\xi}_{z}}{\text{d}z}=\frac{{\it\omega}^{2}}{g}{\it\xi}_{z}.\end{eqnarray}$$

For compatibility of these equations,

(11.57) $$\begin{eqnarray}\frac{gk_{h}^{2}}{{\it\omega}^{2}}=\frac{{\it\omega}^{2}}{g}\quad \Rightarrow \quad {\it\omega}^{2}=\pm gk_{h}.\end{eqnarray}$$

The acceptable solution in which ${\it\xi}_{z}$ decays with depth is

(11.58a,b ) $$\begin{eqnarray}{\it\omega}^{2}=gk_{h},\quad {\it\xi}_{z}\propto \exp (k_{h}z).\end{eqnarray}$$

This is a surface gravity wave known in stellar oscillations as the $f$ mode (fundamental mode). It is vertically evanescent.

The other wave solutions ( $p$ and $g$ modes) can be found analytically in the case of a polytropic atmosphere. ¹⁹ We now eliminate variables in favour of ${\it\Delta}$ . First use the algebraic relation to eliminate ${\it\delta}p$ :

(11.59) $$\begin{eqnarray}{\it\Delta}=-\frac{gk_{h}^{2}}{{\it\omega}^{2}}{\it\xi}_{z}+\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}{\it\Delta}+\frac{\text{d}{\it\xi}_{z}}{\text{d}z},\end{eqnarray}$$

(11.60) $$\begin{eqnarray}-{\it\rho}{\it\omega}^{2}{\it\xi}_{z}={\it\rho}g{\it\Delta}-{\it\rho}g\frac{\text{d}{\it\xi}_{z}}{\text{d}z}+\frac{\text{d}({\it\gamma}p{\it\Delta})}{\text{d}z}.\end{eqnarray}$$

Then eliminate $\text{d}{\it\xi}_{z}/\text{d}z$ :

(11.61) $$\begin{eqnarray}-{\it\rho}{\it\omega}^{2}{\it\xi}_{z}=-{\it\rho}g\left(\frac{gk_{h}^{2}}{{\it\omega}^{2}}{\it\xi}_{z}-\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}{\it\Delta}\right)+\frac{\text{d}({\it\gamma}p{\it\Delta})}{\text{d}z}.\end{eqnarray}$$

Thus we have

(11.62) $$\begin{eqnarray}\frac{\text{d}{\it\xi}_{z}}{\text{d}z}-\frac{gk_{h}^{2}}{{\it\omega}^{2}}{\it\xi}_{z}=\left(1-\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}\right){\it\Delta},\end{eqnarray}$$

(11.63) $$\begin{eqnarray}-{\it\omega}^{2}\left(1-\frac{g^{2}k_{h}^{2}}{{\it\omega}^{4}}\right){\it\xi}_{z}=\frac{1}{{\it\rho}}\frac{\text{d}({\it\rho}v_{s}^{2}{\it\Delta})}{\text{d}z}+v_{s}^{2}\frac{gk_{h}^{2}}{{\it\omega}^{2}}{\it\Delta}.\end{eqnarray}$$

Combine, eliminating ${\it\xi}_{z}$ :

(11.64) $$\begin{eqnarray}\left(\frac{\text{d}}{\text{d}z}-\frac{gk_{h}^{2}}{{\it\omega}^{2}}\right)\left[v_{s}^{2}\frac{\text{d}{\it\Delta}}{\text{d}z}+\frac{1}{{\it\rho}}\frac{\text{d}({\it\rho}v_{s}^{2})}{\text{d}z}{\it\Delta}+v_{s}^{2}\frac{gk_{h}^{2}}{{\it\omega}^{2}}{\it\Delta}\right]+{\it\omega}^{2}\left(1-\frac{g^{2}k_{h}^{2}}{{\it\omega}^{4}}\right)\left(1-\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}\right){\it\Delta}=0.\end{eqnarray}$$

Expand out:

(11.65) $$\begin{eqnarray}\displaystyle & & \displaystyle v_{s}^{2}\frac{\text{d}^{2}{\it\Delta}}{\text{d}z^{2}}+\left[\frac{\text{d}v_{s}^{2}}{\text{d}z}+\frac{1}{{\it\rho}}\frac{\text{d}({\it\rho}v_{s}^{2})}{\text{d}z}\right]\frac{\text{d}{\it\Delta}}{\text{d}z}+\bigg[\frac{\text{d}}{\text{d}z}\left(\frac{1}{{\it\rho}}\frac{\text{d}({\it\rho}v_{s}^{2})}{\text{d}z}\right)+\frac{\text{d}v_{s}^{2}}{\text{d}z}\frac{gk_{h}^{2}}{{\it\omega}^{2}}-\frac{gk_{h}^{2}}{{\it\omega}^{2}}\frac{1}{{\it\rho}}\frac{\text{d}({\it\rho}v_{s}^{2})}{\text{d}z}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\it\omega}^{2}\left(1-\frac{g^{2}k_{h}^{2}}{{\it\omega}^{4}}-\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}\right)\bigg]{\it\Delta}=0.\end{eqnarray}$$

In the case of a polytropic atmosphere, $v_{s}^{2}\propto z$ and ${\it\rho}v_{s}^{2}\propto z^{m+1}$ :

(11.66) $$\begin{eqnarray}v_{s}^{2}\frac{\text{d}^{2}{\it\Delta}}{\text{d}z^{2}}+(m+2)\frac{v_{s}^{2}}{z}\frac{\text{d}{\it\Delta}}{\text{d}z}+\left[-m\frac{gk_{h}^{2}}{{\it\omega}^{2}}\frac{v_{s}^{2}}{z}+{\it\omega}^{2}\left(1-\frac{g^{2}k_{h}^{2}}{{\it\omega}^{4}}-\frac{v_{s}^{2}k_{h}^{2}}{{\it\omega}^{2}}\right)\right]{\it\Delta}=0.\end{eqnarray}$$

In fact, $v_{s}^{2}/z=-{\it\gamma}g/(m+1)$ . Divide through by this factor:

(11.67) $$\begin{eqnarray}z\frac{\text{d}^{2}{\it\Delta}}{\text{d}z^{2}}+(m+2)\frac{\text{d}{\it\Delta}}{\text{d}z}-\left[m\frac{gk_{h}^{2}}{{\it\omega}^{2}}+\frac{(m+1)}{{\it\gamma}g}{\it\omega}^{2}\left(1-\frac{g^{2}k_{h}^{2}}{{\it\omega}^{4}}\right)+k_{h}^{2}z\right]{\it\Delta}=0.\end{eqnarray}$$

Finally,

(11.68) $$\begin{eqnarray}z\frac{\text{d}^{2}{\it\Delta}}{\text{d}z^{2}}+(m+2)\frac{\text{d}{\it\Delta}}{\text{d}z}-(A+k_{h}z)k_{h}{\it\Delta}=0,\end{eqnarray}$$

where

(11.69) $$\begin{eqnarray}A=\frac{(m+1)}{{\it\gamma}}\frac{{\it\omega}^{2}}{gk_{h}}+\left(m-\frac{m+1}{{\it\gamma}}\right)\frac{gk_{h}}{{\it\omega}^{2}}\end{eqnarray}$$

is a dimensionless constant. Let ${\it\Delta}=w(z)\,\text{e}^{k_{h}z}$ :

(11.70) $$\begin{eqnarray}z\frac{\text{d}^{2}w}{\text{d}z^{2}}+(m+2+2k_{h}z)\frac{\text{d}w}{\text{d}z}-(A-m-2)k_{h}w=0.\end{eqnarray}$$

This is related to the confluent hypergeometric equation and has a regular singular point at $z=0$ . Using the method of Frobenius, we seek power-series solutions

(11.71) $$\begin{eqnarray}w=\mathop{\sum }_{r=0}^{\infty }a_{r}z^{{\it\sigma}+r},\end{eqnarray}$$

where ${\it\sigma}$ is a number to be determined and $a_{0}\neq 0$ . The indicial equation is

(11.72) $$\begin{eqnarray}{\it\sigma}({\it\sigma}+m+1)=0\end{eqnarray}$$

and the regular solution has ${\it\sigma}=0$ . The recurrence relation is then

(11.73) $$\begin{eqnarray}\frac{a_{r+1}}{a_{r}}=\frac{(A-m-2-2r)k_{h}}{(r+1)(r+m+2)}.\end{eqnarray}$$

In the case of an infinite series, $a_{r+1}/a_{r}\sim -2k_{h}/r$ as $r\rightarrow \infty$ , so $w$ behaves like $\text{e}^{-2k_{h}z}$ and ${\it\Delta}$ diverges like $\text{e}^{-k_{h}z}$ as $z\rightarrow -\infty$ . Solutions in which ${\it\Delta}$ decays with depth are those for which the series terminates and $w$ is a polynomial. For a polynomial of degree $n-1$ ( $n\geqslant 1$ ),

(11.74) $$\begin{eqnarray}A=2n+m.\end{eqnarray}$$

This gives a quadratic equation for ${\it\omega}^{2}$ :

(11.75) $$\begin{eqnarray}\frac{(m+1)}{{\it\gamma}}\left(\frac{{\it\omega}^{2}}{gk_{h}}\right)^{2}-(2n+m)\left(\frac{{\it\omega}^{2}}{gk_{h}}\right)+\left(m-\frac{m+1}{{\it\gamma}}\right)=0.\end{eqnarray}$$

A negative root for ${\it\omega}^{2}$ exists if and only if $m-(m+1)/{\it\gamma}<0$ , i.e.  $N^{2}<0$ , as expected from Schwarzschild’s criterion for stability.

Figure 19. Dispersion relation, in arbitrary units, for a stably stratified plane-parallel polytropic atmosphere with $m=3$ and ${\it\gamma}=5/3$ . The dashed line is the $f$  mode. Above it are the first ten $p$  modes and below it are the first ten $g$  modes. Each curve is a parabola.

For $n\gg 1$ , the large root is

(11.76) $$\begin{eqnarray}\frac{{\it\omega}^{2}}{gk_{h}}\sim \frac{2n{\it\gamma}}{m+1}\quad (p~\text{modes},{\it\omega}^{2}\propto v_{s}^{2})\end{eqnarray}$$

and the small root is

(11.77) $$\begin{eqnarray}\frac{{\it\omega}^{2}}{gk_{h}}\sim \frac{1}{2n}\left(m-\frac{m+1}{{\it\gamma}}\right)\quad (g~\text{modes},{\it\omega}^{2}\propto N^{2}).\end{eqnarray}$$

The $f$  mode is the ‘trivial’ solution ${\it\Delta}=0$ . $p$  modes (‘ $p$ ’ for pressure) are acoustic waves, which rely on compressibility. $g$  modes are gravity waves, which rely on buoyancy. Typical branches of the dispersion relation are illustrated in figure 19.

In solar-type stars (see figures 20 and 21) the inner part (radiative zone) is convectively stable ( $N^{2}>0$ ) and the outer part (convective zone) is unstable ( $N^{2}<0$ ). However, the convection is so efficient that only a very small entropy gradient is required to sustain the convective heat flux, so $N^{2}$ is very small and negative in the convective zone. Although $g$  modes propagate in the radiative zone at frequencies smaller than $N$ , they cannot reach the surface. Only $f$ and $p$  modes (excited by convection) are observed at the solar surface.

In more massive stars the situation is reversed. Then $f$ , $p$ and $g$  modes can be observed, in principle, at the surface. $g$  modes are particularly well observed in certain classes of white dwarf.

Figure 20. A standard model of the present Sun, up to the photosphere. Density, temperature, gravity and squared buoyancy frequency are plotted versus fractional radius.

Figure 21. The same model plotted on a logarithmic scale. In the convective region where $N^{2}<0$ , the dotted line shows $-N^{2}$ instead.

Related examples: A.22–A.25.

11.5 Tidally forced oscillations

When astrophysical fluid bodies such as stars and planets orbit sufficiently close to one another, they deform each other in ways that can cause irreversible evolution of their spin and orbital motion over astronomical time scales. We consider here some of the simplest aspects of this problem.

Consider a binary star (or star–planet or planet–moon system, etc.) with a circular orbit. Let the orbital separation be $a$ and the orbital (angular) frequency

(11.78) $$\begin{eqnarray}{\it\Omega}_{o}=\left(\frac{GM}{a^{3}}\right)^{1/2},\end{eqnarray}$$

where $M=M_{1}+M_{2}$ is the combined mass of the two bodies. Let $\boldsymbol{R}_{1}(t)$ and $\boldsymbol{R}_{2}(t)$ be the position vectors of the centres of mass of the two bodies, and $\boldsymbol{d}=\boldsymbol{R}_{2}-\boldsymbol{R}_{1}$ their separation.

The gravitational potential due to body 2 (treated as a point mass or spherical mass) at position $\boldsymbol{R}_{1}+\boldsymbol{x}$ within body 1 is

(11.79) $$\begin{eqnarray}\displaystyle -\frac{GM_{2}}{|\boldsymbol{d}-\boldsymbol{x}|}\! & = & \displaystyle \!-GM_{2}\left(|\boldsymbol{d}|^{2}-2\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x}+|\boldsymbol{x}|^{2}\right)^{-1/2}\nonumber\\ \displaystyle & = & \displaystyle -\frac{GM_{2}}{|\boldsymbol{d}|}\left(1-\frac{2\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x}}{|\boldsymbol{d}|^{2}}+\frac{|\boldsymbol{x}|^{2}}{|\boldsymbol{d}|^{2}}\right)^{-1/2}\nonumber\\ \displaystyle \! & = & \displaystyle \!-\frac{GM_{2}}{|\boldsymbol{d}|}\left[1-\frac{1}{2}\left(-\frac{2\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x}}{|\boldsymbol{d}|^{2}}+\frac{|\boldsymbol{x}|^{2}}{|\boldsymbol{d}|^{2}}\right)+\frac{3}{8}\left(-\frac{2\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x}}{|\boldsymbol{d}|^{2}}+\frac{|\boldsymbol{x}|^{2}}{|\boldsymbol{d}|^{2}}\right)^{2}+\cdots \,\right]\nonumber\\ \displaystyle \! & = & \displaystyle \!-\frac{GM_{2}}{|\boldsymbol{d}|}\left[1+\frac{\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x}}{|\boldsymbol{d}|^{2}}+\frac{3(\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x})^{2}-|\boldsymbol{d}|^{2}|\boldsymbol{x}|^{2}}{2|\boldsymbol{d}|^{4}}+O\left(\frac{|\boldsymbol{x}|^{3}}{|\boldsymbol{d}|^{3}}\right)\right].\end{eqnarray}$$

In this Taylor expansion, the term independent of $\boldsymbol{x}$ is a uniform potential that has no effect. The term linear in $\boldsymbol{x}$ gives rise to a uniform acceleration $GM_{2}\boldsymbol{d}/|\boldsymbol{d}|^{3}$ , which causes the orbital motion of body 1. The remaining terms constitute the tidal potential ${\it\Psi}$ ; the quadratic terms written here are the tidal potential in the quadrupolar approximation.

For a circular orbit, the coordinate system can be chosen such that

(11.80) $$\begin{eqnarray}\boldsymbol{d}=(a\cos {\it\Omega}_{o}t,a\sin {\it\Omega}_{o}t,0).\end{eqnarray}$$

Introduce spherical polar coordinates within body 1 such that

(11.81) $$\begin{eqnarray}\boldsymbol{x}=(r\sin {\it\theta}\cos {\it\phi},r\sin {\it\theta}\sin {\it\phi},r\cos {\it\phi}).\end{eqnarray}$$

Then

(11.82) $$\begin{eqnarray}\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{x}=ar\sin {\it\theta}\cos ({\it\phi}-{\it\Omega}_{o}t),\end{eqnarray}$$

(11.83) $$\begin{eqnarray}\displaystyle {\it\Psi}\! & = & \displaystyle \!\frac{GM_{2}r^{2}}{2a^{3}}\left[1-3\sin ^{2}{\it\theta}\cos ^{2}({\it\phi}-{\it\Omega}_{o}t)\right]\nonumber\\ \displaystyle \! & = & \displaystyle \!\frac{GM_{2}r^{2}}{4a^{3}}[2-3\sin ^{2}{\it\theta}-3\sin ^{2}{\it\theta}\cos (2{\it\phi}-2{\it\Omega}_{o}t)].\end{eqnarray}$$

The first two terms are static; the remaining oscillatory part can be written as

(11.84) $$\begin{eqnarray}\text{Re}\left[-\frac{3GM_{2}r^{2}\sin ^{2}{\it\theta}}{4a^{3}}\,\text{e}^{2\text{i}({\it\phi}-{\it\Omega}_{o}t)}\right],\end{eqnarray}$$

which involves the spherical harmonic function $Y_{2}^{2}({\it\theta},{\it\phi})\propto \sin ^{2}{\it\theta}\,\text{e}^{2\text{i}{\it\phi}}$ .

The tidal frequency in a non-rotating frame is $2{\it\Omega}_{o}$ . In a frame rotating with the spin angular velocity ${\it\Omega}_{s}$ of body 1, the tidal frequency is $2({\it\Omega}_{o}-{\it\Omega}_{s})$ , owing to an angular Doppler shift.

If the tidal frequency is sufficiently small, it might be assumed that body 1 responds hydrostatically to the tidal potential. Under this assumption of an equilibrium tide, body 1 is deformed into a spheroid with a tidal bulge that points instantaneously towards body 2, and no tidal torque is exerted.

We can allow for a more general linear response, including dissipation and wave-like disturbances (a dynamical tide) as follows. The most important aspect of the tidally deformed body is its exterior gravitational potential perturbation ${\it\delta}{\it\Phi}$ , because it is only through gravity that the bodies communicate and exchange energy and angular momentum. We write the linear response as

(11.85) $$\begin{eqnarray}{\it\delta}{\it\Phi}=\text{Re}\left[-k\frac{3GM_{2}\sin ^{2}{\it\theta}}{4a^{3}}\frac{R_{1}^{5}}{r^{3}}\,\text{e}^{2\text{i}({\it\phi}-{\it\Omega}_{o}t)}\right],\end{eqnarray}$$

where $R_{1}$ is the radius of body 1 (or some appropriate measure of its radius if it is deformed by its rotation) and $k$ is the potential Love number, a dimensionless complex number that describes the amplitude and phase of the tidal response. Note that ${\it\delta}{\it\Phi}$ involves the same frequency and the same spherical harmonic $Y_{2}^{2}({\it\theta},{\it\phi})$ , but combined with $r^{-3}$ rather than $r^{2}$ to make it a valid solution of Laplace’s equation in the exterior of body 1. The factor of $R_{1}^{5}$ is introduced so that $k$ is dimensionless and measures the ratio of ${\it\delta}{\it\Phi}$ and ${\it\Psi}$ at the surface of body 1.

The imaginary part of $k$ determines the part of the tidal response that is out of phase with the tidal forcing, and which is associated with dissipation and irreversible evolution. The torque acting on the orbit of body 2 is

(11.86) $$\begin{eqnarray}\displaystyle -T\! & = & \displaystyle \!M_{2}r\sin {\it\theta}\left(-\frac{1}{r\sin {\it\theta}}\frac{\partial {\it\delta}{\it\Phi}}{\partial {\it\phi}}\right)\Bigg|_{r=a,\,{\it\theta}={\rm\pi}/2,\,{\it\phi}={\it\Omega}_{o}t}\nonumber\\ \displaystyle \! & = & \displaystyle \!\text{Re}\left(M_{2}k\frac{3GM_{2}}{4a^{3}}\frac{R_{1}^{5}}{a^{3}}2\text{i}\right)\nonumber\\ \displaystyle \! & = & \displaystyle \!-\text{Im}(k)\frac{3GM_{2}^{2}R_{1}^{5}}{2a^{6}}.\end{eqnarray}$$

By Newton’s third law, there is an equal and opposite torque, $+T$ , acting on the spin of body 1.

The orbital angular momentum about the centre of mass is

(11.87) $$\begin{eqnarray}L_{o}={\it\mu}(GMa)^{1/2},\end{eqnarray}$$

where ${\it\mu}=M_{1}M_{2}/M$ is the reduced mass of the system. This result can be obtained by considering

(11.88) $$\begin{eqnarray}M_{1}\left(\frac{M_{2}a}{M}\right)^{2}{\it\Omega}_{o}+M_{2}\left(\frac{M_{1}a}{M}\right)^{2}{\it\Omega}_{o}=\frac{M_{1}M_{2}}{M}(GMa)^{1/2}\end{eqnarray}$$

(see figure 22). It evolves according to

(11.89) $$\begin{eqnarray}\frac{\text{d}L_{o}}{\text{d}t}=-T,\end{eqnarray}$$

which determines the rate of orbital migration:

(11.90) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{1}{2}\frac{M_{1}M_{2}}{M}(GMa)^{1/2}\frac{1}{a}\frac{\text{d}a}{\text{d}t}=-\text{Im}(k)\frac{3GM_{2}^{2}R_{1}^{5}}{2a^{6}}\\ \displaystyle -\frac{1}{a}\frac{\text{d}a}{\text{d}t}=3\,\text{Im}(k)\frac{M_{2}}{M_{1}}\left(\frac{R_{1}}{a}\right)^{5}{\it\Omega}_{o}.\end{array}\right\}\end{eqnarray}$$

The orbital energy

(11.91) $$\begin{eqnarray}E_{o}=-{\it\mu}\frac{GM}{2a}\end{eqnarray}$$

evolves according to

(11.92) $$\begin{eqnarray}\frac{\text{d}E_{o}}{\text{d}t}={\it\mu}\frac{GM}{2a}\frac{1}{a}\frac{\text{d}a}{\text{d}t}=-{\it\Omega}_{o}T.\end{eqnarray}$$

The spin angular momentum $L_{s}=I_{1}{\it\Omega}_{s}$ and spin energy $E_{s}=(I_{1}{\it\Omega}_{s}^{2})/2$ , where $I_{1}$ is the moment of inertia of body 1, evolve according to

(11.93a,b ) $$\begin{eqnarray}\frac{\text{d}L_{s}}{\text{d}t}=T,\quad \frac{\text{d}E_{s}}{\text{d}t}={\it\Omega}_{s}T.\end{eqnarray}$$

The total energy therefore satisfies

(11.94) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}(E_{o}+E_{s})=({\it\Omega}_{s}-{\it\Omega}_{o})T=-D,\end{eqnarray}$$

where $D>0$ is the rate of dissipation of energy. To ensure $D>0$ , the sign of $\text{Im}(k)$ should be the same as the sign of the tidal frequency $2({\it\Omega}_{o}-{\it\Omega}_{s})$ .

Figure 22. A binary star with two components in circular orbital motion about the centre of mass.

In a dissipative spin–orbit coupling, the tidal torque $T$ tries to bring about an equalization of the spin and orbital angular velocities. Its action, mediated by gravity, is comparable to a frictional interaction between differential rotating components in a mechanical system.

In binary stars, and other cases in which the spin angular momentum is small compared to the orbital angular momentum, there is indeed a tendency towards synchronization of the spin with the orbital motion (as the Moon is synchronized with its orbit around the Earth). However, in systems of extreme mass ratio in which the spin of the large body contains most of the angular momentum, the tidal torque instead causes orbital migration away from the synchronous orbit at which ${\it\Omega}={\it\Omega}_{s}$ (figure 23).

Figure 23. Orbital migration away from the synchronous orbit driven by tidal dissipation in a system of extreme mass ratio.

This situation applies to the moons of solar-system planets, most of which migrate outwards, and to extrasolar planets in close orbits around their host stars, where the migration is usually inward and may lead to the destruction of the planet.

11.6 Rotating fluid bodies

Note: in this subsection $(r,{\it\phi},z)$ are cylindrical polar coordinates.

11.6.1 Equilibrium

The equations of ideal gas dynamics in cylindrical polar coordinates are

(11.95) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{D}u_{r}}{\text{D}t}-\frac{u_{{\it\phi}}^{2}}{r}=-\frac{\partial {\it\Phi}}{\partial r}-\frac{1}{{\it\rho}}\frac{\partial p}{\partial r},\\ \displaystyle \frac{\text{D}u_{{\it\phi}}}{\text{D}t}+\frac{u_{r}u_{{\it\phi}}}{r}=-\frac{1}{r}\frac{\partial {\it\Phi}}{\partial {\it\phi}}-\frac{1}{{\it\rho}r}\frac{\partial p}{\partial {\it\phi}},\\ \displaystyle \frac{\text{D}u_{z}}{\text{D}t}=-\frac{\partial {\it\Phi}}{\partial z}-\frac{1}{{\it\rho}}\frac{\partial p}{\partial z},\\ \displaystyle \frac{\text{D}{\it\rho}}{\text{D}t}=-{\it\rho}\left[\frac{1}{r}\frac{\partial }{\partial r}(ru_{r})+\frac{1}{r}\frac{\partial u_{{\it\phi}}}{\partial {\it\phi}}+\frac{\partial u_{z}}{\partial z}\right],\\ \displaystyle \frac{\text{D}p}{\text{D}t}=-{\it\gamma}p\left[\frac{1}{r}\frac{\partial }{\partial r}(ru_{r})+\frac{1}{r}\frac{\partial u_{{\it\phi}}}{\partial {\it\phi}}+\frac{\partial u_{z}}{\partial z}\right],\end{array}\right\}\end{eqnarray}$$

with

(11.96) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}=\frac{\partial }{\partial t}+u_{r}\frac{\partial }{\partial r}+\frac{u_{{\it\phi}}}{r}\frac{\partial }{\partial {\it\phi}}+u_{z}\frac{\partial }{\partial z}.\end{eqnarray}$$

Consider a steady, axisymmetric basic state with density ${\it\rho}(r,z)$ , pressure $p(r,z)$ , gravitational potential ${\it\Phi}(r,z)$ and with differential rotation

(11.97) $$\begin{eqnarray}\boldsymbol{u}=r{\it\Omega}(r,z)\,\boldsymbol{e}_{{\it\phi}}.\end{eqnarray}$$

For equilibrium we require

(11.98) $$\begin{eqnarray}-r{\it\Omega}^{2}\,\boldsymbol{e}_{r}=-\boldsymbol{{\rm\nabla}}{\it\Phi}-\frac{1}{{\it\rho}}\boldsymbol{{\rm\nabla}}p.\end{eqnarray}$$

Take the curl to obtain

(11.99) $$\begin{eqnarray}-r\frac{\partial {\it\Omega}^{2}}{\partial z}\,\boldsymbol{e}_{{\it\phi}}=\boldsymbol{{\rm\nabla}}p\times \boldsymbol{{\rm\nabla}}\left(\frac{1}{{\it\rho}}\right)=\boldsymbol{{\rm\nabla}}T\times \boldsymbol{{\rm\nabla}}s.\end{eqnarray}$$

This is just the vorticity equation in a steady state. It is sometimes called the thermal wind equation. The equilibrium is called barotropic if $\boldsymbol{{\rm\nabla}}p$ is parallel to $\boldsymbol{{\rm\nabla}}{\it\rho}$ , otherwise it is called baroclinic. In a barotropic state the angular velocity is independent of $z$ : ${\it\Omega}={\it\Omega}(r)$ . This is a version of the Taylor–Proudman theorem ²⁰ which states that under certain conditions the velocity in a rotating fluid is independent of height.

We can also write

(11.100) $$\begin{eqnarray}\frac{1}{{\it\rho}}\boldsymbol{{\rm\nabla}}p=\boldsymbol{g}=-\boldsymbol{{\rm\nabla}}{\it\Phi}+r{\it\Omega}^{2}\,\boldsymbol{e}_{r},\end{eqnarray}$$

where $\boldsymbol{g}$ is the effective gravitational acceleration, including the centrifugal force associated with the (non-uniform) rotation.

In a barotropic state with ${\it\Omega}(r)$ we can write

(11.101) $$\begin{eqnarray}\boldsymbol{g}=-\boldsymbol{{\rm\nabla}}{\it\Phi}_{cg},\quad {\it\Phi}_{cg}={\it\Phi}(r,z)+{\it\Psi}(r),\quad {\it\Psi}=-\int r{\it\Omega}^{2}\,\text{d}r.\end{eqnarray}$$

Also, since $p=p({\it\rho})$ in the equilibrium state, we can define the pseudoenthalpy $\tilde{h}({\it\rho})$ such that $\text{d}\tilde{h}=\text{d}p/{\it\rho}$ . An example is a polytropic model for which

(11.102a,b ) $$\begin{eqnarray}p=K{\it\rho}^{1+1/m},\quad \tilde{h}=(m+1)K{\it\rho}^{1/m}.\end{eqnarray}$$

( $\tilde{h}$ equals the true enthalpy only if the equilibrium is homentropic.) The equilibrium condition then reduces to

(11.103) $$\begin{eqnarray}\mathbf{0}=-\boldsymbol{{\rm\nabla}}{\it\Phi}_{cg}-\boldsymbol{{\rm\nabla}}\tilde{h}\end{eqnarray}$$

or

(11.104) $$\begin{eqnarray}{\it\Phi}+{\it\Psi}+\tilde{h}=C=\text{const.}\end{eqnarray}$$

An example of a rapidly and differentially rotating equilibrium is an accretion disc around a central mass $M$ . For a non-self-gravitating disc ${\it\Phi}=-GM(r^{2}+z^{2})^{-1/2}$ . Assume the disc is barotropic and let the arbitrary additive constant in $\tilde{h}$ be defined (as in the polytropic example above) such that $\tilde{h}=0$ at the surfaces $z=\pm H(r)$ of the disc where ${\it\rho}=p=0$ . Then

(11.105) $$\begin{eqnarray}-GM(r^{2}+H^{2})^{-1/2}+{\it\Psi}(r)=C,\end{eqnarray}$$

from which

(11.106) $$\begin{eqnarray}r{\it\Omega}^{2}=-\frac{\text{d}}{\text{d}r}\left[GM(r^{2}+H^{2})^{-1/2}\right].\end{eqnarray}$$

For example, if $H={\it\epsilon}r$ with ${\it\epsilon}=\text{const.}$ being the aspect ratio of the disc, then

(11.107) $$\begin{eqnarray}{\it\Omega}^{2}=(1+{\it\epsilon}^{2})^{-1/2}\frac{GM}{r^{3}}.\end{eqnarray}$$

The thinner the disc is, the closer it is to Keplerian rotation. Once we have found the relation between ${\it\Omega}(r)$ and $H(r)$ , equation (11.104) then determines the spatial distribution of $\tilde{h}$ (and therefore of ${\it\rho}$ and $p$ ) within the disc.

11.6.2 Linear perturbations

The basic state is independent of $t$ and ${\it\phi}$ , allowing us to consider linear perturbations of the form

(11.108) $$\begin{eqnarray}\text{Re}[{\it\delta}u_{r}(r,z)\exp (-\text{i}{\it\omega}t+\text{i}m{\it\phi})],\quad \text{etc.},\end{eqnarray}$$

where $m$ is the azimuthal wavenumber (an integer). The linearized equations in the Cowling approximation are

(11.109) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle -\text{i}\hat{{\it\omega}}{\it\delta}u_{r}-2{\it\Omega}{\it\delta}u_{{\it\phi}}=-\frac{1}{{\it\rho}}\frac{\partial {\it\delta}p}{\partial r}+\frac{{\it\delta}{\it\rho}}{{\it\rho}^{2}}\frac{\partial p}{\partial r},\\ \displaystyle -\text{i}\hat{{\it\omega}}{\it\delta}u_{{\it\phi}}+\frac{1}{r}{\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(r^{2}{\it\Omega})=-\frac{\text{i}m\,{\it\delta}p}{{\it\rho}r},\\ \displaystyle -\text{i}\hat{{\it\omega}}{\it\delta}u_{z}=-\frac{1}{{\it\rho}}\frac{\partial {\it\delta}p}{\partial z}+\frac{{\it\delta}{\it\rho}}{{\it\rho}^{2}}\frac{\partial p}{\partial z},\\ \displaystyle -\text{i}\hat{{\it\omega}}{\it\delta}{\it\rho}+{\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}=-{\it\rho}\left[\frac{1}{r}\frac{\partial }{\partial r}(r\,{\it\delta}u_{r})+\frac{\text{i}m\,{\it\delta}u_{{\it\phi}}}{r}+\frac{\partial {\it\delta}u_{z}}{\partial z}\right],\\ \displaystyle -\text{i}\hat{{\it\omega}}{\it\delta}p+{\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p=-{\it\gamma}p\left[\frac{1}{r}\frac{\partial }{\partial r}(r\,{\it\delta}u_{r})+\frac{\text{i}m\,{\it\delta}u_{{\it\phi}}}{r}+\frac{\partial {\it\delta}u_{z}}{\partial z}\right],\end{array}\right\}\end{eqnarray}$$

where

(11.110) $$\begin{eqnarray}\hat{{\it\omega}}={\it\omega}-m{\it\Omega}\end{eqnarray}$$

is the intrinsic frequency, i.e. the angular frequency of the wave measured in a frame of reference that rotates with the local angular velocity of the fluid.

Eliminate ${\it\delta}u_{{\it\phi}}$ and ${\it\delta}{\it\rho}$ to obtain

(11.111) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (\hat{{\it\omega}}^{2}-A)\,{\it\delta}u_{r}-B\,{\it\delta}u_{z}=-\frac{\text{i}\hat{{\it\omega}}}{{\it\rho}}\left(\frac{\partial {\it\delta}p}{\partial r}-\frac{\partial p}{\partial r}\frac{{\it\delta}p}{{\it\gamma}p}\right)+2{\it\Omega}\frac{\text{i}m\,{\it\delta}p}{{\it\rho}r},\\ \displaystyle -C\,{\it\delta}u_{r}+(\hat{{\it\omega}}^{2}-D)\,{\it\delta}u_{z}=-\frac{\text{i}\hat{{\it\omega}}}{{\it\rho}}\left(\frac{\partial {\it\delta}p}{\partial z}-\frac{\partial p}{\partial z}\frac{{\it\delta}p}{{\it\gamma}p}\right),\end{array}\right\}\end{eqnarray}$$

where

(11.112) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle A=\frac{2{\it\Omega}}{r}\frac{\partial }{\partial r}(r^{2}{\it\Omega})-\frac{1}{{\it\rho}}\frac{\partial p}{\partial r}\left(\frac{1}{{\it\gamma}p}\frac{\partial p}{\partial r}-\frac{1}{{\it\rho}}\frac{\partial {\it\rho}}{\partial r}\right),\\ \displaystyle B=\frac{2{\it\Omega}}{r}\frac{\partial }{\partial z}(r^{2}{\it\Omega})-\frac{1}{{\it\rho}}\frac{\partial p}{\partial r}\left(\frac{1}{{\it\gamma}p}\frac{\partial p}{\partial z}-\frac{1}{{\it\rho}}\frac{\partial {\it\rho}}{\partial z}\right),\\ \displaystyle C=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial z}\left(\frac{1}{{\it\gamma}p}\frac{\partial p}{\partial r}-\frac{1}{{\it\rho}}\frac{\partial {\it\rho}}{\partial r}\right),\\ \displaystyle D=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial z}\left(\frac{1}{{\it\gamma}p}\frac{\partial p}{\partial z}-\frac{1}{{\it\rho}}\frac{\partial {\it\rho}}{\partial z}\right).\end{array}\right\}\end{eqnarray}$$

Note that $A$ , $B$ , $C$ and $D$ involve radial and vertical derivatives of the specific angular momentum $r^{2}{\it\Omega}$ and the specific entropy $s$ . The thermal wind equation implies

(11.113) $$\begin{eqnarray}B=C,\end{eqnarray}$$

so the matrix

(11.114) $$\begin{eqnarray}\unicode[STIX]{x1D648}=\left[\begin{array}{@{}cc@{}}A & B\\ C & D\end{array}\right]=\left[\begin{array}{@{}cc@{}}A & B\\ B & D\end{array}\right]\end{eqnarray}$$

is symmetric.

11.6.3 The Høiland criteria

It can be useful to introduce the Lagrangian displacement ${\bf\xi}$ such that

(11.115) $$\begin{eqnarray}{\rm\Delta}\boldsymbol{u}={\it\delta}\boldsymbol{u}+{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}=\frac{\text{D}{\bf\xi}}{\text{D}t},\end{eqnarray}$$

i.e.

(11.116a-c ) $$\begin{eqnarray}{\it\delta}u_{r}=-\text{i}\hat{{\it\omega}}{\it\xi}_{r},\quad {\it\delta}u_{{\it\phi}}=-\text{i}\hat{{\it\omega}}{\it\xi}_{{\it\phi}}-r{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Omega},\quad {\it\delta}u_{z}=-\text{i}\hat{{\it\omega}}{\it\xi}_{z}.\end{eqnarray}$$

Note that

(11.117) $$\begin{eqnarray}\frac{1}{r}\frac{\partial }{\partial r}(r\,{\it\delta}u_{r})+\frac{\text{i}m\,{\it\delta}u_{{\it\phi}}}{r}+\frac{\partial {\it\delta}u_{z}}{\partial z}=-\text{i}\hat{{\it\omega}}\left[\frac{1}{r}\frac{\partial }{\partial r}(r{\it\xi}_{r})+\frac{\text{i}m{\it\xi}_{{\it\phi}}}{r}+\frac{\partial {\it\xi}_{z}}{\partial z}\right].\end{eqnarray}$$

The linearized equations constitute an eigenvalue problem for ${\it\omega}$ but it is not self-adjoint except when $m=0$ . We specialize to the case $m=0$ (axisymmetric perturbations). Then

(11.118) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle ({\it\omega}^{2}-A){\it\xi}_{r}-B{\it\xi}_{z}=\frac{1}{{\it\rho}}\left(\frac{\partial {\it\delta}p}{\partial r}-\frac{\partial p}{\partial r}\frac{{\it\delta}p}{{\it\gamma}p}\right),\\ \displaystyle -B{\it\xi}_{r}+({\it\omega}^{2}-D){\it\xi}_{z}=\frac{1}{{\it\rho}}\left(\frac{\partial {\it\delta}p}{\partial z}-\frac{\partial p}{\partial z}\frac{{\it\delta}p}{{\it\gamma}p}\right),\end{array}\right\}\end{eqnarray}$$

with

(11.119) $$\begin{eqnarray}{\it\delta}p=-{\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}-{\bf\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p.\end{eqnarray}$$

Multiply the first of (11.118) by ${\it\rho}{\it\xi}_{r}^{\ast }$ and the second by ${\it\rho}{\it\xi}_{z}^{\ast }$ and integrate over the volume $V$ of the fluid (using the boundary condition ${\it\delta}p=0$ ) to obtain

(11.120) $$\begin{eqnarray}\displaystyle {\it\omega}^{2}\int _{V}{\it\rho}(|{\it\xi}_{r}|^{2}+|{\it\xi}_{z}|^{2})\,\text{d}V\! & = & \displaystyle \!\int _{V}\left[{\it\rho}Q({\bf\xi})+{\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\delta}p-\frac{{\it\delta}p}{{\it\gamma}p}{\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p\right]\,\text{d}V\nonumber\\ \displaystyle \! & = & \displaystyle \!\int _{V}\left[{\it\rho}Q({\bf\xi})-\frac{{\it\delta}p}{{\it\gamma}p}({\it\gamma}p\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}^{\ast }+{\bf\xi}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p)\right]\,\text{d}V\nonumber\\ \displaystyle \! & = & \displaystyle \!\int _{V}\left({\it\rho}Q({\bf\xi})+\frac{|{\it\delta}p|^{2}}{{\it\gamma}p}\right)\,\text{d}V,\end{eqnarray}$$

where

(11.121) $$\begin{eqnarray}Q({\bf\xi})=A|{\it\xi}_{r}|^{2}+B({\it\xi}_{r}^{\ast }{\it\xi}_{z}+{\it\xi}_{z}^{\ast }{\it\xi}_{r})+D|{\it\xi}_{z}|^{2}=\left[\begin{array}{@{}cc@{}}{\it\xi}_{r}^{\ast } & {\it\xi}_{z}^{\ast }\end{array}\right]\left[\begin{array}{@{}cc@{}}A & B\\ B & D\end{array}\right]\left[\begin{array}{@{}c@{}}{\it\xi}_{r}\\ {\it\xi}_{z}\end{array}\right]\end{eqnarray}$$

is the (real) Hermitian form associated with the matrix $\unicode[STIX]{x1D648}$ .

Note that this integral involves only the meridional components of the displacement. If we had not made the Cowling approximation there would be the usual negative definite contribution to ${\it\omega}^{2}$ from self-gravitation.

The above integral relation therefore shows that ${\it\omega}^{2}$ is real, and a variational property ensures that instability to axisymmetric perturbations occurs if and only if the integral on the right-hand side can be made negative by a suitable trial displacement. If $Q$ is positive definite then ${\it\omega}^{2}>0$ and we have stability. Now the characteristic equation of the matrix $\unicode[STIX]{x1D648}$ is

(11.122) $$\begin{eqnarray}{\it\lambda}^{2}-(A+D){\it\lambda}+AD-B^{2}=0.\end{eqnarray}$$

The eigenvalues ${\it\lambda}_{\pm }$ are both positive if and only if

(11.123a,b ) $$\begin{eqnarray}A+D>0\quad \text{and}\quad AD-B^{2}>0.\end{eqnarray}$$

If these conditions are satisfied throughout the fluid then $Q>0$ , which implies ${\it\omega}^{2}>0$ , so the fluid is stable to axisymmetric perturbations (neglecting self-gravitation). These conditions are also necessary for stability. If one of the eigenvalues is negative in some region of the meridional plane, then a trial displacement can be found which is localized in that region, has ${\it\delta}p=0$ and $Q<0$ , implying instability. (By choosing ${\bf\xi}$ in the correct direction and tuning $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\xi}$ appropriately, it is possible to arrange for ${\it\delta}p$ to vanish.)

Using $\ell =r^{2}{\it\Omega}$ (specific angular momentum) and $s=c_{p}({\it\gamma}^{-1}\ln p-\ln {\it\rho})+\text{const.}$ (specific entropy) for a perfect ideal gas, we have

(11.124) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \displaystyle A=\frac{1}{r^{3}}\frac{\partial \ell ^{2}}{\partial r}-\frac{g_{r}}{c_{p}}\frac{\partial s}{\partial r},\\ \displaystyle B=\frac{1}{r^{3}}\frac{\partial \ell ^{2}}{\partial z}-\frac{g_{r}}{c_{p}}\frac{\partial s}{\partial z}=-\frac{g_{z}}{c_{p}}\frac{\partial s}{\partial r},\\ \displaystyle D=-\frac{g_{z}}{c_{p}}\frac{\partial s}{\partial z},\end{array}\right\}\end{eqnarray}$$

so the two conditions become

(11.125) $$\begin{eqnarray}\frac{1}{r^{3}}\frac{\partial \ell ^{2}}{\partial r}-\frac{1}{c_{p}}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s>0\end{eqnarray}$$

and

(11.126) $$\begin{eqnarray}-g_{z}\left(\frac{\partial \ell ^{2}}{\partial r}\frac{\partial s}{\partial z}-\frac{\partial \ell ^{2}}{\partial z}\frac{\partial s}{\partial r}\right)>0.\end{eqnarray}$$

These are the Høiland stability criteria ²¹ .

(If the criteria are marginally satisfied a further investigation may be required.)

Consider first the non-rotating case $\ell =0$ . The first criterion reduces to the Schwarzschild criterion for convective stability,

(11.127) $$\begin{eqnarray}-\frac{1}{c_{p}}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s\equiv N^{2}>0.\end{eqnarray}$$

In the homentropic case $s=\text{const.}$ (which is a barotropic model) they reduce to the Rayleigh criterion ²² for centrifugal (inertial) stability,

(11.128) $$\begin{eqnarray}\frac{\text{d}\ell ^{2}}{\text{d}r}>0,\end{eqnarray}$$

which states that the specific angular momentum should increase with $r$ for stability.

The second Høiland criterion is equivalent to

(11.129) $$\begin{eqnarray}(\boldsymbol{e}_{r}\times (-\boldsymbol{g}))\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\ell ^{2}\times \boldsymbol{{\rm\nabla}}s)>0.\end{eqnarray}$$

In other words the vectors $\boldsymbol{e}_{r}\times (-\boldsymbol{g})$ and $\boldsymbol{{\rm\nabla}}\ell ^{2}\times \boldsymbol{{\rm\nabla}}s$ should be parallel (rather than antiparallel). In a rotating star in which the specific entropy increases outwards, for stability we require that the specific angular momentum should increase with $r$ on each surface of constant entropy.

Related example: A.26.