2  Magnetic fields

In MHD, we need to augment our fluid equations with a new variable, \(\Bb(\xb,t)\), called the magnetic field.

The S.I. unit of magnetic field is the Tesla (symbol \(\mathrm{T}\)). Typical orders of magnitude are \(10^{-3}\,\mathrm{T}\) for a fridge magnet, \(10^{-5}\,\mathrm{T}\) for the Earth’s magnetic field, or \(10^{-1}\,\mathrm{T}\) in a sunspot.

2.1 The solenoidal condition

Our first fundamental law is that \[ \boxed{\oint_S\Bb\cdot\dS = 0 \quad\textrm{for any closed surface $S$}.} \tag{2.1}\] This is the statement that no isolated magnetic monopoles can exist.

There are several consequences:

  1. If \(\Bb\) is differentiable, then it must satisfy Maxwell’s solenoidal condition \[ \nabla\cdot\Bb = 0. \tag{2.2}\]

  2. If \(\Bb\) is differentiable throughout some volume \(V\) in which Equation 2.1 is satisfied, then there exists a single-valued vector potential \(\Ab\) such that \(\Bb=\nabla\times\Ab\).

Proof of (i). Let \(V\) be the volume enclosed by \(S\) in Equation 2.1. Then, provided \(\Bb\) is differentiable in \(V\), the Divergence Theorem gives \(\int_V\nabla\cdot\Bb\dV = 0.\) Since Equation 2.1 holds for any \(S\), we must have \(\nabla\cdot\Bb=0\) in any region where \(\Bb\) is differentiable.

Proof of (ii). Proving this in full generality is outside the scope of this course. For \(V=\mathbb{R}^3\), a suitable \(\Ab\) is given by the Biot-Savart integral \[ \Ab(\xb) = \int_V\nabla G(\xb,\xb') \times \Bb(\xb')\dV', \quad\textrm{where } G(\xb,\xb') = \frac{1}{4\pi|\xb-\xb'|}. \tag{2.3}\] To see this, take the curl. For shorthand, write \(\Bb'\equiv\Bb(\xb')\). Then \[\begin{align} \nabla\times\Ab(\xb) &= \int_V\nabla\times\big(\nabla G\times\Bb'\big)\dV'\\ &= \int_V\nabla\times\nabla\times\big(G\Bb'\big)\dV' \quad \textrm{[since $\Bb'$ is independent of $\xb$]}\\ &= \int_V\Big( \nabla\big[\nabla\cdot(G\Bb')\big]- \Bb'\Delta G\Big)\dV' \quad \textrm{[using (F3) on Formula Sheet]}\\ &=\nabla\int_V\nabla\cdot(G\Bb')\dV' - \int_V\Delta G\Bb'\dV'\\ &= -\nabla\int_V\nabla'\cdot(G\Bb')\dV'- \int_V\Delta G\Bb'\dV' \quad \textrm{[using $\nabla'G=-\nabla G$ and $\nabla'\cdot\Bb'=0$]}\\ &= - \int_V\Delta G\Bb'\dV' \quad\textrm{[using $V=\mathbb{R}_3$]}. \end{align}\] In the last step we assumed that \(|\Bb'|\to 0\) sufficiently fast as \(|\xb'|\to\infty\), and applied the Divergence Theorem to the first integral. Now \(G\) is the Green’s function for the Poisson equation on \(\mathbb{R}^3\), so \[ \Delta G(\xb,\xb') = -\delta(\xb-\xb'), \tag{2.4}\] and hence \[ \nabla\times\Ab(\xb) = \int_V\Bb(\xb')\delta(\xb-\xb')\dV' = \Bb(\xb). \]

Computing \(\Ab\) from \(\Bb\) is analogous to determining a velocity field \(\ub\) from the vorticity \(\omb\).

There are more practical methods of computing \(\Ab\) (some on the problem sheet).

The vector potential is not unique, for if \(\Ab\) is a vector potential for some \(\Bb\), then so is \(\Ab'=\Ab + \nabla\chi\) for any sufficiently differentiable scalar function \(\chi\).

If \(\Bb\) is two dimensional, with no \(z\)-component and no dependence on \(z\), then we have the convenient representation \[ \Bb = \nabla\times\big(A(x,y)\eb_z\big). \tag{2.5}\] Then \(\nabla\cdot\Bb=0\) automatically and \[ \Bb = \nabla A\times\eb_z = \ddy{A}{y}\eb_x - \ddy{A}{x}\eb_y. \tag{2.6}\] The function \(A(x,y)\) is called a flux function because the difference \(A(x_1,y_1)-A(x_2,y_2)\) gives the magnetic flux (per unit height in \(z\)) between these two points.

We want \(\displaystyle\ddy{A}{y}=B_0y\), so \(\displaystyle A(x,y)=\int_0^yB_0y'\,\mathrm{d}y' = \frac{B_0y^2}{2} + f(x)\) for some function \(f(x)\).

Differentiating gives \(\displaystyle\ddy{A}{x} = f'(x)\), which we want equal to \(-B_0x\). Thus \(\displaystyle f(x)= -\frac{B_0x^2}{2} + c\), and \(\displaystyle A(x,y) = \frac{B_0}{2}(y^2 - x^2) + c\).

Consider now the magnetic flux per unit height, \(\Phi(x)\), through a vertical surface between \((0,0)\) and the point \((x,0)\).

The normal to this surface is \(\eb_y\), so from \(\Bb\) we calculate \[ \textrm{[flux per unit height]} = \int_0^1\int_0^x\eb_y\cdot\Bb(x',0,z')\,\mathrm{d}x'\mathrm{d}z' = B_0\int_0^xx'\,\mathrm{d}x' = \frac{B_0x^2}{2}. \] Indeed the flux function gives \(\displaystyle A(0,0)-A(x,0)= \frac{B_0x^2}{2}\).

A magnetic field line (sometimes called line of force) is a curve \(\xb(s)\) that is everywhere tangent to \(\Bb\), so satisfies \[ \dds{\xb} = \frac{\Bb(\xb(s))}{|\Bb(\xb(s))|}. \tag{2.7}\] In this parametrisation, \(\displaystyle\left|\dds{\xb}\right|=1\), so \(s\) represents arclength along the field line from \(\xb(0)\).

Magnetic field lines are analogous to streamlines of a velocity field.

Field lines are curves \(\xb(s) = x(s)\eb_x + y(s)\eb_y + z(s)\eb_z\) satisfying \[ \dds{x} = \frac{B_x(\xb(s))}{|\Bb(\xb(s))|}, \quad \dds{y} = \frac{B_y(\xb(s))}{|\Bb(\xb(s))|}, \quad \dds{z} = 0. \] The third equation shows that field lines lie in planes of constant \(z\). We can eliminate \(s\) by combining the first two equations, which gives \[ \frac{\mathrm{d}x}{\mathrm{d}y} = \frac{B_x(\xb(s))}{B_y(\xb(s))} = \frac{y}{x} \quad \implies \int x\,\mathrm{d}x = \int y\,\mathrm{d}y \quad \implies y^2-x^2 = \textrm{constant}. \]

Notice that the field lines in this example are curves of constant \(A\). This is generally true for magnetic fields of the form Equation 2.5 because \(\Bb\cdot\nabla A = 0\).

The solenoidal condition \(\nabla\cdot\Bb=0\) implies that a magnetic field line cannot end at a point (except where \(\Bb=\bfzero\)), but must do one of three things:

  1. form a closed curve;
  2. extend to infinity (as in the above example);
  3. stay within a finite region and fill space “ergodically” without closing.

The so-called Arnold-Beltrami-Childress field has the form \[ \Bb = \big(A\sin z + C\cos y\big)\eb_x + \big(B\sin x + A\cos z\big)\eb_y + \big(C\sin y + B\cos x\big)\eb_z, \] on the triply-periodic domain \(x,y,z\in [-\pi,\pi)\). Despite the simple expression for \(\Bb\), the field line pattern is very complex and contains both regular and ergodic regions.

The plot below shows numerically computed field lines for the case \(A=1\), \(B=\sqrt{2/3}\), \(C=\sqrt{1/3}\), found by solving Equation 2.7. The region \(x<0\) is largely “regular” while \(x>0\) is largely ergodic.

2.2 Ampère’s Law

Our second fundamental law is the relation between magnetic field and electric current, which takes the form of Ampère’s Law \[ \boxed{\oint_\gamma\Bb\cdot\dl = \mu_0\int_S\Jb\cdot\dS \quad \textrm{for any closed curve $\gamma$ spanned by surface $S$}.} \tag{2.8}\] Here \(\Jb(\xb,t)\) is the current density, and \(\mu_0\) is a physical constant called the vacuum permeability.

The magnetic field generated by the current in a wire.

Figure 2.1: Ampère’s Law is the general version of the basic physical observation that a current flowing through a wire generates an azimuthal magnetic field around it.

The S.I. unit of current density is Amps per square metre, \(\mathrm{A}\,\mathrm{m}^{-2}\). In S.I. units, the vacuum permeability is \(\mu_0\approx 4\pi\times 10^{-7}\,\mathrm{T}\mathrm{A}^{-1}\mathrm{m} = 4\pi\times 10^{-7}\,\mathrm{H}\mathrm{m}^{-1}\).

Since the problem is axisymmetric, we work in cylindrical coordinates \((r,\theta,z)\), and assume that \(\Bb=\Bb(r)\). For such a magnetic field, \(\nabla\cdot\Bb=0\) requires [formula sheet (F9)] that \[ \frac{1}{r}\ddy{}{r}(rB_r) = 0 \quad \iff \quad B_r = \frac{a}{r}. \] To avoid a singularity at \(r=0\), we must set \(a=0\) so \(B_r=0\).

Now we can use Equation 2.8 to calculate \(B_\theta\) and \(B_z\). For \(B_\theta\), take \(\gamma\) to be a circle \(r=\textrm{constant}\), so Equation 2.8 gives \[ \int_0^{2\pi}B_\theta(r)r\,\mathrm{d}\theta = \mu_0 I \quad \iff \quad B_\theta(r) = \frac{\mu_0I}{2\pi r}. \] For \(B_z\), take \(\gamma\) to be a rectangle in a plane \(\theta=\textrm{constant}\), with \(r\in[0,R]\) and \(z\in[0,1]\). Then Equation 2.8 with \(B_r=0\) gives \[ B_z(R) - B_z(0) = 0, \] from which we conclude that \(B_z=b\) (constant). So \[ \Bb = \frac{\mu_0I}{2\pi r}\eb_\theta + b\eb_z. \] The fact that \(\Bb\) is not uniquely defined is a consequence of the fact that it is an integral of \(\Jb\).

One could also compute \(\Bb\) from \(\Jb\) using the Biot-Savart integral, as when computing \(\Ab\) from \(\Bb\) (see problem sheet).

Note that Equation 2.8 neglects the “displacement current” which Maxwell famously added to model non-stationary situations. This is justified provided that the fluid motions of interest are much slower than the speed of light. Thus MHD is an entirely classical (non-relativistic) theory, and in particular electromagnetic waves have been filtered out.

Again, we can derive the differential form of Equation 2.8 in a region where \(\Bb\) and \(\Jb\) are differentiable. Applying Stokes’ Theorem shows that \(\oint_\gamma\Bb\cdot\dl = \int_S\nabla\times\Bb\cdot\dS\). Since Equation 2.8 holds for any closed curve \(\gamma\), we must have that \[ \nabla\times\Bb = \mu_0\Jb. \tag{2.9}\] It follows that \(\nabla\cdot\Jb=0\).

2.3 Potential fields

Notice that we can have a non-zero magnetic field even when \(\Jb=\bfzero\), providing that \(\nabla\times\Bb=\bfzero\). Such magnetic fields are called potential fields, because they can be written as \(\Bb=\nabla\Phi\) for some scalar potential \(\Phi(\xb)\).

Example – a uniform magnetic field \(\Bb=b\eb_z\) (from the previous example).

Clearly \(\nabla\times\Bb=\bfzero\) so this magnetic field is potential, with \(\Phi = bz + c\) for an arbitrary constant \(c\). It extends to infinity.

Notice that \(\nabla\cdot\Bb=0\) implies that \(\nabla\cdot\nabla\Phi=0\), so \[ \Delta\Phi = 0 \quad \textrm{(Laplace's equation)}. \tag{2.10}\]

Physically, such a uniform magnetic field with no generating current is not possible, and (as we will see) would have infinite magnetic energy. However, it is often used as a local model, for example in the theory of MHD waves…

If the domain is not the whole of \(\mathbb{R}^3\), we can have non-trivial solutions to Equation 2.10.

You should already be familiar with solution of Laplace’s equation by separation of variables. Let \(\Bb=\nabla\Phi\) and try the ansatz \(\Phi(x,z) = X(x)Y(y)\). Then \[ \Delta\Phi = 0 \qquad \iff \qquad \frac{X''}{X} + \frac{Y''}{Y} = 0 \qquad \iff \qquad \frac{X''}{X} = - \frac{Y''}{Y} = -k^2, \] where \(-k^2\) is a constant. I choose it to be negative so that the solutions are trigonometric functions that can satisfy the boundary conditions \(B_x=\partial\Phi/\partial x\) at \(x=0,1\). So \(X(x) = a\cos(k x)\) with \(k=n\pi\) for \(n\in\mathbb{Z}\). The \(y\)-equation then gives \[ \frac{Y''}{Y} = (n\pi)^2. \] Solutions are growing or decaying exponentials, but the boundary condition as \(y\to\infty\) implies that \(Y(y) = b\mathrm{e}^{-n\pi y}\). So combining the constants \(a\) and \(b\) for each term, the general solution is \[ \Phi(x,y) = \sum_{n=0}^{\infty}c_n\cos(n\pi x)\mathrm{e}^{-n\pi y} \] with corresponding magnetic field \[\begin{align} B_x &= \frac{\partial\Phi}{\partial x} = -\sum_{n=1}^\infty n\pi c_n\sin(n\pi x) \mathrm{e}^{-n\pi y},\\ B_y &= \frac{\partial\Phi}{\partial y} = -\sum_{n=1}^\infty n\pi c_n\cos(n\pi x) \mathrm{e}^{-n\pi y}. \end{align}\]

To specify a particular solution, one would need to determine the coefficients \(c_n\), for example by specifying the boundary condition \(B_y(x,0)\) and using the usual Fourier series approach.

For the simple case \(c_k = -\delta_{1k}\), we have \(\Bb = \pi\big(\sin(\pi x)\eb_x + \cos(\pi x)\eb_y\big)\mathrm{e}^{-\pi y}\) and you can show that the field lines are the curves \(\sin(\pi x)\mathrm{e}^{-\pi y} = \textrm{const.}\) (check this!) The plot shows this magnetic field (where the domain is unbounded in \(z\)).

In fact, the potential field in a bounded region is uniquely determined by the boundary conditions.

Proposition 2.1 For a volume \(V\) with boundary condition \(\nb\cdot\Bb=g(\xb)\) on \(\partial V\) there is at most one solution \(\Bb=\nabla\Phi\) with \(\Delta\Phi=0\) in \(V\).

Proof. Suppose potential fields \(\Bb_1=\nabla\Phi_1\) and \(\Bb_2=\nabla\Phi_2\) both satisfy the boundary condition, \(\nb\cdot\Bb_1|_{\partial V}=\nb\cdot\Bb_2|_{\partial V}=0\). Define \(\Bb_{12}=\Bb_1 - \Bb_2\), so that \(\nb\cdot\Bb_{12}|_{\partial V}=0\).

Now \[\begin{align} \int_V|\Bb_{12}|^2\dV &= \int_V\nabla(\Phi_1-\Phi_2)\cdot\Bb_{12}\dV\\ &= \int_V \Big[\nabla\cdot\Big((\Phi_1-\Phi_2)\Bb_{12}\Big) - (\Phi_1-\Phi_2)\nabla\cdot\Bb_{12}\Big]\dV \quad \textrm{[by formula sheet (F1)]}\\ &= \oint_{\partial V}(\Phi_1-\Phi_2)\,\nb\cdot\Bb_{12}\,\mathrm{d}S \quad \textrm{[since $\nabla\cdot\Bb_{12}=0$]}\\ &= 0 \quad \textrm{[using $\nb\cdot\Bb_{12}|_{\partial V}=0$]}. \end{align}\] Hence \(\Bb_{12}\equiv\bfzero\) so \(\Bb_1\equiv\Bb_2\).

The proof also works in the case of a semi-infinite region like the previous example, provided \(|\Bb|\to 0\) as \(|\xb|\to\infty\) so that the final integral still vanishes.

A simple model for the magnetic field above the Earth’s surface \(r=1\) is to assume \(\Jb=\bfzero\) for \(r>1\) and specify \(B_r(1,\theta,\phi)=2\cos\theta\) on the boundary \(r=1\). We assume also that \(|\Bb|\to 0\) as \(r\to\infty\). This gives a boundary value problem for \(\Bb=\nabla\Phi\) in \(r>1\).

In spherical coordinates, \(\Delta\Phi=0\) is given by formula (F17). Assuming no \(\phi\)-dependence (for simplicity), \[ \frac{1}{r}\frac{\partial^2}{\partial r^2}(r\Phi) + \frac{1}{r^2\sin\theta}\ddy{}{\theta}\left(\sin\theta\ddy{\Phi}{\theta}\right) = 0. \] Again we can solve by separation of variables. The ansatz \(\Phi(r,\theta)=R(r)\Theta(\theta)\) gives \[ \frac{r}{R}\frac{\mathrm{d}^2}{\mathrm{d}r^2}(rR) = -\frac{1}{\Theta\sin\theta}\frac{\mathrm{d}}{\mathrm{d}\theta}\left(\sin\theta\frac{\mathrm{d}\Theta}{\mathrm{d}\theta}\right) = \alpha. \] where \(\alpha\) is constant. Taking first the \(\theta\) equation, it is simplest to define \(s=\cos\theta\) so that \(\displaystyle \frac{\mathrm{d}}{\mathrm{d}\theta} = -\sin\theta\frac{\mathrm{d}}{\mathrm{d}s}\) and the equation becomes \[ -\dds{}\left((1-s^2)\dds{\Theta}\right) = \alpha\Theta \quad \iff \quad (1-s^2)\frac{\mathrm{d}^2\Theta}{\mathrm{d}s^2} -2s\dds{\Theta} + \alpha\Theta = 0, \] which is the Legendre equation [equation (F21) on the formula sheet]. We need the solution to be defined at \(s=\pm 1\), since these are the north and south poles \(\theta=0,\pi\). So we must take \(\alpha=\ell(\ell+1)\) for \(\ell\in\mathbb{Z}\), so that the solutions are Legendre polynomials \(\Theta(s) = P_\ell(s)\).

The \(r\)-equation has the form \[ r\frac{\mathrm{d}}{\mathrm{d}r}\left(r\frac{\mathrm{d}R}{\mathrm{d}r} + R\right) - \alpha R = 0 \quad \iff \quad r^2\frac{\mathrm{d}^2R}{\mathrm{d}r^2} + 2r\frac{\mathrm{d}R}{\mathrm{d}r} - \ell(\ell+1)R = 0. \] Trying \(R=r^\beta\) gives \[\begin{align} \beta(\beta-1)r^\beta + 2\beta r^\beta - \ell(\ell+1) r^\beta = 0 \quad &\iff \beta^2 + \beta - \ell(\ell+1) = 0\\ &\iff (\beta - \ell)(\beta + \ell + 1) = 0. \end{align}\] So the solution for given \(\ell\) has the form \(R(r) = cr^\ell + dr^{-\ell-1}\).

The general solution is then \[ \Phi(r,\theta) = \sum_{l=0}^\infty(c_\ell r^\ell + d_\ell r^{-\ell-1})P_\ell(\cos\theta). \]

Now we apply the boundary conditions. For the solution to decay as \(r\to\infty\), we need \(c_\ell=0\) for all \(\ell\). We then determine the \(d_\ell\) using the lower boundary condition. We need \[ B_r(1,\theta,\phi) = \left.\frac{\partial\Phi}{\partial r}\right|_{r=1} = -\sum_{\ell=0}^\infty (\ell+1)d_\ell P_\ell(\cos\theta) = 2\cos\theta. \tag{2.11}\] Since \(\cos\theta = P_1(\cos\theta)\), we see that \[ d_\ell=\begin{cases} -1 & \ell=1,\\ 0 & \textrm{otherwise}. \end{cases} \] So the solution in \(r>1\) is \[ \Bb = \frac{2\cos\theta}{r^3}\evr + \frac{\sin\theta}{r^3}\evt. \]

This is a magnetic dipole with its axis aligned with the north pole:

…we can determine the coefficients \(d_\ell\) similar to Fourier series, using the orthogonality relation \[ \int_{-1}^1P_\ell(x)P_{\ell'}(x)\,\mathrm{d}x = \frac{2}{2\ell + 1}\delta_{\ell\ell'}. \] [This corresponds to the particular normalization where \(P_\ell(1)=1\) for all \(\ell\).]

In our example, multiplying Equation 2.11 by \(P_1(\cos\theta)\) and integrating over \(\cos\theta\) gives \[ -2d_1\frac{2}{3} = \frac{2}{3}(2) \quad \implies d_1 = -1, \quad \textrm{with $d_\ell=0$ for $\ell\neq 1$.} \]

Separation of variables also works if we allow \(\Phi\) to depend on all three coordinates \(r\), \(\theta\), \(\phi\). Such solutions are used routinely to model the magnetic field in the Sun’s atmosphere.

A potential field extrapolation from real solar data.

Figure 2.2: A potential field extrapolation using boundary data from the NASA Solar Dynamics Observatory satellite. The black-white shading on the Sun’s surface shows the measured \(B_r\) (determined by telescope observations through splitting of spectral lines by the Zeeman effect), where the strong bipolar patches are called active regions. You can find an up-to-date model at http://suntoday.lmsal.com.

We will see later that potential fields are equilibria of the full MHD equations, and play an important role because they are the minimum energy solutions.