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 approximationFootnote ¹ 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.

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 equationsFootnote ² 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 fieldsFootnote ³ 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 wavesFootnote ⁴ (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 isorotationFootnote ⁵ .

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.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 Reference Chandrasekhar1969).

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 equationFootnote ⁷ 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 relationsFootnote ⁸

(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.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.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.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 (Reference Parker1958)Footnote ¹³ 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 (Reference Bondi1952)Footnote ¹⁴ 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.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 fieldFootnote ¹⁵ . 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 assistanceFootnote ¹⁶ . 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.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.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ä frequencyFootnote ¹⁷ (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 criterionFootnote ¹⁸ 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)