Skip to main content Accessibility help



  • Access
  • Cited by 13



MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Astrophysical fluid dynamics
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Astrophysical fluid dynamics
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Astrophysical fluid dynamics
        Available formats
Export citation


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

1 Introduction

1.1 Areas of application

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

Important areas of application include:

  1. (i) Instabilities in astrophysical fluids

  2. (ii) Convection

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

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

  5. (v) Astrophysical dynamos

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

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

  8. (viii) Tidal disruption and stellar collisions

  9. (ix) Supernovae

  10. (x) Planetary nebulae

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

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

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

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

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

  16. (xvi) Planetary atmospheric dynamics

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

  18. (xviii) Cosmology and structure formation

1.2 Theoretical varieties

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

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

  1. (i) HD: hydrodynamics

  2. (ii) MHD: magnetohydrodynamics

  3. (iii) RHD: radiation hydrodynamics

  4. (iv) RMHD: radiation magnetohydrodynamics

  5. (v) GRHD: general relativistic hydrodynamics

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

1.3 Characteristic features

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

2 Ideal gas dynamics

2.1 Fluid variables

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

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

2.2 Eulerian and Lagrangian viewpoints

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

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

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

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

2.3 Material points and structures

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

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

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

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

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

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

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

They therefore evolve according to

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

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

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

2.4 Equation of mass conservation

The equation of mass conservation,

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

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

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

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

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

2.5 Equation of motion

The equation of motion,

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

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

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

Viscous forces are neglected in ideal gas dynamics.

2.6 Poisson’s equation

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

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

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

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

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

A non-self-gravitating fluid is one of negligible mass for which ${\it\Phi}_{\mathit{int}}$ can be neglected. More generally, the Cowling approximation 1 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}$$


(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}$$


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

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

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

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

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

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

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

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

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

2.9.2 Derivation of fluid equations

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

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

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

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

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

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

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

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

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

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

2.9.3 Validity of a fluid approach

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

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

Some typical order-of-magnitude estimates:

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

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

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

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

3 Ideal magnetohydrodynamics

3.1 Elementary derivation of the MHD equations

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

3.1.1 Galilean electromagnetism

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

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

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

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

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

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

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

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

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

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

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

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

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

These relations correspond to the limit of the Lorentz transformation for electromagnetic fields 3 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}$$


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

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

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

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

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

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

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

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

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

The combination

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

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

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

is known as the plasma beta.

3.1.4 Self-consistency of approximations

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

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

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

3.1.5 Summary of the MHD equations

The full set of ideal MHD equations is

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

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

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

3.2 Physical interpretation of MHD

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

3.2.1 Kinematics of the magnetic field

The ideal induction equation

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

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

One way to show this result is to use the identity

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

to write the induction equation in the form

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

and use the equation of mass conservation,

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

to obtain

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

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

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

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

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

Precisely the same equation as the ideal induction equation,

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

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

Related examples: A.5, A.6.

3.2.2 The Lorentz force

The Lorentz force per unit volume,

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

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

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

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

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

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

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

showing the magnetic tension and pressure.

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

3.2.3 Differential rotation and torsional Alfvén waves

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

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

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

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

The ideal induction equation reduces to (exercise)

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

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

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

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

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.7A.9.

4 Conservation laws, symmetries and hyperbolic structure

4.1 Introduction

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

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

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

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

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

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

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

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

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

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

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

4.2 Synthesis of the total energy equation

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

Kinetic energy:

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

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

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

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

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

Sum of these three:

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

The last term can be rewritten as

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

Using mass conservation:

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

Magnetic energy:

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

Total energy:

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

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

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

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

The total energy equation is then

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

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

4.3 Other conservation laws in ideal MHD

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

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

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

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

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

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

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

In ideal MHD, therefore, magnetic helicity is conserved:

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

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

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

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

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

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

The cross-helicity in a volume $V$ is

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

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

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

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

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

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

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

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

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

Related examples: A.10, A.11.

4.4 Symmetries of the equations

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

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

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

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

  4. (iv) Galilean transformations.

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

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

    (4.28a,b ) $$\begin{eqnarray}\boldsymbol{x}\mapsto {\it\lambda}\,\boldsymbol{x},\quad t\mapsto {\it\mu}\,t,\end{eqnarray}$$
    (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}$$


(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}$$


(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}$$


(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}$$


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

is related to the inertia tensor of the system,

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

is a kinetic energy tensor and

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

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

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

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

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

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

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

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

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

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

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

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

5 Linear waves in homogeneous media

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

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

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

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

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

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

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

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

The linearized equation of motion is

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

where the perturbation of total pressure is

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

The gravitational potential perturbation satisfies the linearized Poisson equation

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

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

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


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

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

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

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

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

where we have changed the sign and omitted the tilde.

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

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

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

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

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

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

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

For Alfvén waves, therefore,

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

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

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

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

and then

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

These equations can be written in the form

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

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

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

which simplifies to

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

The two solutions

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

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

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

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

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

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

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

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

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

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

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

We can interpret:

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

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

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

Related example: A.12.

6 Nonlinear waves, shocks and other discontinuities

6.1 One-dimensional gas dynamics

6.1.1 Riemann’s analysis

The equations of mass conservation and motion in one dimension are

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

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

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


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

We add and subtract to obtain

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

Define the two Riemann invariants

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

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

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

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

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

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

6.1.2 Method of characteristics

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

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

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

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

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

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

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

Figure 6. Domains of dependence and of influence.

6.1.3 A simple wave

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

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

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

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

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

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

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

Figure 7. Formation of a shock from intersecting characteristics.

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

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

Indeed, the formal solution of the inviscid Burgers equation is

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

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

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

6.2 General analysis of simple nonlinear waves

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

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

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

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

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

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

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

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

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

the inviscid Burgers equation.

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

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

6.3 Shocks and other discontinuities

6.3.1 Jump conditions

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

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

Figure 9. A shock front in its rest frame.

Consider any equation in conservative form

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

For a stationary solution in one dimension,

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

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

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

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

From mass conservation:

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

From momentum conservation:

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

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

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

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

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

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

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

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

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

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

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

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

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

6.3.2 Non-magnetic shocks

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

(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}$$


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

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

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

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

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

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

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

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

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

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

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

6.3.3 Oblique shocks

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

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

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

6.3.4 Other discontinuities

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

6.3.5 MHD shocks and discontinuities

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

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

6.3.6 The Riemann problem

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

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

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

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

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

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